Functions_MatVectComplex.cxx
1 // Copyright (C) 2001-2011 Marc Duruflé
2 // Copyright (C) 2001-2011 Vivien Mallet
3 //
4 // This file is part of the linear-algebra library Seldon,
5 // http://seldon.sourceforge.net/.
6 //
7 // Seldon is free software; you can redistribute it and/or modify it under the
8 // terms of the GNU Lesser General Public License as published by the Free
9 // Software Foundation; either version 2.1 of the License, or (at your option)
10 // any later version.
11 //
12 // Seldon is distributed in the hope that it will be useful, but WITHOUT ANY
13 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15 // more details.
16 //
17 // You should have received a copy of the GNU Lesser General Public License
18 // along with Seldon. If not, see http://www.gnu.org/licenses/.
19 
20 
21 #ifndef SELDON_FILE_FUNCTIONS_MATVECT_COMPLEX_CXX
22 
23 /*
24  Functions defined in this file:
25  (storage RowComplexSparse, ArrayRowComplexSparse, etc)
26 
27  alpha.M*X + beta.Y -> Y
28  MltAdd(alpha, M, X, beta, Y)
29 
30  SOR(A, X, B, omega, nb_iter, type_ssor)
31 
32 */
33 
34 namespace Seldon
35 {
36 
37  /*************
38  * MltVector *
39  *************/
40 
41 
42  template <class T1, class Prop1, class Allocator1,
43  class T2, class Storage2, class Allocator2,
44  class T4, class Storage4, class Allocator4>
45  void MltVector(const Matrix<T1, Prop1, RowComplexSparse, Allocator1>& M,
46  const Vector<T2, Storage2, Allocator2>& X,
47  Vector<T4, Storage4, Allocator4>& Y)
48  {
49  int i; long j;
50 
51  int ma = M.GetM();
52 
53 #ifdef SELDON_CHECK_DIMENSIONS
54  CheckDim(M, X, Y, "MltAdd(alpha, M, X, beta, Y)");
55 #endif
56 
57  typename ClassComplexType<T1>::Treal rzero(0);
58  T1 zero(0, 0);
59  T1 temp;
60 
61  long* real_ptr = M.GetRealPtr();
62  long* imag_ptr = M.GetImagPtr();
63  int* real_ind = M.GetRealInd();
64  int* imag_ind = M.GetImagInd();
65  typename Matrix<T1, Prop1, RowComplexSparse, Allocator1>::pointer
66  real_data = M.GetRealData();
67  typename Matrix<T1, Prop1, RowComplexSparse, Allocator1>::pointer
68  imag_data = M.GetImagData();
69 
70  for (i = 0; i < ma; i++)
71  {
72  temp = zero;
73  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
74  temp += real_data[j] * X(real_ind[j]);
75 
76  Y(i) = temp;
77  }
78 
79  for (i = 0; i < ma; i++)
80  {
81  temp = zero;
82  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
83  temp += T1(rzero, imag_data[j]) * X(imag_ind[j]);
84 
85  Y(i) += temp;
86  }
87  }
88 
89 
90  template <class T1, class Prop1, class Allocator1,
91  class T2, class Storage2, class Allocator2,
92  class T4, class Storage4, class Allocator4>
93  void MltVector(const Matrix<T1, Prop1, ColComplexSparse, Allocator1>& M,
94  const Vector<T2, Storage2, Allocator2>& X,
95  Vector<T4, Storage4, Allocator4>& Y)
96  {
97  int i; long j;
98 
99 #ifdef SELDON_CHECK_DIMENSIONS
100  CheckDim(M, X, Y, "MltAdd(alpha, M, X, beta, Y)");
101 #endif
102 
103  typename ClassComplexType<T1>::Treal rzero(0);
104  long* real_ptr = M.GetRealPtr();
105  long* imag_ptr = M.GetImagPtr();
106  int* real_ind = M.GetRealInd();
107  int* imag_ind = M.GetImagInd();
108  typename Matrix<T1, Prop1, ColComplexSparse, Allocator1>::pointer
109  real_data = M.GetRealData();
110  typename Matrix<T1, Prop1, ColComplexSparse, Allocator1>::pointer
111  imag_data = M.GetImagData();
112 
113  Y.Fill(0);
114  for (i = 0; i < M.GetN(); i++)
115  {
116  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
117  Y(real_ind[j]) += real_data[j] * X(i);
118 
119  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
120  Y(imag_ind[j]) += T1(rzero, imag_data[j]) * X(i);
121  }
122  }
123 
124 
125  /*** Symmetric complex sparse matrices ***/
126 
127 
128  template <class T1, class Prop1, class Allocator1,
129  class T2, class Storage2, class Allocator2,
130  class T4, class Storage4, class Allocator4>
131  void MltVector(const Matrix<T1, Prop1, RowSymComplexSparse, Allocator1>& M,
132  const Vector<T2, Storage2, Allocator2>& X,
133  Vector<T4, Storage4, Allocator4>& Y)
134  {
135  int ma = M.GetM();
136 
137 #ifdef SELDON_CHECK_DIMENSIONS
138  CheckDim(M, X, Y, "MltAdd(alpha, M, X, beta, Y)");
139 #endif
140 
141  int i; long j;
142  typename ClassComplexType<T1>::Treal rzero(0);
143  T1 zero(0, 0);
144  T1 temp;
145 
146  long* real_ptr = M.GetRealPtr();
147  long* imag_ptr = M.GetImagPtr();
148  int* real_ind = M.GetRealInd();
149  int* imag_ind = M.GetImagInd();
150  typename Matrix<T1, Prop1, RowSymComplexSparse, Allocator1>::pointer
151  real_data = M.GetRealData();
152  typename Matrix<T1, Prop1, RowSymComplexSparse, Allocator1>::pointer
153  imag_data = M.GetImagData();
154 
155  Y.Fill(0);
156  for (i = 0; i < ma; i++)
157  {
158  temp = zero;
159  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
160  {
161  temp += real_data[j] * X(real_ind[j]);
162  if (real_ind[j] != i)
163  Y(real_ind[j]) += real_data[j] * X(i);
164  }
165 
166  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
167  {
168  temp += T1(rzero, imag_data[j]) * X(imag_ind[j]);
169  if (imag_ind[j] != i)
170  Y(imag_ind[j]) += T1(rzero, imag_data[j]) * X(i);
171  }
172 
173  Y(i) += temp;
174  }
175  }
176 
177 
178  template <class T1, class Prop1, class Allocator1,
179  class T2, class Storage2, class Allocator2,
180  class T4, class Storage4, class Allocator4>
181  void MltVector(const Matrix<T1, Prop1, ColSymComplexSparse, Allocator1>& M,
182  const Vector<T2, Storage2, Allocator2>& X,
183  Vector<T4, Storage4, Allocator4>& Y)
184  {
185  int ma = M.GetM();
186 
187 #ifdef SELDON_CHECK_DIMENSIONS
188  CheckDim(M, X, Y, "MltAdd(alpha, M, X, beta, Y)");
189 #endif
190 
191  int i; long j;
192  typename ClassComplexType<T1>::Treal rzero(0);
193  T1 zero(0, 0);
194  T1 temp;
195 
196  long* real_ptr = M.GetRealPtr();
197  long* imag_ptr = M.GetImagPtr();
198  int* real_ind = M.GetRealInd();
199  int* imag_ind = M.GetImagInd();
200  typename Matrix<T1, Prop1, ColSymComplexSparse, Allocator1>::pointer
201  real_data = M.GetRealData();
202  typename Matrix<T1, Prop1, ColSymComplexSparse, Allocator1>::pointer
203  imag_data = M.GetImagData();
204 
205  Y.Fill(0);
206  for (i = 0; i<ma; i++)
207  {
208  temp = zero;
209  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
210  {
211  temp += real_data[j] * X(real_ind[j]);
212  if (real_ind[j] != i)
213  Y(real_ind[j]) += real_data[j] * X(i);
214  }
215 
216  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
217  {
218  temp += T1(rzero, imag_data[j]) * X(imag_ind[j]);
219  if (imag_ind[j] != i)
220  Y(imag_ind[j]) += T1(rzero, imag_data[j]) * X(i);
221  }
222 
223  Y(i) += temp;
224  }
225  }
226 
227 
228  /*** Complex sparse matrices, *Trans ***/
229 
230 
231  template <class T1, class Prop1, class Allocator1,
232  class T2, class Storage2, class Allocator2,
233  class T4, class Storage4, class Allocator4>
234  void MltVector(const SeldonTranspose& Trans,
235  const Matrix<T1, Prop1, RowComplexSparse, Allocator1>& M,
236  const Vector<T2, Storage2, Allocator2>& X,
237  Vector<T4, Storage4, Allocator4>& Y)
238  {
239  if (Trans.NoTrans())
240  {
241  MltVector(M, X, Y);
242  return;
243  }
244 
245  int i; long j;
246 
247  int ma = M.GetM();
248 
249 #ifdef SELDON_CHECK_DIMENSIONS
250  CheckDim(Trans, M, X, Y, "MltAdd(alpha, SeldonTrans, M, X, beta, Y)");
251 #endif
252 
253  typename ClassComplexType<T1>::Treal rzero(0);
254  long* real_ptr = M.GetRealPtr();
255  long* imag_ptr = M.GetImagPtr();
256  int* real_ind = M.GetRealInd();
257  int* imag_ind = M.GetImagInd();
258  typename Matrix<T1, Prop1, RowComplexSparse, Allocator1>::pointer
259  real_data = M.GetRealData();
260  typename Matrix<T1, Prop1, RowComplexSparse, Allocator1>::pointer
261  imag_data = M.GetImagData();
262 
263  Y.Fill(0);
264 
265  if (Trans.Trans())
266  {
267  for (i = 0; i < ma; i++)
268  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
269  Y(real_ind[j]) += real_data[j] * X(i);
270 
271  for (i = 0; i < ma; i++)
272  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
273  Y(imag_ind[j]) += T1(rzero, imag_data[j]) * X(i);
274  }
275  else
276  {
277  for (i = 0; i < ma; i++)
278  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
279  Y(real_ind[j]) += real_data[j] * X(i);
280 
281  for (i = 0; i < ma; i++)
282  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
283  Y(imag_ind[j]) += T1(rzero, - imag_data[j]) * X(i);
284  }
285  }
286 
287 
288  template <class T1, class Prop1, class Allocator1,
289  class T2, class Storage2, class Allocator2,
290  class T4, class Storage4, class Allocator4>
291  void MltVector(const SeldonTranspose& Trans,
292  const Matrix<T1, Prop1, ColComplexSparse, Allocator1>& M,
293  const Vector<T2, Storage2, Allocator2>& X,
294  Vector<T4, Storage4, Allocator4>& Y)
295  {
296  if (Trans.NoTrans())
297  {
298  MltVector(M, X, Y);
299  return;
300  }
301 
302  int i; long j;
303 
304 #ifdef SELDON_CHECK_DIMENSIONS
305  CheckDim(Trans, M, X, Y, "MltAdd(alpha, SeldonTrans, M, X, beta, Y)");
306 #endif
307 
308  T4 temp;
309  typename ClassComplexType<T1>::Treal rzero(0);
310 
311  long* real_ptr = M.GetRealPtr();
312  long* imag_ptr = M.GetImagPtr();
313  int* real_ind = M.GetRealInd();
314  int* imag_ind = M.GetImagInd();
315  typename Matrix<T1, Prop1, ColComplexSparse, Allocator1>::pointer
316  real_data = M.GetRealData();
317  typename Matrix<T1, Prop1, ColComplexSparse, Allocator1>::pointer
318  imag_data = M.GetImagData();
319 
320  if (Trans.Trans())
321  {
322  for (i = 0; i < M.GetN(); i++)
323  {
324  SetComplexZero(temp);
325  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
326  temp += real_data[j] * X(real_ind[j]);
327 
328  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
329  temp += T1(rzero, imag_data[j]) * X(imag_ind[j]);
330 
331  Y(i) = temp;
332  }
333  }
334  else
335  {
336  for (i = 0; i < M.GetN(); i++)
337  {
338  SetComplexZero(temp);
339  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
340  temp += real_data[j] * X(real_ind[j]);
341 
342  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
343  temp += T1(rzero, -imag_data[j]) * X(imag_ind[j]);
344 
345  Y(i) = temp;
346  }
347  }
348  }
349 
350 
351  /*** Symmetric complex sparse matrices, *Trans ***/
352 
353 
354  template <class T1, class Prop1, class Allocator1,
355  class T2, class Storage2, class Allocator2,
356  class T4, class Storage4, class Allocator4>
357  void MltVector(const SeldonTranspose& Trans,
358  const Matrix<T1, Prop1, RowSymComplexSparse, Allocator1>& M,
359  const Vector<T2, Storage2, Allocator2>& X,
360  Vector<T4, Storage4, Allocator4>& Y)
361  {
362  if (!Trans.ConjTrans())
363  {
364  MltVector(M, X, Y);
365  return;
366  }
367 
368  int ma = M.GetM();
369 
370 #ifdef SELDON_CHECK_DIMENSIONS
371  CheckDim(Trans, M, X, Y, "MltAdd(alpha, SeldonConjTrans, M, X, beta, Y)");
372 #endif
373 
374  int i; long j;
375  T1 zero(0, 0);
376  T1 temp;
377  typename ClassComplexType<T1>::Treal rzero(0);
378 
379  long* real_ptr = M.GetRealPtr();
380  long* imag_ptr = M.GetImagPtr();
381  int* real_ind = M.GetRealInd();
382  int* imag_ind = M.GetImagInd();
383  typename Matrix<T1, Prop1, RowSymComplexSparse, Allocator1>::pointer
384  real_data = M.GetRealData();
385  typename Matrix<T1, Prop1, RowSymComplexSparse, Allocator1>::pointer
386  imag_data = M.GetImagData();
387 
388  Y.Fill(0);
389  for (i = 0; i < ma; i++)
390  {
391  temp = zero;
392  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
393  {
394  temp += real_data[j] * X(real_ind[j]);
395  if (real_ind[j] != i)
396  Y(real_ind[j]) += real_data[j] * X(i);
397  }
398 
399  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
400  {
401  temp += T1(rzero, - imag_data[j]) * X(imag_ind[j]);
402  if (imag_ind[j] != i)
403  Y(imag_ind[j]) += T1(rzero, - imag_data[j]) * X(i);
404  }
405 
406  Y(i) += temp;
407  }
408  }
409 
410 
411  template <class T1, class Prop1, class Allocator1,
412  class T2, class Storage2, class Allocator2,
413  class T4, class Storage4, class Allocator4>
414  void MltVector(const SeldonTranspose& Trans,
415  const Matrix<T1, Prop1, ColSymComplexSparse, Allocator1>& M,
416  const Vector<T2, Storage2, Allocator2>& X,
417  Vector<T4, Storage4, Allocator4>& Y)
418  {
419  if (!Trans.ConjTrans())
420  {
421  MltVector(M, X, Y);
422  return;
423  }
424 
425  int ma = M.GetM();
426 
427 #ifdef SELDON_CHECK_DIMENSIONS
428  CheckDim(Trans, M, X, Y, "MltAdd(alpha, SeldonConjTrans, M, X, beta, Y)");
429 #endif
430 
431  int i; long j;
432  T1 zero(0, 0);
433  T1 temp;
434  typename ClassComplexType<T1>::Treal rzero(0);
435 
436  long* real_ptr = M.GetRealPtr();
437  long* imag_ptr = M.GetImagPtr();
438  int* real_ind = M.GetRealInd();
439  int* imag_ind = M.GetImagInd();
440  typename Matrix<T1, Prop1, ColSymComplexSparse, Allocator1>::pointer
441  real_data = M.GetRealData();
442  typename Matrix<T1, Prop1, ColSymComplexSparse, Allocator1>::pointer
443  imag_data = M.GetImagData();
444 
445  Y.Fill(0);
446  for (i = 0; i < ma; i++)
447  {
448  temp = zero;
449  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
450  {
451  temp += real_data[j] * X(real_ind[j]);
452  if (real_ind[j] != i)
453  Y(real_ind[j]) += real_data[j] * X(i);
454  }
455 
456  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
457  {
458  temp += T1(rzero, - imag_data[j]) * X(imag_ind[j]);
459  if (imag_ind[j] != i)
460  Y(imag_ind[j]) += T1(rzero, - imag_data[j]) * X(i);
461  }
462 
463  Y(i) += temp;
464  }
465  }
466 
467 
468  /*** ArrayRowSymComplexSparse ***/
469 
470 
471  template<class T1, class T2, class T3,
472  class Allocator1, class Allocator2, class Allocator3>
473  void MltVector(const Matrix<T1, Symmetric,
474  ArrayRowSymComplexSparse, Allocator1>& A,
475  const Vector<T2, VectFull, Allocator2>& B,
476  Vector<T3, VectFull, Allocator3>& C)
477  {
478  C.Fill(0);
479 
480  int m = A.GetM(), n, p;
481  typename ClassComplexType<T1>::Treal val;
482  for (int i = 0 ; i < m ; i++)
483  {
484  n = A.GetRealRowSize(i);
485  for (int k = 0; k < n ; k++)
486  {
487  p = A.IndexReal(i, k);
488  val = A.ValueReal(i, k);
489  if (p == i)
490  C(i) += val * B(i);
491  else
492  {
493  C(i) += val * B(p);
494  C(p) += val * B(i);
495  }
496  }
497 
498  n = A.GetImagRowSize(i);
499  for (int k = 0; k < n ; k++)
500  {
501  p = A.IndexImag(i, k);
502  val = A.ValueImag(i, k);
503  if (p == i)
504  C(i) += T1(0, val) * B(i);
505  else
506  {
507  C(i) += T1(0, val) * B(p);
508  C(p) += T1(0, val) * B(i);
509  }
510  }
511  }
512  }
513 
514 
515  template<class T1, class T2, class T3,
516  class Allocator1, class Allocator2, class Allocator3>
517  void MltVector(const SeldonTranspose& Trans, const Matrix<T1, Symmetric,
518  ArrayRowSymComplexSparse, Allocator1>& A,
519  const Vector<T2, VectFull, Allocator2>& B,
520  Vector<T3, VectFull, Allocator3>& C)
521  {
522  if (!Trans.ConjTrans())
523  {
524  MltVector(A, B, C);
525  return;
526  }
527 
528  C.Fill(0);
529  int m = A.GetM(),n,p;
530  typename ClassComplexType<T1>::Treal val;
531  for (int i = 0 ; i < m ; i++)
532  {
533  n = A.GetRealRowSize(i);
534  for (int k = 0; k < n ; k++)
535  {
536  p = A.IndexReal(i, k);
537  val = A.ValueReal(i, k);
538  if (p == i)
539  C(i) += val * B(i);
540  else
541  {
542  C(i) += val * B(p);
543  C(p) += val * B(i);
544  }
545  }
546  n = A.GetImagRowSize(i);
547  for (int k = 0; k < n ; k++)
548  {
549  p = A.IndexImag(i, k);
550  val = A.ValueImag(i, k);
551  if (p == i)
552  C(i) -= T1(0, val) * B(i);
553  else
554  {
555  C(i) -= T1(0, val) * B(p);
556  C(p) -= T1(0, val) * B(i);
557  }
558  }
559  }
560  }
561 
562 
563  template<class T1, class T2, class T3,
564  class Allocator1, class Allocator2, class Allocator3>
565  void MltVector(const Matrix<T1, Symmetric,
566  ArrayColSymComplexSparse, Allocator1>& A,
567  const Vector<T2, VectFull, Allocator2>& B,
568  Vector<T3, VectFull, Allocator3>& C)
569  {
570  C.Fill(0);
571  int m = A.GetM(), n, p;
572  typename ClassComplexType<T1>::Treal val;
573  for (int i = 0 ; i < m ; i++)
574  {
575  n = A.GetRealColumnSize(i);
576  for (int k = 0; k < n ; k++)
577  {
578  p = A.IndexReal(i, k);
579  val = A.ValueReal(i, k);
580  if (p == i)
581  C(i) += val * B(i);
582  else
583  {
584  C(i) += val * B(p);
585  C(p) += val * B(i);
586  }
587  }
588 
589  n = A.GetImagColumnSize(i);
590  for (int k = 0; k < n ; k++)
591  {
592  p = A.IndexImag(i, k);
593  val = A.ValueImag(i, k);
594  if (p == i)
595  C(i) += T1(0, val) * B(i);
596  else
597  {
598  C(i) += T1(0, val) * B(p);
599  C(p) += T1(0, val) * B(i);
600  }
601  }
602  }
603  }
604 
605 
606  template<class T1, class T2, class T3,
607  class Allocator1, class Allocator2, class Allocator3>
608  void MltVector(const SeldonTranspose& Trans, const Matrix<T1, Symmetric,
609  ArrayColSymComplexSparse, Allocator1>& A,
610  const Vector<T2, VectFull, Allocator2>& B,
611  Vector<T3, VectFull, Allocator3>& C)
612  {
613  if (!Trans.ConjTrans())
614  {
615  MltVector(A, B, C);
616  return;
617  }
618 
619  C.Fill(0);
620 
621  int m = A.GetM(), n, p;
622  typename ClassComplexType<T1>::Treal val;
623  for (int i = 0 ; i < m ; i++)
624  {
625  n = A.GetRealColumnSize(i);
626  for (int k = 0; k < n ; k++)
627  {
628  p = A.IndexReal(i, k);
629  val = A.ValueReal(i, k);
630  if (p == i)
631  C(i) += val * B(i);
632  else
633  {
634  C(i) += val * B(p);
635  C(p) += val * B(i);
636  }
637  }
638 
639  n = A.GetImagColumnSize(i);
640  for (int k = 0; k < n ; k++)
641  {
642  p = A.IndexImag(i, k);
643  val = A.ValueImag(i, k);
644  if (p == i)
645  C(i) -= T1(0, val) * B(i);
646  else
647  {
648  C(i) -= T1(0, val) * B(p);
649  C(p) -= T1(0, val) * B(i);
650  }
651  }
652  }
653  }
654 
655 
656  /*** ArrayRowComplexSparse ***/
657 
658 
659  template<class T1, class T2, class T3,
660  class Allocator1, class Allocator2, class Allocator3>
661  void MltVector(const Matrix<T1, General,
662  ArrayRowComplexSparse, Allocator1>& A,
663  const Vector<T2, VectFull, Allocator2>& B,
664  Vector<T3, VectFull, Allocator3>& C)
665  {
666  int m = A.GetM(), n;
667  T3 temp;
668  for (int i = 0 ; i < m ; i++)
669  {
670  SetComplexZero(temp);
671  n = A.GetRealRowSize(i);
672  for (int k = 0; k < n ; k++)
673  temp += A.ValueReal(i, k)*B(A.IndexReal(i, k));
674 
675  n = A.GetImagRowSize(i);
676  for (int k = 0; k < n ; k++)
677  temp += T1(0, A.ValueImag(i, k)) * B(A.IndexImag(i, k));
678 
679  C(i) = temp;
680  }
681  }
682 
683 
684  template<class T1, class T2, class T3,
685  class Allocator1, class Allocator2, class Allocator3>
686  void MltVector(const SeldonTranspose& Trans, const Matrix<T1, General,
687  ArrayRowComplexSparse, Allocator1>& A,
688  const Vector<T2, VectFull, Allocator2>& B,
689  Vector<T3, VectFull, Allocator3>& C)
690  {
691  if (Trans.NoTrans())
692  {
693  MltVector(A, B, C);
694  return;
695  }
696 
697  C.Fill(0);
698 
699  int m = A.GetM(), n, p;
700  typename ClassComplexType<T1>::Treal val;
701 
702  if (Trans.Trans())
703  {
704  for (int i = 0 ; i < m ; i++)
705  {
706  n = A.GetRealRowSize(i);
707  for (int k = 0; k < n ; k++)
708  {
709  p = A.IndexReal(i, k);
710  val = A.ValueReal(i, k);
711  C(p) += val * B(i);
712  }
713 
714  n = A.GetImagRowSize(i);
715  for (int k = 0; k < n ; k++)
716  {
717  p = A.IndexImag(i, k);
718  val = A.ValueImag(i, k);
719  C(p) += T1(0, val) * B(i);
720  }
721  }
722  }
723  else
724  {
725  for (int i = 0 ; i < m ; i++)
726  {
727  n = A.GetRealRowSize(i);
728  for (int k = 0; k < n ; k++)
729  {
730  p = A.IndexReal(i, k);
731  val = A.ValueReal(i, k);
732  C(p) += val * B(i);
733  }
734 
735  n = A.GetImagRowSize(i);
736  for (int k = 0; k < n ; k++)
737  {
738  p = A.IndexImag(i, k);
739  val = A.ValueImag(i, k);
740  C(p) -= T1(0, val) * B(i);
741  }
742  }
743  }
744  }
745 
746 
747  template<class T1, class T2, class T3,
748  class Allocator1, class Allocator2, class Allocator3>
749  void MltVector(const Matrix<T1, General,
750  ArrayColComplexSparse, Allocator1>& A,
751  const Vector<T2, VectFull, Allocator2>& B,
752  Vector<T3, VectFull, Allocator3>& C)
753  {
754  C.Fill(0);
755 
756  int n, p;
757  typename ClassComplexType<T1>::Treal val;
758  for (int i = 0 ; i < A.GetN(); i++)
759  {
760  n = A.GetRealColumnSize(i);
761  for (int k = 0; k < n ; k++)
762  {
763  p = A.IndexReal(i, k);
764  val = A.ValueReal(i, k);
765  C(p) += val * B(i);
766  }
767 
768  n = A.GetImagColumnSize(i);
769  for (int k = 0; k < n ; k++)
770  {
771  p = A.IndexImag(i, k);
772  val = A.ValueImag(i, k);
773  C(p) += T1(0, val) * B(i);
774  }
775  }
776  }
777 
778 
779  template<class T1, class T2, class T3,
780  class Allocator1, class Allocator2, class Allocator3>
781  void MltVector(const SeldonTranspose& Trans, const Matrix<T1, General,
782  ArrayColComplexSparse, Allocator1>& A,
783  const Vector<T2, VectFull, Allocator2>& B,
784  Vector<T3, VectFull, Allocator3>& C)
785  {
786  if (Trans.NoTrans())
787  {
788  MltVector(A, B, C);
789  return;
790  }
791 
792  T3 temp; int n;
793 
794  if (Trans.Trans())
795  {
796  for (int i = 0; i < A.GetN(); i++)
797  {
798  SetComplexZero(temp);
799  n = A.GetRealColumnSize(i);
800  for (int k = 0; k < n ; k++)
801  temp += A.ValueReal(i, k) * B(A.IndexReal(i, k));
802 
803  n = A.GetImagColumnSize(i);
804  for (int k = 0; k < n ; k++)
805  temp += T1(0, A.ValueImag(i, k)) * B(A.IndexImag(i, k));
806 
807  C(i) = temp;
808  }
809  }
810  else
811  {
812  for (int i = 0; i < A.GetN(); i++)
813  {
814  SetComplexZero(temp);
815  n = A.GetRealColumnSize(i);
816  for (int k = 0; k < n ; k++)
817  temp += A.ValueReal(i, k) * B(A.IndexReal(i, k));
818 
819  n = A.GetImagColumnSize(i);
820  for (int k = 0; k < n ; k++)
821  temp -= T1(0, A.ValueImag(i, k)) * B(A.IndexImag(i, k));
822 
823  C(i) = temp;
824  }
825  }
826  }
827 
828 
829  /****************
830  * MltAddVector *
831  ****************/
832 
833 
834  /*** Complex sparse matrices ***/
835 
836 
837  template <class T0,
838  class T1, class Prop1, class Allocator1,
839  class T2, class Storage2, class Allocator2,
840  class T3,
841  class T4, class Storage4, class Allocator4>
842  void MltAddVector(const T0& alpha,
843  const Matrix<T1, Prop1, RowComplexSparse, Allocator1>& M,
844  const Vector<T2, Storage2, Allocator2>& X,
845  const T3& beta, Vector<T4, Storage4, Allocator4>& Y)
846  {
847  int i; long j;
848 
849  int ma = M.GetM();
850 
851 #ifdef SELDON_CHECK_DIMENSIONS
852  CheckDim(M, X, Y, "MltAdd(alpha, M, X, beta, Y)");
853 #endif
854 
855  Mlt(beta, Y);
856 
857  typename ClassComplexType<T1>::Treal rzero(0);
858  T1 zero(0, 0);
859  T1 temp;
860 
861  long* real_ptr = M.GetRealPtr();
862  long* imag_ptr = M.GetImagPtr();
863  int* real_ind = M.GetRealInd();
864  int* imag_ind = M.GetImagInd();
865  typename Matrix<T1, Prop1, RowComplexSparse, Allocator1>::pointer
866  real_data = M.GetRealData();
867  typename Matrix<T1, Prop1, RowComplexSparse, Allocator1>::pointer
868  imag_data = M.GetImagData();
869 
870  for (i = 0; i < ma; i++)
871  {
872  temp = zero;
873  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
874  temp += real_data[j] * X(real_ind[j]);
875  Y(i) += alpha * temp;
876  }
877 
878  for (i = 0; i < ma; i++)
879  {
880  temp = zero;
881  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
882  temp += T1(rzero, imag_data[j]) * X(imag_ind[j]);
883  Y(i) += alpha * temp;
884  }
885  }
886 
887 
888  template <class T0,
889  class T1, class Prop1, class Allocator1,
890  class T2, class Storage2, class Allocator2,
891  class T3,
892  class T4, class Storage4, class Allocator4>
893  void MltAddVector(const T0& alpha,
894  const Matrix<T1, Prop1, ColComplexSparse, Allocator1>& M,
895  const Vector<T2, Storage2, Allocator2>& X,
896  const T3& beta, Vector<T4, Storage4, Allocator4>& Y)
897  {
898  int i; long j;
899 
900 #ifdef SELDON_CHECK_DIMENSIONS
901  CheckDim(M, X, Y, "MltAdd(alpha, M, X, beta, Y)");
902 #endif
903 
904  Mlt(beta, Y);
905 
906  typename ClassComplexType<T1>::Treal rzero(0);
907  long* real_ptr = M.GetRealPtr();
908  long* imag_ptr = M.GetImagPtr();
909  int* real_ind = M.GetRealInd();
910  int* imag_ind = M.GetImagInd();
911  typename Matrix<T1, Prop1, ColComplexSparse, Allocator1>::pointer
912  real_data = M.GetRealData();
913  typename Matrix<T1, Prop1, ColComplexSparse, Allocator1>::pointer
914  imag_data = M.GetImagData();
915 
916  for (i = 0; i < M.GetN(); i++)
917  {
918  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
919  Y(real_ind[j]) += alpha * real_data[j] * X(i);
920 
921  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
922  Y(imag_ind[j]) += alpha * T1(rzero, imag_data[j]) * X(i);
923  }
924  }
925 
926 
927  /*** Symmetric complex sparse matrices ***/
928 
929 
930  template <class T0,
931  class T1, class Prop1, class Allocator1,
932  class T2, class Storage2, class Allocator2,
933  class T3,
934  class T4, class Storage4, class Allocator4>
935  void MltAddVector(const T0& alpha,
936  const Matrix<T1, Prop1, RowSymComplexSparse, Allocator1>& M,
937  const Vector<T2, Storage2, Allocator2>& X,
938  const T3& beta, Vector<T4, Storage4, Allocator4>& Y)
939  {
940  int ma = M.GetM();
941 
942 #ifdef SELDON_CHECK_DIMENSIONS
943  CheckDim(M, X, Y, "MltAdd(alpha, M, X, beta, Y)");
944 #endif
945 
946  Mlt(beta, Y);
947 
948  int i; long j;
949  typename ClassComplexType<T1>::Treal rzero(0);
950  T1 zero(0, 0);
951  T1 temp;
952 
953  long* real_ptr = M.GetRealPtr();
954  long* imag_ptr = M.GetImagPtr();
955  int* real_ind = M.GetRealInd();
956  int* imag_ind = M.GetImagInd();
957  typename Matrix<T1, Prop1, RowSymComplexSparse, Allocator1>::pointer
958  real_data = M.GetRealData();
959  typename Matrix<T1, Prop1, RowSymComplexSparse, Allocator1>::pointer
960  imag_data = M.GetImagData();
961 
962  for (i = 0; i < ma; i++)
963  {
964  temp = zero;
965  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
966  {
967  temp += real_data[j] * X(real_ind[j]);
968  if (real_ind[j] != i)
969  Y(real_ind[j]) += alpha * real_data[j] * X(i);
970  }
971 
972  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
973  {
974  temp += T1(rzero, imag_data[j]) * X(imag_ind[j]);
975  if (imag_ind[j] != i)
976  Y(imag_ind[j]) += alpha * T1(rzero, imag_data[j]) * X(i);
977  }
978 
979  Y(i) += alpha * temp;
980  }
981  }
982 
983 
984  template <class T0,
985  class T1, class Prop1, class Allocator1,
986  class T2, class Storage2, class Allocator2,
987  class T3,
988  class T4, class Storage4, class Allocator4>
989  void MltAddVector(const T0& alpha,
990  const Matrix<T1, Prop1, ColSymComplexSparse, Allocator1>& M,
991  const Vector<T2, Storage2, Allocator2>& X,
992  const T3& beta, Vector<T4, Storage4, Allocator4>& Y)
993  {
994  int ma = M.GetM();
995 
996 #ifdef SELDON_CHECK_DIMENSIONS
997  CheckDim(M, X, Y, "MltAdd(alpha, M, X, beta, Y)");
998 #endif
999 
1000  Mlt(beta, Y);
1001 
1002  int i; long j;
1003  typename ClassComplexType<T1>::Treal rzero(0);
1004  T1 zero(0, 0);
1005  T1 temp;
1006 
1007  long* real_ptr = M.GetRealPtr();
1008  long* imag_ptr = M.GetImagPtr();
1009  int* real_ind = M.GetRealInd();
1010  int* imag_ind = M.GetImagInd();
1011  typename Matrix<T1, Prop1, ColSymComplexSparse, Allocator1>::pointer
1012  real_data = M.GetRealData();
1013  typename Matrix<T1, Prop1, ColSymComplexSparse, Allocator1>::pointer
1014  imag_data = M.GetImagData();
1015 
1016  for (i = 0; i<ma; i++)
1017  {
1018  temp = zero;
1019  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
1020  {
1021  temp += real_data[j] * X(real_ind[j]);
1022  if (real_ind[j] != i)
1023  Y(real_ind[j]) += alpha * real_data[j] * X(i);
1024  }
1025 
1026  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
1027  {
1028  temp += T1(rzero, imag_data[j]) * X(imag_ind[j]);
1029  if (imag_ind[j] != i)
1030  Y(imag_ind[j]) += alpha * T1(rzero, imag_data[j]) * X(i);
1031  }
1032 
1033  Y(i) += alpha * temp;
1034  }
1035  }
1036 
1037 
1038  /*** Complex sparse matrices, *Trans ***/
1039 
1040 
1041  template <class T0,
1042  class T1, class Prop1, class Allocator1,
1043  class T2, class Storage2, class Allocator2,
1044  class T3,
1045  class T4, class Storage4, class Allocator4>
1046  void MltAddVector(const T0& alpha,
1047  const SeldonTranspose& Trans,
1048  const Matrix<T1, Prop1, RowComplexSparse, Allocator1>& M,
1049  const Vector<T2, Storage2, Allocator2>& X,
1050  const T3& beta, Vector<T4, Storage4, Allocator4>& Y)
1051  {
1052  if (Trans.NoTrans())
1053  {
1054  MltAddVector(alpha, M, X, beta, Y);
1055  return;
1056  }
1057 
1058  int i; long j;
1059 
1060  int ma = M.GetM();
1061 
1062 #ifdef SELDON_CHECK_DIMENSIONS
1063  CheckDim(Trans, M, X, Y, "MltAdd(alpha, SeldonTrans, M, X, beta, Y)");
1064 #endif
1065 
1066  Mlt(beta, Y);
1067 
1068  typename ClassComplexType<T1>::Treal rzero(0);
1069  long* real_ptr = M.GetRealPtr();
1070  long* imag_ptr = M.GetImagPtr();
1071  int* real_ind = M.GetRealInd();
1072  int* imag_ind = M.GetImagInd();
1073  typename Matrix<T1, Prop1, RowComplexSparse, Allocator1>::pointer
1074  real_data = M.GetRealData();
1075  typename Matrix<T1, Prop1, RowComplexSparse, Allocator1>::pointer
1076  imag_data = M.GetImagData();
1077 
1078  if (Trans.Trans())
1079  {
1080  for (i = 0; i < ma; i++)
1081  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
1082  Y(real_ind[j]) += alpha * real_data[j] * X(i);
1083 
1084  for (i = 0; i < ma; i++)
1085  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
1086  Y(imag_ind[j]) += alpha * T1(rzero, imag_data[j]) * X(i);
1087  }
1088  else
1089  {
1090  for (i = 0; i < ma; i++)
1091  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
1092  Y(real_ind[j]) += alpha * real_data[j] * X(i);
1093 
1094  for (i = 0; i < ma; i++)
1095  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
1096  Y(imag_ind[j]) += alpha * T1(rzero, - imag_data[j]) * X(i);
1097  }
1098  }
1099 
1100 
1101  template <class T0,
1102  class T1, class Prop1, class Allocator1,
1103  class T2, class Storage2, class Allocator2,
1104  class T3,
1105  class T4, class Storage4, class Allocator4>
1106  void MltAddVector(const T0& alpha,
1107  const SeldonTranspose& Trans,
1108  const Matrix<T1, Prop1, ColComplexSparse, Allocator1>& M,
1109  const Vector<T2, Storage2, Allocator2>& X,
1110  const T3& beta, Vector<T4, Storage4, Allocator4>& Y)
1111  {
1112  if (Trans.NoTrans())
1113  {
1114  MltAddVector(alpha, M, X, beta, Y);
1115  return;
1116  }
1117 
1118  int i; long j;
1119 
1120 #ifdef SELDON_CHECK_DIMENSIONS
1121  CheckDim(Trans, M, X, Y, "MltAdd(alpha, SeldonTrans, M, X, beta, Y)");
1122 #endif
1123 
1124  Mlt(beta, Y);
1125 
1126  T4 temp;
1127  typename ClassComplexType<T1>::Treal rzero(0);
1128 
1129  long* real_ptr = M.GetRealPtr();
1130  long* imag_ptr = M.GetImagPtr();
1131  int* real_ind = M.GetRealInd();
1132  int* imag_ind = M.GetImagInd();
1133  typename Matrix<T1, Prop1, ColComplexSparse, Allocator1>::pointer
1134  real_data = M.GetRealData();
1135  typename Matrix<T1, Prop1, ColComplexSparse, Allocator1>::pointer
1136  imag_data = M.GetImagData();
1137 
1138  if (Trans.Trans())
1139  {
1140  for (i = 0; i < M.GetN(); i++)
1141  {
1142  SetComplexZero(temp);
1143  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
1144  temp += real_data[j] * X(real_ind[j]);
1145 
1146  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
1147  temp += T1(rzero, imag_data[j]) * X(imag_ind[j]);
1148 
1149  Y(i) += alpha * temp;
1150  }
1151  }
1152  else
1153  {
1154  for (i = 0; i < M.GetN(); i++)
1155  {
1156  SetComplexZero(temp);
1157  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
1158  temp += real_data[j] * X(real_ind[j]);
1159 
1160  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
1161  temp += T1(rzero, -imag_data[j]) * X(imag_ind[j]);
1162 
1163  Y(i) += alpha * temp;
1164  }
1165  }
1166  }
1167 
1168 
1169  /*** Symmetric complex sparse matrices, *Trans ***/
1170 
1171 
1172  template <class T0,
1173  class T1, class Prop1, class Allocator1,
1174  class T2, class Storage2, class Allocator2,
1175  class T3,
1176  class T4, class Storage4, class Allocator4>
1177  void MltAddVector(const T0& alpha,
1178  const SeldonTranspose& Trans,
1179  const Matrix<T1, Prop1, RowSymComplexSparse, Allocator1>& M,
1180  const Vector<T2, Storage2, Allocator2>& X,
1181  const T3& beta, Vector<T4, Storage4, Allocator4>& Y)
1182  {
1183  if (!Trans.ConjTrans())
1184  {
1185  MltAddVector(alpha, M, X, beta, Y);
1186  return;
1187  }
1188 
1189  int ma = M.GetM();
1190 
1191 #ifdef SELDON_CHECK_DIMENSIONS
1192  CheckDim(Trans, M, X, Y, "MltAdd(alpha, SeldonConjTrans, M, X, beta, Y)");
1193 #endif
1194 
1195  Mlt(beta, Y);
1196 
1197  int i; long j;
1198  T1 zero(0, 0);
1199  T1 temp;
1200  typename ClassComplexType<T1>::Treal rzero(0);
1201 
1202  long* real_ptr = M.GetRealPtr();
1203  long* imag_ptr = M.GetImagPtr();
1204  int* real_ind = M.GetRealInd();
1205  int* imag_ind = M.GetImagInd();
1206  typename Matrix<T1, Prop1, RowSymComplexSparse, Allocator1>::pointer
1207  real_data = M.GetRealData();
1208  typename Matrix<T1, Prop1, RowSymComplexSparse, Allocator1>::pointer
1209  imag_data = M.GetImagData();
1210 
1211  for (i = 0; i < ma; i++)
1212  {
1213  temp = zero;
1214  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
1215  {
1216  temp += real_data[j] * X(real_ind[j]);
1217  if (real_ind[j] != i)
1218  Y(real_ind[j]) += alpha * real_data[j] * X(i);
1219  }
1220 
1221  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
1222  {
1223  temp += T1(rzero, - imag_data[j]) * X(imag_ind[j]);
1224  if (imag_ind[j] != i)
1225  Y(imag_ind[j]) += alpha * T1(rzero, - imag_data[j]) * X(i);
1226  }
1227 
1228  Y(i) += alpha * temp;
1229  }
1230  }
1231 
1232 
1233  template <class T0,
1234  class T1, class Prop1, class Allocator1,
1235  class T2, class Storage2, class Allocator2,
1236  class T3,
1237  class T4, class Storage4, class Allocator4>
1238  void MltAddVector(const T0& alpha,
1239  const SeldonTranspose& Trans,
1240  const Matrix<T1, Prop1, ColSymComplexSparse, Allocator1>& M,
1241  const Vector<T2, Storage2, Allocator2>& X,
1242  const T3& beta, Vector<T4, Storage4, Allocator4>& Y)
1243  {
1244  if (!Trans.ConjTrans())
1245  {
1246  MltAddVector(alpha, M, X, beta, Y);
1247  return;
1248  }
1249 
1250  int ma = M.GetM();
1251 
1252 #ifdef SELDON_CHECK_DIMENSIONS
1253  CheckDim(Trans, M, X, Y, "MltAdd(alpha, SeldonConjTrans, M, X, beta, Y)");
1254 #endif
1255 
1256  Mlt(beta, Y);
1257 
1258  int i; long j;
1259  T1 zero(0, 0);
1260  T1 temp;
1261  typename ClassComplexType<T1>::Treal rzero(0);
1262 
1263  long* real_ptr = M.GetRealPtr();
1264  long* imag_ptr = M.GetImagPtr();
1265  int* real_ind = M.GetRealInd();
1266  int* imag_ind = M.GetImagInd();
1267  typename Matrix<T1, Prop1, ColSymComplexSparse, Allocator1>::pointer
1268  real_data = M.GetRealData();
1269  typename Matrix<T1, Prop1, ColSymComplexSparse, Allocator1>::pointer
1270  imag_data = M.GetImagData();
1271 
1272  for (i = 0; i < ma; i++)
1273  {
1274  temp = zero;
1275  for (j = real_ptr[i]; j < real_ptr[i + 1]; j++)
1276  {
1277  temp += real_data[j] * X(real_ind[j]);
1278  if (real_ind[j] != i)
1279  Y(real_ind[j]) += alpha * real_data[j] * X(i);
1280  }
1281 
1282  for (j = imag_ptr[i]; j < imag_ptr[i + 1]; j++)
1283  {
1284  temp += T1(rzero, - imag_data[j]) * X(imag_ind[j]);
1285  if (imag_ind[j] != i)
1286  Y(imag_ind[j]) += alpha * T1(rzero, - imag_data[j]) * X(i);
1287  }
1288 
1289  Y(i) += alpha * temp;
1290  }
1291  }
1292 
1293 
1294  /*** ArrayRowSymComplexSparse ***/
1295 
1296 
1297  template<class T0, class T1, class T2, class T3, class T4,
1298  class Allocator1, class Allocator2, class Allocator3>
1299  void MltAddVector(const T0& alpha, const Matrix<T1, Symmetric,
1300  ArrayRowSymComplexSparse, Allocator1>& A,
1301  const Vector<T2, VectFull, Allocator2>& B,
1302  const T4& beta,
1303  Vector<T3, VectFull, Allocator3>& C)
1304  {
1305  T4 zero; T0 one;
1306  SetComplexZero(zero);
1307  SetComplexOne(one);
1308 
1309  if (beta == zero)
1310  C.Fill(0);
1311  else
1312  Mlt(beta, C);
1313 
1314  int m = A.GetM(), n, p;
1315  typename ClassComplexType<T1>::Treal val;
1316  T3 val_cplx;
1317  if (alpha == one)
1318  {
1319  for (int i = 0 ; i < m ; i++)
1320  {
1321  n = A.GetRealRowSize(i);
1322  for (int k = 0; k < n ; k++)
1323  {
1324  p = A.IndexReal(i, k);
1325  val = A.ValueReal(i, k);
1326  if (p == i)
1327  C(i) += val * B(i);
1328  else
1329  {
1330  C(i) += val * B(p);
1331  C(p) += val * B(i);
1332  }
1333  }
1334  n = A.GetImagRowSize(i);
1335  for (int k = 0; k < n ; k++)
1336  {
1337  p = A.IndexImag(i, k);
1338  val = A.ValueImag(i, k);
1339  if (p == i)
1340  C(i) += T1(0, val) * B(i);
1341  else
1342  {
1343  C(i) += T1(0, val) * B(p);
1344  C(p) += T1(0, val) * B(i);
1345  }
1346  }
1347  }
1348  }
1349  else // alpha != 1.
1350  {
1351  for (int i = 0 ; i < m ; i++)
1352  {
1353  n = A.GetRealRowSize(i);
1354  for (int k = 0; k < n ; k++)
1355  {
1356  p = A.IndexReal(i, k);
1357  val_cplx = alpha * T1(A.ValueReal(i, k), 0);
1358  if (p == i)
1359  C(i) += val_cplx * B(i);
1360  else
1361  {
1362  C(i) += val_cplx * B(p);
1363  C(p) += val_cplx * B(i);
1364  }
1365  }
1366  n = A.GetImagRowSize(i);
1367  for (int k = 0; k < n ; k++)
1368  {
1369  p = A.IndexImag(i, k);
1370  val_cplx = alpha * T1(0, A.ValueImag(i, k));
1371  if (p == i)
1372  C(i) += val_cplx * B(i);
1373  else
1374  {
1375  C(i) += val_cplx * B(p);
1376  C(p) += val_cplx * B(i);
1377  }
1378  }
1379  }
1380  }
1381  }
1382 
1383 
1384  template<class T0, class T1, class T2, class T3, class T4,
1385  class Allocator1, class Allocator2, class Allocator3>
1386  void MltAddVector(const T0& alpha,
1387  const SeldonTranspose& Trans, const Matrix<T1, Symmetric,
1388  ArrayRowSymComplexSparse, Allocator1>& A,
1389  const Vector<T2, VectFull, Allocator2>& B,
1390  const T4& beta,
1391  Vector<T3, VectFull, Allocator3>& C)
1392  {
1393  if (!Trans.ConjTrans())
1394  {
1395  MltAddVector(alpha, A, B, beta, C);
1396  return;
1397  }
1398 
1399  T4 zero; T0 one;
1400  SetComplexZero(zero);
1401  SetComplexOne(one);
1402 
1403  if (beta == zero)
1404  C.Fill(0);
1405  else
1406  Mlt(beta, C);
1407 
1408  int m = A.GetM(),n,p;
1409  typename ClassComplexType<T1>::Treal val;
1410  T3 val_cplx;
1411  if (alpha == one)
1412  {
1413  for (int i = 0 ; i < m ; i++)
1414  {
1415  n = A.GetRealRowSize(i);
1416  for (int k = 0; k < n ; k++)
1417  {
1418  p = A.IndexReal(i, k);
1419  val = A.ValueReal(i, k);
1420  if (p == i)
1421  C(i) += val * B(i);
1422  else
1423  {
1424  C(i) += val * B(p);
1425  C(p) += val * B(i);
1426  }
1427  }
1428  n = A.GetImagRowSize(i);
1429  for (int k = 0; k < n ; k++)
1430  {
1431  p = A.IndexImag(i, k);
1432  val = A.ValueImag(i, k);
1433  if (p == i)
1434  C(i) -= T1(0, val) * B(i);
1435  else
1436  {
1437  C(i) -= T1(0, val) * B(p);
1438  C(p) -= T1(0, val) * B(i);
1439  }
1440  }
1441  }
1442  }
1443  else
1444  {
1445  // alpha different from 1
1446  for (int i = 0 ; i < m ; i++)
1447  {
1448  n = A.GetRealRowSize(i);
1449  for (int k = 0; k < n ; k++)
1450  {
1451  p = A.IndexReal(i, k);
1452  val_cplx = alpha * T1(A.ValueReal(i, k), 0);
1453  if (p == i)
1454  C(i) += val_cplx * B(i);
1455  else
1456  {
1457  C(i) += val_cplx * B(p);
1458  C(p) += val_cplx * B(i);
1459  }
1460  }
1461  n = A.GetImagRowSize(i);
1462  for (int k = 0; k < n ; k++)
1463  {
1464  p = A.IndexImag(i, k);
1465  val_cplx = alpha * T1(0, A.ValueImag(i, k));
1466  if (p == i)
1467  C(i) -= val_cplx * B(i);
1468  else
1469  {
1470  C(i) -= val_cplx * B(p);
1471  C(p) -= val_cplx * B(i);
1472  }
1473  }
1474  }
1475  }
1476  }
1477 
1478 
1479  template<class T0, class T1, class T2, class T3, class T4,
1480  class Allocator1, class Allocator2, class Allocator3>
1481  void MltAddVector(const T0& alpha, const Matrix<T1, Symmetric,
1482  ArrayColSymComplexSparse, Allocator1>& A,
1483  const Vector<T2, VectFull, Allocator2>& B,
1484  const T4& beta,
1485  Vector<T3, VectFull, Allocator3>& C)
1486  {
1487  T4 zero; T0 one;
1488  SetComplexZero(zero);
1489  SetComplexOne(one);
1490 
1491  if (beta == zero)
1492  C.Fill(0);
1493  else
1494  Mlt(beta, C);
1495 
1496  int m = A.GetM(), n, p;
1497  typename ClassComplexType<T1>::Treal val;
1498  T3 val_cplx;
1499  if (alpha == one)
1500  {
1501  for (int i = 0 ; i < m ; i++)
1502  {
1503  n = A.GetRealColumnSize(i);
1504  for (int k = 0; k < n ; k++)
1505  {
1506  p = A.IndexReal(i, k);
1507  val = A.ValueReal(i, k);
1508  if (p == i)
1509  C(i) += val * B(i);
1510  else
1511  {
1512  C(i) += val * B(p);
1513  C(p) += val * B(i);
1514  }
1515  }
1516  n = A.GetImagColumnSize(i);
1517  for (int k = 0; k < n ; k++)
1518  {
1519  p = A.IndexImag(i, k);
1520  val = A.ValueImag(i, k);
1521  if (p == i)
1522  C(i) += T1(0, val) * B(i);
1523  else
1524  {
1525  C(i) += T1(0, val) * B(p);
1526  C(p) += T1(0, val) * B(i);
1527  }
1528  }
1529  }
1530  }
1531  else // alpha != 1.
1532  {
1533  for (int i = 0 ; i < m ; i++)
1534  {
1535  n = A.GetRealColumnSize(i);
1536  for (int k = 0; k < n ; k++)
1537  {
1538  p = A.IndexReal(i, k);
1539  val_cplx = alpha * T1(A.ValueReal(i, k), 0);
1540  if (p == i)
1541  C(i) += val_cplx * B(i);
1542  else
1543  {
1544  C(i) += val_cplx * B(p);
1545  C(p) += val_cplx * B(i);
1546  }
1547  }
1548  n = A.GetImagColumnSize(i);
1549  for (int k = 0; k < n ; k++)
1550  {
1551  p = A.IndexImag(i, k);
1552  val_cplx = alpha * T1(0, A.ValueImag(i, k));
1553  if (p == i)
1554  C(i) += val_cplx * B(i);
1555  else
1556  {
1557  C(i) += val_cplx * B(p);
1558  C(p) += val_cplx * B(i);
1559  }
1560  }
1561  }
1562  }
1563  }
1564 
1565 
1566  template<class T0, class T1, class T2, class T3, class T4,
1567  class Allocator1, class Allocator2, class Allocator3>
1568  void MltAddVector(const T0& alpha,
1569  const SeldonTranspose& Trans, const Matrix<T1, Symmetric,
1570  ArrayColSymComplexSparse, Allocator1>& A,
1571  const Vector<T2, VectFull, Allocator2>& B,
1572  const T4& beta,
1573  Vector<T3, VectFull, Allocator3>& C)
1574  {
1575  if (!Trans.ConjTrans())
1576  {
1577  MltAddVector(alpha, A, B, beta, C);
1578  return;
1579  }
1580 
1581  T4 zero; T0 one;
1582  SetComplexZero(zero);
1583  SetComplexOne(one);
1584 
1585  if (beta == zero)
1586  C.Fill(0);
1587  else
1588  Mlt(beta, C);
1589 
1590  int m = A.GetM(),n,p;
1591  typename ClassComplexType<T1>::Treal val;
1592  T3 val_cplx;
1593  if (alpha == one)
1594  {
1595  for (int i = 0 ; i < m ; i++)
1596  {
1597  n = A.GetRealColumnSize(i);
1598  for (int k = 0; k < n ; k++)
1599  {
1600  p = A.IndexReal(i, k);
1601  val = A.ValueReal(i, k);
1602  if (p == i)
1603  C(i) += val * B(i);
1604  else
1605  {
1606  C(i) += val * B(p);
1607  C(p) += val * B(i);
1608  }
1609  }
1610  n = A.GetImagColumnSize(i);
1611  for (int k = 0; k < n ; k++)
1612  {
1613  p = A.IndexImag(i, k);
1614  val = A.ValueImag(i, k);
1615  if (p == i)
1616  C(i) -= T1(0, val) * B(i);
1617  else
1618  {
1619  C(i) -= T1(0, val) * B(p);
1620  C(p) -= T1(0, val) * B(i);
1621  }
1622  }
1623  }
1624  }
1625  else
1626  {
1627  // alpha different from 1
1628  for (int i = 0 ; i < m ; i++)
1629  {
1630  n = A.GetRealColumnSize(i);
1631  for (int k = 0; k < n ; k++)
1632  {
1633  p = A.IndexReal(i, k);
1634  val_cplx = alpha * T1(A.ValueReal(i, k), 0);
1635  if (p == i)
1636  C(i) += val_cplx * B(i);
1637  else
1638  {
1639  C(i) += val_cplx * B(p);
1640  C(p) += val_cplx * B(i);
1641  }
1642  }
1643  n = A.GetImagColumnSize(i);
1644  for (int k = 0; k < n ; k++)
1645  {
1646  p = A.IndexImag(i, k);
1647  val_cplx = alpha * T1(0, A.ValueImag(i, k));
1648  if (p == i)
1649  C(i) -= val_cplx * B(i);
1650  else
1651  {
1652  C(i) -= val_cplx * B(p);
1653  C(p) -= val_cplx * B(i);
1654  }
1655  }
1656  }
1657  }
1658  }
1659 
1660 
1661  /*** ArrayRowComplexSparse ***/
1662 
1663 
1664  template<class T0, class T1, class T2, class T3, class T4,
1665  class Allocator1, class Allocator2, class Allocator3>
1666  void MltAddVector(const T0& alpha, const Matrix<T1, General,
1667  ArrayRowComplexSparse, Allocator1>& A,
1668  const Vector<T2, VectFull, Allocator2>& B,
1669  const T4& beta,
1670  Vector<T3, VectFull, Allocator3>& C)
1671  {
1672  T4 zero; T0 one;
1673  SetComplexZero(zero);
1674  SetComplexOne(one);
1675 
1676  if (beta == zero)
1677  C.Fill(0);
1678  else
1679  Mlt(beta, C);
1680 
1681  int m = A.GetM(), n, p;
1682  typename ClassComplexType<T1>::Treal val;
1683  T3 val_cplx;
1684  if (alpha == one)
1685  {
1686  for (int i = 0 ; i < m ; i++)
1687  {
1688  n = A.GetRealRowSize(i);
1689  for (int k = 0; k < n ; k++)
1690  {
1691  p = A.IndexReal(i, k);
1692  val = A.ValueReal(i, k);
1693  C(i) += val * B(p);
1694  }
1695  n = A.GetImagRowSize(i);
1696  for (int k = 0; k < n ; k++)
1697  {
1698  p = A.IndexImag(i, k);
1699  val = A.ValueImag(i, k);
1700  C(i) += T1(0, val) * B(p);
1701  }
1702  }
1703  }
1704  else
1705  {
1706  // alpha different from 1
1707  for (int i = 0 ; i < m ; i++)
1708  {
1709  n = A.GetRealRowSize(i);
1710  for (int k = 0; k < n ; k++)
1711  {
1712  p = A.IndexReal(i, k);
1713  val_cplx = alpha * T1(A.ValueReal(i, k), 0);
1714  C(i) += val_cplx * B(p);
1715  }
1716  n = A.GetImagRowSize(i);
1717  for (int k = 0; k < n ; k++)
1718  {
1719  p = A.IndexImag(i, k);
1720  val_cplx = alpha * T1(0, A.ValueImag(i, k));
1721  C(i) += val_cplx * B(p);
1722  }
1723  }
1724  }
1725  }
1726 
1727 
1728  template<class T0, class T1, class T2, class T3, class T4,
1729  class Allocator1, class Allocator2, class Allocator3>
1730  void MltAddVector(const T0& alpha,
1731  const SeldonTranspose& Trans, const Matrix<T1, General,
1732  ArrayRowComplexSparse, Allocator1>& A,
1733  const Vector<T2, VectFull, Allocator2>& B,
1734  const T4& beta,
1735  Vector<T3, VectFull, Allocator3>& C)
1736  {
1737  if (Trans.NoTrans())
1738  {
1739  MltAddVector(alpha, A, B, beta, C);
1740  return;
1741  }
1742 
1743  T4 zero; T0 one;
1744  SetComplexZero(zero);
1745  SetComplexOne(one);
1746 
1747  if (beta == zero)
1748  C.Fill(0);
1749  else
1750  Mlt(beta, C);
1751 
1752  int m = A.GetM(),n,p;
1753  typename ClassComplexType<T1>::Treal val;
1754  T3 val_cplx;
1755 
1756  if (Trans.Trans())
1757  {
1758  if (alpha == one)
1759  {
1760  for (int i = 0 ; i < m ; i++)
1761  {
1762  n = A.GetRealRowSize(i);
1763  for (int k = 0; k < n ; k++)
1764  {
1765  p = A.IndexReal(i, k);
1766  val = A.ValueReal(i, k);
1767  C(p) += val * B(i);
1768  }
1769  n = A.GetImagRowSize(i);
1770  for (int k = 0; k < n ; k++)
1771  {
1772  p = A.IndexImag(i, k);
1773  val = A.ValueImag(i, k);
1774  C(p) += T1(0, val) * B(i);
1775  }
1776  }
1777  }
1778  else
1779  {
1780  // alpha different from 1
1781  for (int i = 0 ; i < m ; i++)
1782  {
1783  n = A.GetRealRowSize(i);
1784  for (int k = 0; k < n ; k++)
1785  {
1786  p = A.IndexReal(i, k);
1787  val_cplx = alpha * T1(A.ValueReal(i, k), 0);
1788  C(p) += val_cplx * B(i);
1789  }
1790 
1791  n = A.GetImagRowSize(i);
1792  for (int k = 0; k < n ; k++)
1793  {
1794  p = A.IndexImag(i, k);
1795  val_cplx = alpha * T1(0, A.ValueImag(i, k));
1796  C(p) += val_cplx * B(i);
1797  }
1798  }
1799  }
1800  }
1801  else
1802  {
1803  if (alpha == one)
1804  {
1805  for (int i = 0 ; i < m ; i++)
1806  {
1807  n = A.GetRealRowSize(i);
1808  for (int k = 0; k < n ; k++)
1809  {
1810  p = A.IndexReal(i, k);
1811  val = A.ValueReal(i, k);
1812  C(p) += val * B(i);
1813  }
1814  n = A.GetImagRowSize(i);
1815  for (int k = 0; k < n ; k++)
1816  {
1817  p = A.IndexImag(i, k);
1818  val = A.ValueImag(i, k);
1819  C(p) -= T1(0, val) * B(i);
1820  }
1821  }
1822  }
1823  else
1824  {
1825  // alpha different from 1
1826  for (int i = 0 ; i < m ; i++)
1827  {
1828  n = A.GetRealRowSize(i);
1829  for (int k = 0; k < n ; k++)
1830  {
1831  p = A.IndexReal(i, k);
1832  val_cplx = alpha * T1(A.ValueReal(i, k), 0);
1833  C(p) += val_cplx * B(i);
1834  }
1835  n = A.GetImagRowSize(i);
1836  for (int k = 0; k < n ; k++)
1837  {
1838  p = A.IndexImag(i, k);
1839  val_cplx = alpha * T1(0, -A.ValueImag(i, k));
1840  C(p) += val_cplx * B(i);
1841  }
1842  }
1843  }
1844  }
1845  }
1846 
1847 
1848  template<class T0, class T1, class T2, class T3, class T4,
1849  class Allocator1, class Allocator2, class Allocator3>
1850  void MltAddVector(const T0& alpha, const Matrix<T1, General,
1851  ArrayColComplexSparse, Allocator1>& A,
1852  const Vector<T2, VectFull, Allocator2>& B,
1853  const T4& beta,
1854  Vector<T3, VectFull, Allocator3>& C)
1855  {
1856  T4 zero; T0 one;
1857  SetComplexZero(zero);
1858  SetComplexOne(one);
1859 
1860  if (beta == zero)
1861  C.Fill(0);
1862  else
1863  Mlt(beta, C);
1864 
1865  int n, p;
1866  typename ClassComplexType<T1>::Treal val;
1867  T3 val_cplx;
1868  if (alpha == one)
1869  {
1870  for (int i = 0 ; i < A.GetN(); i++)
1871  {
1872  n = A.GetRealColumnSize(i);
1873  for (int k = 0; k < n ; k++)
1874  {
1875  p = A.IndexReal(i, k);
1876  val = A.ValueReal(i, k);
1877  C(p) += val * B(i);
1878  }
1879  n = A.GetImagColumnSize(i);
1880  for (int k = 0; k < n ; k++)
1881  {
1882  p = A.IndexImag(i, k);
1883  val = A.ValueImag(i, k);
1884  C(p) += T1(0, val) * B(i);
1885  }
1886  }
1887  }
1888  else
1889  {
1890  // alpha different from 1
1891  for (int i = 0; i < A.GetN(); i++)
1892  {
1893  n = A.GetRealColumnSize(i);
1894  for (int k = 0; k < n ; k++)
1895  {
1896  p = A.IndexReal(i, k);
1897  val_cplx = alpha * T1(A.ValueReal(i, k), 0);
1898  C(p) += val_cplx * B(i);
1899  }
1900  n = A.GetImagColumnSize(i);
1901  for (int k = 0; k < n ; k++)
1902  {
1903  p = A.IndexImag(i, k);
1904  val_cplx = alpha * T1(0, A.ValueImag(i, k));
1905  C(p) += val_cplx * B(i);
1906  }
1907  }
1908  }
1909  }
1910 
1911 
1912  template<class T0, class T1, class T2, class T3, class T4,
1913  class Allocator1, class Allocator2, class Allocator3>
1914  void MltAddVector(const T0& alpha,
1915  const SeldonTranspose& Trans, const Matrix<T1, General,
1916  ArrayColComplexSparse, Allocator1>& A,
1917  const Vector<T2, VectFull, Allocator2>& B,
1918  const T4& beta,
1919  Vector<T3, VectFull, Allocator3>& C)
1920  {
1921  if (Trans.NoTrans())
1922  {
1923  MltAddVector(alpha, A, B, beta, C);
1924  return;
1925  }
1926 
1927  T4 zero; T0 one;
1928  SetComplexZero(zero);
1929  SetComplexOne(one);
1930 
1931  if (beta == zero)
1932  C.Fill(0);
1933  else
1934  Mlt(beta, C);
1935 
1936  int n, p;
1937  typename ClassComplexType<T1>::Treal val;
1938  T3 val_cplx;
1939 
1940  if (Trans.Trans())
1941  {
1942  if (alpha == one)
1943  {
1944  for (int i = 0; i < A.GetN(); i++)
1945  {
1946  n = A.GetRealColumnSize(i);
1947  for (int k = 0; k < n ; k++)
1948  {
1949  p = A.IndexReal(i, k);
1950  val = A.ValueReal(i, k);
1951  C(i) += val * B(p);
1952  }
1953 
1954  n = A.GetImagColumnSize(i);
1955  for (int k = 0; k < n ; k++)
1956  {
1957  p = A.IndexImag(i, k);
1958  val = A.ValueImag(i, k);
1959  C(i) += T1(0, val) * B(p);
1960  }
1961  }
1962  }
1963  else
1964  {
1965  // alpha different from 1
1966  for (int i = 0; i < A.GetN(); i++)
1967  {
1968  n = A.GetRealColumnSize(i);
1969  for (int k = 0; k < n ; k++)
1970  {
1971  p = A.IndexReal(i, k);
1972  val_cplx = alpha * T1(A.ValueReal(i, k), 0);
1973  C(i) += val_cplx * B(p);
1974  }
1975 
1976  n = A.GetImagColumnSize(i);
1977  for (int k = 0; k < n ; k++)
1978  {
1979  p = A.IndexImag(i, k);
1980  val_cplx = alpha * T1(0, A.ValueImag(i, k));
1981  C(i) += val_cplx * B(p);
1982  }
1983  }
1984  }
1985  }
1986  else
1987  {
1988  if (alpha == one)
1989  {
1990  for (int i = 0; i < A.GetN(); i++)
1991  {
1992  n = A.GetRealColumnSize(i);
1993  for (int k = 0; k < n; k++)
1994  {
1995  p = A.IndexReal(i, k);
1996  val = A.ValueReal(i, k);
1997  C(i) += val * B(p);
1998  }
1999  n = A.GetImagColumnSize(i);
2000  for (int k = 0; k < n; k++)
2001  {
2002  p = A.IndexImag(i, k);
2003  val = A.ValueImag(i, k);
2004  C(i) -= T1(0, val) * B(p);
2005  }
2006  }
2007  }
2008  else
2009  {
2010  // alpha different from 1
2011  for (int i = 0 ; i < A.GetN(); i++)
2012  {
2013  n = A.GetRealColumnSize(i);
2014  for (int k = 0; k < n ; k++)
2015  {
2016  p = A.IndexReal(i, k);
2017  val_cplx = alpha * T1(A.ValueReal(i, k), 0);
2018  C(i) += val_cplx * B(p);
2019  }
2020  n = A.GetImagColumnSize(i);
2021  for (int k = 0; k < n; k++)
2022  {
2023  p = A.IndexImag(i, k);
2024  val_cplx = alpha * T1(0, -A.ValueImag(i, k));
2025  C(i) += val_cplx * B(p);
2026  }
2027  }
2028  }
2029  }
2030  }
2031 
2032 
2034  // SOR //
2035 
2036 
2038 
2045  template <class T0, class Prop0, class Allocator0,
2046  class T1, class Storage1, class Allocator1,
2047  class T2, class Storage2, class Allocator2, class T3>
2048  void SorVector(const Matrix<T0, Prop0, RowComplexSparse, Allocator0>& A,
2049  Vector<complex<T2>, Storage2, Allocator2>& X,
2050  const Vector<complex<T1>, Storage1, Allocator1>& B,
2051  const T3& omega, int iter, int type_ssor)
2052  {
2053  complex<T1> temp, zero; T3 one;
2054  SetComplexZero(zero);
2055  SetComplexOne(one);
2056 
2057  int ma = A.GetM();
2058 
2059 #ifdef SELDON_CHECK_BOUNDS
2060  int na = A.GetN();
2061  if (na != ma)
2062  throw WrongDim("SOR", "Matrix must be squared.");
2063 
2064  if (ma != X.GetM() || ma != B.GetM())
2065  throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
2066 #endif
2067 
2068  long* ptr_real = A.GetRealPtr();
2069  long* ptr_imag = A.GetImagPtr();
2070  int* ind_real = A.GetRealInd();
2071  int* ind_imag = A.GetImagInd();
2072 
2073  typename ClassComplexType<T0>::Treal* data_real = A.GetRealData();
2074  typename ClassComplexType<T0>::Treal* data_imag = A.GetImagData();
2075 
2076  T0 ajj;
2077 
2078  if (type_ssor % 2 == 0)
2079  for (int i = 0; i < iter; i++)
2080  for (int j = 0; j < ma; j++)
2081  {
2082  temp = zero;
2083  ajj = zero;
2084  for (long k = ptr_real[j]; k < ptr_real[j+1]; k++)
2085  {
2086  if (ind_real[k] != j)
2087  temp += data_real[k] * X(ind_real[k]);
2088  else
2089  SetComplexReal(data_real[k], ajj);
2090  }
2091 
2092  for (long k = ptr_imag[j]; k < ptr_imag[j+1]; k++)
2093  {
2094  if (ind_imag[k] != j)
2095  temp += complex<T1>(0, data_imag[k]) * X(ind_imag[k]);
2096  else
2097  ajj += complex<T1>(0, data_imag[k]);
2098  }
2099 
2100 #ifdef SELDON_CHECK_BOUNDS
2101  if (ajj == zero)
2102  throw WrongArgument("SOR", "Matrix must contain"
2103  " a non-null diagonal");
2104 #endif
2105 
2106  X(j) = (one-omega) * X(j) + omega * (B(j) - temp) / ajj;
2107  }
2108 
2109  if (type_ssor % 3 == 0)
2110  for (int i = 0; i < iter; i++)
2111  for (int j = ma-1; j >= 0; j--)
2112  {
2113  temp = zero;
2114  ajj = zero;
2115  for (long k = ptr_real[j]; k < ptr_real[j+1]; k++)
2116  {
2117  if (ind_real[k] != j)
2118  temp += data_real[k] * X(ind_real[k]);
2119  else
2120  SetComplexReal(data_real[k], ajj);
2121  }
2122 
2123  for (long k = ptr_imag[j]; k < ptr_imag[j+1]; k++)
2124  {
2125  if (ind_imag[k] != j)
2126  temp += complex<T1>(0, data_imag[k]) * X(ind_imag[k]);
2127  else
2128  ajj += complex<T1>(0, data_imag[k]);
2129  }
2130 
2131  X(j) = (one-omega) * X(j) + omega * (B(j) - temp) / ajj;
2132  }
2133  }
2134 
2135 
2137 
2144  template <class T0, class Prop0, class Allocator0,
2145  class T1, class Storage1, class Allocator1,
2146  class T2, class Storage2, class Allocator2, class T3>
2147  void SorVector(const Matrix<T0, Prop0, ArrayRowComplexSparse, Allocator0>& A,
2148  Vector<complex<T2>, Storage2, Allocator2>& X,
2149  const Vector<complex<T1>, Storage1, Allocator1>& B,
2150  const T3& omega, int iter, int type_ssor)
2151  {
2152  complex<T1> temp, zero; T3 one;
2153  SetComplexZero(zero);
2154  SetComplexOne(one);
2155 
2156  int ma = A.GetM();
2157 
2158 #ifdef SELDON_CHECK_BOUNDS
2159  int na = A.GetN();
2160  if (na != ma)
2161  throw WrongDim("SOR", "Matrix must be squared.");
2162 
2163  if (ma != X.GetM() || ma != B.GetM())
2164  throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
2165 #endif
2166 
2167  T0 ajj;
2168 
2169  if (type_ssor%2 == 0)
2170  for (int i = 0; i < iter; i++)
2171  for (int j = 0; j < ma; j++)
2172  {
2173  temp = zero;
2174  ajj = zero;
2175  for (int k = 0; k < A.GetRealRowSize(j); k++)
2176  {
2177  if (A.IndexReal(j, k) != j)
2178  temp += A.ValueReal(j,k) * X(A.IndexReal(j,k));
2179  else
2180  ajj += A.ValueReal(j, k);
2181  }
2182 
2183  for (int k = 0; k < A.GetImagRowSize(j); k++)
2184  {
2185  if (A.IndexImag(j, k) != j)
2186  temp += T0(0, A.ValueImag(j,k))
2187  * X(A.IndexImag(j, k));
2188  else
2189  ajj += T0(0, A.ValueImag(j,k));
2190  }
2191 
2192 #ifdef SELDON_CHECK_BOUNDS
2193  if (ajj == zero)
2194  throw WrongArgument("SOR", "Matrix must contain"
2195  " a non-null diagonal");
2196 #endif
2197 
2198  X(j) = (one-omega) * X(j) + omega * (B(j) - temp) / ajj;
2199  }
2200 
2201  if (type_ssor % 3 == 0)
2202  for (int i = 0; i < iter; i++)
2203  for (int j = ma-1; j >= 0; j--)
2204  {
2205  temp = zero;
2206  ajj = zero;
2207  for (int k = 0; k < A.GetRealRowSize(j); k++)
2208  {
2209  if (A.IndexReal(j, k) != j)
2210  temp += A.ValueReal(j,k) * X(A.IndexReal(j,k));
2211  else
2212  ajj += A.ValueReal(j, k);
2213  }
2214 
2215  for (int k = 0; k < A.GetImagRowSize(j); k++)
2216  {
2217  if (A.IndexImag(j, k) != j)
2218  temp += T0(0, A.ValueImag(j,k))
2219  * X(A.IndexImag(j, k));
2220  else
2221  ajj += T0(0, A.ValueImag(j,k));
2222  }
2223 
2224  X(j) = (one-omega) * X(j) + omega * (B(j) - temp) / ajj;
2225  }
2226  }
2227 
2228 
2230 
2237  template <class T0, class Prop0, class Allocator0,
2238  class T1, class Storage1, class Allocator1,
2239  class T2, class Storage2, class Allocator2, class T3>
2240  void SorVector(const Matrix<T0, Prop0, ColComplexSparse, Allocator0>& A,
2241  Vector<complex<T2>, Storage2, Allocator2>& X,
2242  const Vector<complex<T1>, Storage1, Allocator1>& B,
2243  const T3& omega,int iter, int type_ssor)
2244  {
2245  complex<T1> zero; T3 one;
2246  SetComplexZero(zero);
2247  SetComplexOne(one);
2248 
2249  int ma = A.GetM();
2250 
2251 #ifdef SELDON_CHECK_BOUNDS
2252  int na = A.GetN();
2253  if (na != ma)
2254  throw WrongDim("SOR", "Matrix must be squared.");
2255 
2256  if (ma != X.GetM() || ma != B.GetM())
2257  throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
2258 #endif
2259 
2260  long* ptr_real = A.GetRealPtr();
2261  long* ind_real = A.GetRealInd();
2262  typename ClassComplexType<T0>::Treal* data_real = A.GetRealData();
2263 
2264  int* ptr_imag = A.GetImagPtr();
2265  int* ind_imag = A.GetImagInd();
2266  typename ClassComplexType<T0>::Treal* data_imag = A.GetImagData();
2267 
2268  // Let us consider the following splitting : A = D - L - U
2269  // D diagonal of A
2270  // L lower part of A
2271  // U upper part of A
2272 
2273  // Forward sweep
2274  // (D/omega - L) X^{n+1/2} = (U + (1-omega)/omega D) X^n + B
2275  T0 ajj;
2276  T3 coef = (one - omega) / omega;
2277  if (type_ssor % 2 == 0)
2278  for (int i = 0; i < iter; i++)
2279  {
2280  // First we compute X = (U + (1-omega)/omega D) X + B
2281  for (int j = 0; j < ma; j++)
2282  {
2283  long kr = ptr_real[j];
2284  while ( (kr < ptr_real[j+1]) && (ind_real[kr] < j))
2285  {
2286  X(ind_real[kr]) -= data_real[kr]*X(j);
2287  kr++;
2288  }
2289 
2290  if ( (kr < ptr_real[j+1]) && (ind_real[kr] == j))
2291  SetComplexReal(data_real[kr], ajj);
2292 
2293  long ki = ptr_imag[j];
2294  while ( (ki < ptr_imag[j+1]) && (ind_imag[ki] < j))
2295  {
2296  X(ind_imag[ki]) -= T0(0, data_imag[ki])*X(j);
2297  ki++;
2298  }
2299 
2300  if ( (ki < ptr_imag[j+1]) && (ind_imag[ki] == j))
2301  ajj += T0(0, data_imag[ki]);
2302 
2303 #ifdef SELDON_CHECK_BOUNDS
2304  if (ajj == zero)
2305  throw WrongArgument("SOR", "Matrix must contain"
2306  " a non-null diagonal");
2307 #endif
2308 
2309  X(j) = B(j) + coef * ajj * X(j);
2310  }
2311 
2312  // Then we solve (D/omega - L) X = X
2313  for (int j = 0; j < ma; j++)
2314  {
2315  long kr = ptr_real[j];
2316  while ( (kr < ptr_real[j+1]) && (ind_real[kr] < j))
2317  kr++;
2318 
2319  if ( (kr < ptr_real[j+1]) && (ind_real[kr] == j))
2320  {
2321  SetComplexReal(data_real[kr], ajj);
2322  kr++;
2323  }
2324 
2325  long ki = ptr_imag[j];
2326  while ( (ki < ptr_imag[j+1]) && (ind_imag[ki] < j))
2327  ki++;
2328 
2329  if ( (ki < ptr_imag[j+1]) && (ind_imag[ki] == j))
2330  {
2331  ajj += T0(0, data_imag[ki]);
2332  ki++;
2333  }
2334 
2335  X(j) *= omega/ajj;
2336  while (kr < ptr_real[j+1])
2337  {
2338  X(ind_real[kr]) -= data_real[kr]*X(j);
2339  kr++;
2340  }
2341 
2342  while (ki < ptr_imag[j+1])
2343  {
2344  X(ind_imag[ki]) -= T0(0, data_imag[ki])*X(j);
2345  ki++;
2346  }
2347  }
2348  }
2349 
2350  // Backward sweep.
2351  // (D/omega - U) X^{n+1} = (L + (1-omega)/omega D) X^{n+1/2} + B
2352  if (type_ssor % 3 == 0)
2353  for (int i = 0; i < iter; i++)
2354  {
2355  // First we compute X = (L + (1-omega)/omega D) X + B.
2356  for (int j = ma-1; j >= 0; j--)
2357  {
2358  long kr = ptr_real[j+1]-1;
2359  while ( (kr >= ptr_real[j]) && (ind_real[kr] > j) )
2360  {
2361  X(ind_real[kr]) -= data_real[kr]*X(j);
2362  kr--;
2363  }
2364 
2365  if ( (kr >= ptr_real[j]) && (ind_real[kr] == j))
2366  SetComplexReal(data_real[kr], ajj);
2367 
2368  long ki = ptr_imag[j+1]-1;
2369  while ( (ki >= ptr_imag[j]) && (ind_imag[ki] > j) )
2370  {
2371  X(ind_imag[ki]) -= T0(0, data_imag[ki])*X(j);
2372  ki--;
2373  }
2374 
2375  if ( (ki >= ptr_imag[j]) && (ind_imag[ki] == j))
2376  ajj += T0(0, data_imag[ki]);
2377 
2378 #ifdef SELDON_CHECK_BOUNDS
2379  if (ajj == zero)
2380  throw WrongArgument("SOR", "Matrix must contain"
2381  " a non-null diagonal");
2382 #endif
2383 
2384  X(j) = B(j) + coef * ajj * X(j);
2385  }
2386 
2387  // Then we solve (D/omega - U) X = X.
2388  for (int j = ma-1; j >= 0; j--)
2389  {
2390  long kr = ptr_real[j+1]-1;
2391  while ( (kr >= ptr_real[j]) && (ind_real[kr] > j) )
2392  kr--;
2393 
2394  if ( (kr >= ptr_real[j]) && (ind_real[kr] == j))
2395  {
2396  ajj = T0(data_real[kr], 0);
2397  kr--;
2398  }
2399 
2400  long ki = ptr_imag[j+1]-1;
2401  while ( (ki >= ptr_imag[j]) && (ind_imag[ki] > j) )
2402  ki--;
2403 
2404  if ( (ki >= ptr_imag[j]) && (ind_imag[ki] == j))
2405  {
2406  ajj += T0(0, data_imag[ki]);
2407  ki--;
2408  }
2409 
2410  X(j) *= omega/ajj;
2411  while (kr >= ptr_real[j])
2412  {
2413  X(ind_real[kr]) -= data_real[kr]*X(j);
2414  kr--;
2415  }
2416 
2417  while (ki >= ptr_imag[j])
2418  {
2419  X(ind_imag[ki]) -= T0(0, data_imag[ki])*X(j);
2420  ki--;
2421  }
2422  }
2423  }
2424  }
2425 
2426 
2428 
2435  template <class T0, class Prop0, class Allocator0,
2436  class T1, class Storage1, class Allocator1,
2437  class T2, class Storage2, class Allocator2, class T3>
2438  void SorVector(const Matrix<T0, Prop0, ArrayColComplexSparse, Allocator0>& A,
2439  Vector<complex<T2>, Storage2, Allocator2>& X,
2440  const Vector<complex<T1>, Storage1, Allocator1>& B,
2441  const T3& omega, int iter, int type_ssor)
2442  {
2443  complex<T1> zero; T3 one;
2444  SetComplexZero(zero);
2445  SetComplexOne(one);
2446 
2447  int ma = A.GetM();
2448 
2449 #ifdef SELDON_CHECK_BOUNDS
2450  int na = A.GetN();
2451  if (na != ma)
2452  throw WrongDim("SOR", "Matrix must be squared.");
2453 
2454  if (ma != X.GetM() || ma != B.GetM())
2455  throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
2456 #endif
2457 
2458  // Let us consider the following splitting : A = D - L - U
2459  // D diagonal of A
2460  // L lower part of A
2461  // U upper part of A
2462 
2463  // Forward sweep
2464  // (D/omega - L) X^{n+1/2} = (U + (1-omega)/omega D) X^n + B
2465  T0 ajj;
2466  T3 coef = (one - omega) / omega;
2467  if (type_ssor % 2 == 0)
2468  for (int i = 0; i < iter; i++)
2469  {
2470  // First we compute X = (U + (1-omega)/omega D) X + B
2471  for (int j = 0; j < ma; j++)
2472  {
2473  long kr = 0;
2474  while ( (kr < A.GetRealColumnSize(j)) && (A.IndexReal(j, kr) < j))
2475  {
2476  X(A.IndexReal(j, kr)) -= A.ValueReal(j, kr) * X(j);
2477  kr++;
2478  }
2479 
2480  if ( (kr < A.GetRealColumnSize(j)) && (A.IndexReal(j, kr) == j))
2481  ajj = T0(A.ValueReal(j, kr), 0);
2482 
2483  long ki = 0;
2484  while ( (ki < A.GetImagColumnSize(j)) && (A.IndexImag(j, ki) < j))
2485  {
2486  X(A.IndexImag(j, ki)) -= T0(0, A.ValueImag(j, ki))*X(j);
2487  ki++;
2488  }
2489 
2490  if ( (ki < A.GetImagColumnSize(j)) && (A.IndexImag(j, ki) == j) )
2491  ajj += T0(0, A.ValueImag(j, ki));
2492 
2493 #ifdef SELDON_CHECK_BOUNDS
2494  if (ajj == zero)
2495  throw WrongArgument("SOR", "Matrix must contain"
2496  " a non-null diagonal");
2497 #endif
2498 
2499  X(j) = B(j) + coef * ajj * X(j);
2500  }
2501 
2502  // Then we solve (D/omega - L) X = X
2503  for (int j = 0; j < ma; j++)
2504  {
2505  long kr = 0;
2506  while ( (kr < A.GetRealColumnSize(j)) && (A.IndexReal(j, kr) < j))
2507  kr++;
2508 
2509  if ( (kr < A.GetRealColumnSize(j)) && (A.IndexReal(j, kr) == j) )
2510  {
2511  ajj = T0(A.ValueReal(j, kr), 0);
2512  kr++;
2513  }
2514 
2515  long ki = 0;
2516  while ( (ki < A.GetImagColumnSize(j)) && (A.IndexImag(j, ki) < j))
2517  ki++;
2518 
2519  if ( (ki < A.GetImagColumnSize(j)) && (A.IndexImag(j, ki) == j))
2520  {
2521  ajj += T0(0, A.ValueImag(j, ki));
2522  ki++;
2523  }
2524 
2525  X(j) *= omega/ajj;
2526  while (kr < A.GetRealColumnSize(j))
2527  {
2528  X(A.IndexReal(j, kr)) -= A.ValueReal(j, kr)*X(j);
2529  kr++;
2530  }
2531 
2532  while (ki < A.GetImagColumnSize(j))
2533  {
2534  X(A.IndexImag(j, ki)) -= T0(0, A.ValueImag(j, ki))*X(j);
2535  ki++;
2536  }
2537  }
2538  }
2539 
2540  // Backward sweep.
2541  // (D/omega - U) X^{n+1} = (L + (1-omega)/omega D) X^{n+1/2} + B
2542  if (type_ssor % 3 == 0)
2543  for (int i = 0; i < iter; i++)
2544  {
2545  // First we compute X = (L + (1-omega)/omega D) X + B.
2546  for (int j = ma-1; j >= 0; j--)
2547  {
2548  long kr = A.GetRealColumnSize(j)-1;
2549  while ( (kr >= 0) && (A.IndexReal(j, kr) > j) )
2550  {
2551  X(A.IndexReal(j, kr)) -= A.ValueReal(j, kr)*X(j);
2552  kr--;
2553  }
2554 
2555  if ( (kr >= 0) && (A.IndexReal(j, kr) == j))
2556  ajj = T0(A.ValueReal(j, kr), 0);
2557 
2558  long ki = A.GetImagColumnSize(j)-1;
2559  while ( (ki >= 0) && (A.IndexImag(j, ki) > j) )
2560  {
2561  X(A.IndexImag(j, ki)) -= T0(0, A.ValueImag(j, ki))*X(j);
2562  ki--;
2563  }
2564 
2565  if ( (ki >= 0) && (A.IndexImag(j, ki) == j))
2566  ajj += T0(0, A.ValueImag(j, ki));
2567 
2568 #ifdef SELDON_CHECK_BOUNDS
2569  if (ajj == zero)
2570  throw WrongArgument("SOR", "Matrix must contain"
2571  " a non-null diagonal");
2572 #endif
2573 
2574  X(j) = B(j) + coef * ajj * X(j);
2575  }
2576 
2577  // Then we solve (D/omega - U) X = X.
2578  for (int j = ma-1; j >= 0; j--)
2579  {
2580  long kr = A.GetRealColumnSize(j)-1;
2581  while ( (kr >= 0) && (A.IndexReal(j, kr) > j) )
2582  kr--;
2583 
2584  if ( (kr >= 0) && (A.IndexReal(j, kr) == j))
2585  {
2586  ajj = T0(A.ValueReal(j, kr), 0);
2587  kr--;
2588  }
2589 
2590  long ki = A.GetImagColumnSize(j)-1;
2591  while ( (ki >= 0) && (A.IndexImag(j, ki) > j) )
2592  ki--;
2593 
2594  if ( (ki >= 0) && (A.IndexImag(j, ki) == j))
2595  {
2596  ajj += T0(0, A.ValueImag(j, ki));
2597  ki--;
2598  }
2599 
2600  X(j) *= omega/ajj;
2601  while (kr >= 0)
2602  {
2603  X(A.IndexReal(j, kr)) -= A.ValueReal(j, kr)*X(j);
2604  kr--;
2605  }
2606 
2607  while (ki >= 0)
2608  {
2609  X(A.IndexImag(j, ki)) -= T0(0, A.ValueImag(j, ki))*X(j);
2610  ki--;
2611  }
2612  }
2613  }
2614  }
2615 
2616 
2618 
2625  template <class T0, class Prop0, class Allocator0,
2626  class T1, class Storage1, class Allocator1,
2627  class T2, class Storage2, class Allocator2, class T3>
2628  void SorVector(const Matrix<T0, Prop0, RowSymComplexSparse, Allocator0>& A,
2629  Vector<complex<T2>, Storage2, Allocator2>& X,
2630  const Vector<complex<T1>, Storage1, Allocator1>& B,
2631  const T3& omega, int iter, int type_ssor)
2632  {
2633  complex<T1> temp, zero; T3 one;
2634  SetComplexZero(zero);
2635  SetComplexOne(one);
2636 
2637  int ma = A.GetM();
2638 
2639 #ifdef SELDON_CHECK_BOUNDS
2640  int na = A.GetN();
2641  if (na != ma)
2642  throw WrongDim("SOR", "Matrix must be squared.");
2643 
2644  if (ma != X.GetM() || ma != B.GetM())
2645  throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
2646 #endif
2647 
2648  long* ptr_real = A.GetRealPtr();
2649  int* ind_real = A.GetRealInd();
2650  typename ClassComplexType<T0>::Treal* data_real = A.GetRealData();
2651 
2652  long* ptr_imag = A.GetImagPtr();
2653  int* ind_imag = A.GetImagInd();
2654  typename ClassComplexType<T0>::Treal* data_imag = A.GetImagData();
2655 
2656  T0 ajj;
2657 
2658  // Let us consider the following splitting : A = D - L - U
2659  // D diagonal of A
2660  // L lower part of A
2661  // U upper part of A, A is symmetric, so L = U^t
2662 
2663  // Forward sweep
2664  // (D/omega - L) X^{n+1/2} = (U + (1-omega)/omega D) X^n + B
2665  T3 coef = (T3(1) - omega) / omega;
2666  if (type_ssor % 2 == 0)
2667  for (int i = 0; i < iter; i++)
2668  {
2669  // First we do X = (U + (1-omega)/omega D) X + B
2670  for (int j = 0; j < ma; j++)
2671  {
2672  temp = zero;
2673  ajj = zero;
2674  long kr = ptr_real[j];
2675  if ((kr < ptr_real[j+1]) && (ind_real[kr] == j))
2676  ajj = T0(data_real[kr++], 0);
2677 
2678  long ki = ptr_imag[j];
2679  if ((ki < ptr_imag[j+1]) && (ind_imag[ki] == j))
2680  ajj += T0(0, data_imag[ki++]);
2681 
2682 #ifdef SELDON_CHECK_BOUNDS
2683  if ( ajj == zero)
2684  throw WrongArgument("SOR", "Matrix must contain"
2685  " a non-null diagonal");
2686 #endif
2687 
2688  for (long k = kr; k < ptr_real[j+1]; k++)
2689  temp += data_real[k] * X(ind_real[k]);
2690 
2691  for (long k = ki; k < ptr_imag[j+1]; k++)
2692  temp += T0(0, data_imag[k]) * X(ind_imag[k]);
2693 
2694  X(j) = coef * ajj * X(j) + B(j) - temp;
2695  }
2696 
2697  // Then we solve (D/omega - L) X = X
2698  for (int j = 0; j < ma; j++)
2699  {
2700  ajj = zero;
2701  long kr = ptr_real[j];
2702  if ((kr < ptr_real[j+1]) && (ind_real[kr] == j))
2703  ajj = T0(data_real[kr++], 0);
2704 
2705  long ki = ptr_imag[j];
2706  if ((ki < ptr_imag[j+1]) && (ind_imag[ki] == j))
2707  ajj += T0(0, data_imag[ki++]);
2708 
2709  X(j) *= omega / ajj;
2710  for (long k = kr; k < ptr_real[j+1]; k++)
2711  X(ind_real[k]) -= data_real[k]*X(j);
2712 
2713  for (long k = ki; k < ptr_imag[j+1]; k++)
2714  X(ind_imag[k]) -= T0(0, data_imag[k])*X(j);
2715  }
2716  }
2717 
2718 
2719  // Backward sweep.
2720  // (D/omega - U) X^{n+1} = (L + (1-omega)/omega D) X^{n+1/2} + B
2721  if (type_ssor % 3 == 0)
2722  for (int i = 0; i < iter; i++)
2723  {
2724  // First we compute X = (L + (1-omega)/omega D) X + B.
2725  for (int j = ma-1; j >= 0; j--)
2726  {
2727  ajj = zero;
2728  long kr = ptr_real[j];
2729  if ((kr < ptr_real[j+1]) && (ind_real[kr] == j))
2730  ajj = T0(data_real[kr++], 0);
2731 
2732  long ki = ptr_imag[j];
2733  if ((ki < ptr_imag[j+1]) && (ind_imag[ki] == j))
2734  ajj += T0(0, data_imag[ki++]);
2735 
2736  for (long k = kr; k < ptr_real[j+1]; k++)
2737  X(ind_real[k]) -= data_real[k]*X(j);
2738 
2739  for (long k = ki; k < ptr_imag[j+1]; k++)
2740  X(ind_imag[k]) -= T0(0, data_imag[k])*X(j);
2741 
2742  X(j) = B(j) + coef * ajj * X(j);
2743  }
2744 
2745  // Then we solve (D/omega - U) X = X.
2746  for (int j = ma-1; j >= 0; j--)
2747  {
2748  temp = zero;
2749  ajj = zero;
2750  long kr = ptr_real[j];
2751  if ((kr < ptr_real[j+1]) && (ind_real[kr] == j))
2752  ajj = T0(data_real[kr++], 0);
2753 
2754  long ki = ptr_imag[j];
2755  if ((ki < ptr_imag[j+1]) && (ind_imag[ki] == j))
2756  ajj += T0(0, data_imag[ki++]);
2757 
2758  for (long k = kr; k < ptr_real[j+1]; k++)
2759  temp += data_real[k]*X(ind_real[k]);
2760 
2761  for (long k = ki; k < ptr_imag[j+1]; k++)
2762  temp += T0(0, data_imag[k])*X(ind_imag[k]);
2763 
2764  X(j) = (X(j) - temp) * omega / ajj;
2765  }
2766  }
2767  }
2768 
2769 
2771 
2778  template <class T0, class Prop0, class Allocator0,
2779  class T1, class Storage1, class Allocator1,
2780  class T2, class Storage2, class Allocator2, class T3>
2781  void SorVector(const Matrix<T0, Prop0, ArrayRowSymComplexSparse, Allocator0>& A,
2782  Vector<complex<T2>, Storage2, Allocator2>& X,
2783  const Vector<complex<T1>, Storage1, Allocator1>& B,
2784  const T3& omega, int iter, int type_ssor)
2785  {
2786  complex<T1> temp, zero; T3 one;
2787  SetComplexZero(zero);
2788  SetComplexOne(one);
2789 
2790  int ma = A.GetM();
2791 
2792 #ifdef SELDON_CHECK_BOUNDS
2793  int na = A.GetN();
2794  if (na != ma)
2795  throw WrongDim("SOR", "Matrix must be squared.");
2796 
2797  if (ma != X.GetM() || ma != B.GetM())
2798  throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
2799 #endif
2800 
2801  T0 ajj;
2802 
2803  // Let us consider the following splitting : A = D - L - U
2804  // D diagonal of A
2805  // L lower part of A
2806  // U upper part of A, A is symmetric, so L = U^t
2807 
2808  // Forward sweep
2809  // (D/omega - L) X^{n+1/2} = (U + (1-omega)/omega D) X^n + B
2810  T3 coef = (one - omega) / omega;
2811  if (type_ssor % 2 == 0)
2812  for (int i = 0; i < iter; i++)
2813  {
2814  // First we do X = (U + (1-omega)/omega D) X + B
2815  for (int j = 0; j < ma; j++)
2816  {
2817  temp = zero;
2818  ajj = zero;
2819  int kr = 0;
2820  if ((kr < A.GetRealRowSize(j)) && (A.IndexReal(j, kr) == j))
2821  ajj = T0(A.ValueReal(j, kr++), 0);
2822 
2823  int ki = 0;
2824  if ((ki < A.GetImagRowSize(j)) && (A.IndexImag(j, ki) == j))
2825  ajj += T0(0, A.ValueImag(j, ki++));
2826 
2827 #ifdef SELDON_CHECK_BOUNDS
2828  if ( ajj == zero)
2829  throw WrongArgument("SOR", "Matrix must contain"
2830  " a non-null diagonal");
2831 #endif
2832 
2833  for (int k = kr; k < A.GetRealRowSize(j); k++)
2834  temp += A.ValueReal(j, k) * X(A.IndexReal(j, k));
2835 
2836  for (int k = ki; k < A.GetImagRowSize(j); k++)
2837  temp += T0(0, A.ValueImag(j, k)) * X(A.IndexImag(j, k));
2838 
2839  X(j) = coef * ajj * X(j) + B(j) - temp;
2840  }
2841 
2842  // Then we solve (D/omega - L) X = X
2843  for (int j = 0; j < ma; j++)
2844  {
2845  ajj = zero;
2846  int kr = 0;
2847  if ((kr < A.GetRealRowSize(j)) && (A.IndexReal(j, kr) == j))
2848  ajj = T0(A.ValueReal(j, kr++), 0);
2849 
2850  int ki = 0;
2851  if ((ki < A.GetImagRowSize(j)) && (A.IndexImag(j, ki) == j))
2852  ajj += T0(0, A.ValueImag(j, ki++));
2853 
2854  X(j) *= omega / ajj;
2855  for (int k = kr; k < A.GetRealRowSize(j); k++)
2856  X(A.IndexReal(j, k)) -= A.ValueReal(j, k)*X(j);
2857 
2858  for (int k = ki; k < A.GetImagRowSize(j); k++)
2859  X(A.IndexImag(j, k)) -= T0(0, A.ValueImag(j, k))*X(j);
2860  }
2861  }
2862 
2863 
2864  // Backward sweep.
2865  // (D/omega - U) X^{n+1} = (L + (1-omega)/omega D) X^{n+1/2} + B
2866  if (type_ssor % 3 == 0)
2867  for (int i = 0; i < iter; i++)
2868  {
2869  // First we compute X = (L + (1-omega)/omega D) X + B.
2870  for (int j = ma-1; j >= 0; j--)
2871  {
2872  ajj = zero;
2873  int kr = 0;
2874  if ((kr < A.GetRealRowSize(j)) && (A.IndexReal(j, kr) == j))
2875  ajj = T0(A.ValueReal(j, kr++), 0);
2876 
2877  int ki = 0;
2878  if ((ki < A.GetImagRowSize(j)) && (A.IndexImag(j, ki) == j))
2879  ajj += T0(0, A.ValueImag(j, ki++));
2880 
2881  for (int k = kr; k < A.GetRealRowSize(j); k++)
2882  X(A.IndexReal(j, k)) -= A.ValueReal(j, k)*X(j);
2883 
2884  for (int k = ki; k < A.GetImagRowSize(j); k++)
2885  X(A.IndexImag(j, k)) -= T0(0, A.ValueImag(j, k))*X(j);
2886 
2887  X(j) = B(j) + coef * ajj * X(j);
2888  }
2889 
2890  // Then we solve (D/omega - U) X = X.
2891  for (int j = ma-1; j >= 0; j--)
2892  {
2893  temp = zero;
2894  ajj = zero;
2895  int kr = 0;
2896  if ((kr < A.GetRealRowSize(j)) && (A.IndexReal(j, kr) == j))
2897  ajj = T0(A.ValueReal(j, kr++), 0);
2898 
2899  int ki = 0;
2900  if ((ki < A.GetImagRowSize(j)) && (A.IndexImag(j, ki) == j))
2901  ajj += T0(0, A.ValueImag(j, ki++));
2902 
2903  for (int k = kr; k < A.GetRealRowSize(j); k++)
2904  temp += A.ValueReal(j, k)*X(A.IndexReal(j, k));
2905 
2906  for (int k = ki; k < A.GetImagRowSize(j); k++)
2907  temp += T0(0, A.ValueImag(j, k))*X(A.IndexImag(j, k));
2908 
2909  X(j) = (X(j) - temp) * omega / ajj;
2910  }
2911  }
2912  }
2913 
2914 
2916 
2923  template <class T0, class Prop0, class Allocator0,
2924  class T1, class Storage1, class Allocator1,
2925  class T2, class Storage2, class Allocator2, class T3>
2926  void SorVector(const Matrix<T0, Prop0, ColSymComplexSparse, Allocator0>& A,
2927  Vector<complex<T2>, Storage2, Allocator2>& X,
2928  const Vector<complex<T1>, Storage1, Allocator1>& B,
2929  const T3& omega, int iter, int type_ssor)
2930  {
2931  complex<T1> temp, zero; T3 one;
2932  SetComplexZero(zero);
2933  SetComplexOne(one);
2934 
2935  int ma = A.GetM();
2936 
2937 #ifdef SELDON_CHECK_BOUNDS
2938  int na = A.GetN();
2939  if (na != ma)
2940  throw WrongDim("SOR", "Matrix must be squared.");
2941 
2942  if (ma != X.GetM() || ma != B.GetM())
2943  throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
2944 #endif
2945 
2946  long* ptr_real = A.GetRealPtr();
2947  int* ind_real = A.GetRealInd();
2948  typename ClassComplexType<T0>::Treal* data_real = A.GetRealData();
2949 
2950  long* ptr_imag = A.GetImagPtr();
2951  int* ind_imag = A.GetImagInd();
2952  typename ClassComplexType<T0>::Treal* data_imag = A.GetImagData();
2953 
2954  T0 ajj;
2955 
2956  // Let us consider the following splitting : A = D - L - U
2957  // D diagonal of A
2958  // L lower part of A
2959  // U upper part of A, A is symmetric, so L = U^t
2960 
2961  // Forward sweep
2962  // (D/omega - L) X^{n+1/2} = (U + (1-omega)/omega D) X^n + B
2963  T3 coef = (one - omega) / omega;
2964  if (type_ssor % 2 == 0)
2965  for (int i = 0; i < iter; i++)
2966  {
2967  // First we do X = (U + (1-omega)/omega D) X + B
2968  for (int j = 0; j < ma; j++)
2969  {
2970  ajj = zero;
2971  long kr = ptr_real[j+1]-1;
2972  if ((kr >= ptr_real[j]) && (ind_real[kr] == j))
2973  ajj = T0(data_real[kr--], 0);
2974 
2975  long ki = ptr_imag[j+1]-1;
2976  if ((ki >= ptr_imag[j]) && (ind_imag[ki] == j))
2977  ajj += T0(0, data_imag[ki--]);
2978 
2979 #ifdef SELDON_CHECK_BOUNDS
2980  if ( ajj == zero)
2981  throw WrongArgument("SOR", "Matrix must contain"
2982  " a non-null diagonal");
2983 #endif
2984 
2985  for (long k = ptr_real[j]; k <= kr; k++)
2986  X(ind_real[k]) -= data_real[k] * X(j);
2987 
2988  for (long k = ptr_imag[j]; k <= ki; k++)
2989  X(ind_imag[k]) -= T0(0, data_imag[k]) * X(j);
2990 
2991  X(j) = coef * ajj * X(j) + B(j);
2992  }
2993 
2994  // Then we solve (D/omega - L) X = X
2995  for (int j = 0; j < ma; j++)
2996  {
2997  ajj = zero;
2998  long kr = ptr_real[j+1]-1;
2999  if ((kr >= ptr_real[j]) && (ind_real[kr] == j))
3000  ajj = T0(data_real[kr--], 0);
3001 
3002  long ki = ptr_imag[j+1]-1;
3003  if ((ki >= ptr_imag[j]) && (ind_imag[ki] == j))
3004  ajj += T0(0, data_imag[ki--]);
3005 
3006  for (long k = ptr_real[j]; k <= kr; k++)
3007  X(j) -= data_real[k] * X(ind_real[k]);
3008 
3009  for (long k = ptr_imag[j]; k <= ki; k++)
3010  X(j) -= T0(0, data_imag[k]) * X(ind_imag[k]);
3011 
3012  X(j) *= omega / ajj;
3013  }
3014  }
3015 
3016 
3017  // Backward sweep.
3018  // (D/omega - U) X^{n+1} = (L + (1-omega)/omega D) X^{n+1/2} + B
3019  if (type_ssor % 3 == 0)
3020  for (int i = 0; i < iter; i++)
3021  {
3022  // First we compute X = (L + (1-omega)/omega D) X + B.
3023  for (int j = ma-1; j >= 0; j--)
3024  {
3025  temp = zero;
3026  ajj = zero;
3027  long kr = ptr_real[j+1]-1;
3028  if ((kr >= ptr_real[j]) && (ind_real[kr] == j))
3029  ajj = T0(data_real[kr--], 0);
3030 
3031  long ki = ptr_imag[j+1]-1;
3032  if ((ki >= ptr_imag[j]) && (ind_imag[ki] == j))
3033  ajj += T0(0, data_imag[ki--]);
3034 
3035  for (long k = ptr_real[j]; k <= kr; k++)
3036  temp -= data_real[k] * X(ind_real[k]);
3037 
3038  for (long k = ptr_imag[j]; k <= ki; k++)
3039  temp -= T0(0, data_imag[k]) * X(ind_imag[k]);
3040 
3041  X(j) = B(j) + coef * ajj * X(j) + temp;
3042  }
3043 
3044  // Then we solve (D/omega - U) X = X.
3045  for (int j = ma-1; j >= 0; j--)
3046  {
3047  temp = zero;
3048  ajj = zero;
3049  long kr = ptr_real[j+1]-1;
3050  if ((kr >= ptr_real[j]) && (ind_real[kr] == j))
3051  ajj = T0(data_real[kr--], 0);
3052 
3053  long ki = ptr_imag[j+1]-1;
3054  if ((ki >= ptr_imag[j]) && (ind_imag[ki] == j))
3055  ajj += T0(0, data_imag[ki--]);
3056 
3057  X(j) *= omega / ajj;
3058  for (long k = ptr_real[j]; k <= kr; k++)
3059  X(ind_real[k]) -= data_real[k] * X(j);
3060 
3061  for (long k = ptr_imag[j]; k <= ki; k++)
3062  X(ind_imag[k]) -= T0(0, data_imag[k]) * X(j);
3063  }
3064  }
3065  }
3066 
3067 
3069 
3076  template <class T0, class Prop0, class Allocator0,
3077  class T1, class Storage1, class Allocator1,
3078  class T2, class Storage2, class Allocator2, class T3>
3079  void SorVector(const Matrix<T0, Prop0, ArrayColSymComplexSparse, Allocator0>& A,
3080  Vector<complex<T2>, Storage2, Allocator2>& X,
3081  const Vector<complex<T1>, Storage1, Allocator1>& B,
3082  const T3& omega, int iter, int type_ssor)
3083  {
3084  complex<T1> temp, zero; T3 one;
3085  SetComplexZero(zero);
3086  SetComplexOne(one);
3087 
3088  int ma = A.GetM();
3089 
3090 #ifdef SELDON_CHECK_BOUNDS
3091  int na = A.GetN();
3092  if (na != ma)
3093  throw WrongDim("SOR", "Matrix must be squared.");
3094 
3095  if (ma != X.GetM() || ma != B.GetM())
3096  throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
3097 #endif
3098 
3099  T0 ajj;
3100 
3101  // Let us consider the following splitting : A = D - L - U
3102  // D diagonal of A
3103  // L lower part of A
3104  // U upper part of A, A is symmetric, so L = U^t
3105 
3106  // Forward sweep
3107  // (D/omega - L) X^{n+1/2} = (U + (1-omega)/omega D) X^n + B
3108  T3 coef = (T3(1) - omega) / omega;
3109  if (type_ssor % 2 == 0)
3110  for (int i = 0; i < iter; i++)
3111  {
3112  // First we do X = (U + (1-omega)/omega D) X + B
3113  for (int j = 0; j < ma; j++)
3114  {
3115  ajj = zero;
3116  int kr = A.GetRealColumnSize(j)-1;
3117  if ((kr >= 0) && (A.IndexReal(j, kr) == j))
3118  ajj = T0(A.ValueReal(j, kr--), 0);
3119 
3120  int ki = A.GetImagColumnSize(j)-1;
3121  if ((ki >= 0) && (A.IndexImag(j, ki) == j))
3122  ajj += T0(0, A.ValueImag(j, ki--));
3123 
3124 #ifdef SELDON_CHECK_BOUNDS
3125  if ( ajj == zero)
3126  throw WrongArgument("SOR", "Matrix must contain"
3127  " a non-null diagonal");
3128 #endif
3129 
3130  for (int k = 0; k <= kr; k++)
3131  X(A.IndexReal(j, k)) -= A.ValueReal(j, k) * X(j);
3132 
3133  for (int k = 0; k <= ki; k++)
3134  X(A.IndexImag(j, k)) -= T0(0, A.ValueImag(j, k)) * X(j);
3135 
3136  X(j) = coef * ajj * X(j) + B(j);
3137  }
3138 
3139  // Then we solve (D/omega - L) X = X
3140  for (int j = 0; j < ma; j++)
3141  {
3142  ajj = zero;
3143  int kr = A.GetRealColumnSize(j)-1;
3144  if ((kr >= 0) && (A.IndexReal(j, kr) == j))
3145  ajj = T0(A.ValueReal(j, kr--), 0);
3146 
3147  int ki = A.GetImagColumnSize(j)-1;
3148  if ((ki >= 0) && (A.IndexImag(j, ki) == j))
3149  ajj += T0(0, A.ValueImag(j, ki--));
3150 
3151  for (int k = 0; k <= kr; k++)
3152  X(j) -= A.ValueReal(j, k) * X(A.IndexReal(j, k));
3153 
3154  for (int k = 0; k <= ki; k++)
3155  X(j) -= T0(0, A.ValueImag(j, k)) * X(A.IndexImag(j, k));
3156 
3157  X(j) *= omega / ajj;
3158  }
3159  }
3160 
3161 
3162  // Backward sweep.
3163  // (D/omega - U) X^{n+1} = (L + (1-omega)/omega D) X^{n+1/2} + B
3164  if (type_ssor % 3 == 0)
3165  for (int i = 0; i < iter; i++)
3166  {
3167  // First we compute X = (L + (1-omega)/omega D) X + B.
3168  for (int j = ma-1; j >= 0; j--)
3169  {
3170  temp = zero;
3171  ajj = zero;
3172  int kr = A.GetRealColumnSize(j)-1;
3173  if ((kr >= 0) && (A.IndexReal(j, kr) == j))
3174  ajj = T0(A.ValueReal(j, kr--), 0);
3175 
3176  int ki = A.GetImagColumnSize(j)-1;
3177  if ((ki >= 0) && (A.IndexImag(j, ki) == j))
3178  ajj += T0(0, A.ValueImag(j, ki--));
3179 
3180  for (int k = 0; k <= kr; k++)
3181  temp -= A.ValueReal(j, k) * X(A.IndexReal(j, k));
3182 
3183  for (int k = 0; k <= ki; k++)
3184  temp -= T0(0, A.ValueImag(j, k)) * X(A.IndexImag(j, k));
3185 
3186  X(j) = B(j) + coef * ajj * X(j) + temp;
3187  }
3188 
3189  // Then we solve (D/omega - U) X = X.
3190  for (int j = ma-1; j >= 0; j--)
3191  {
3192  temp = zero;
3193  ajj = zero;
3194  int kr = A.GetRealColumnSize(j)-1;
3195  if ((kr >= 0) && (A.IndexReal(j, kr) == j))
3196  ajj = T0(A.ValueReal(j, kr--), 0);
3197 
3198  int ki = A.GetImagColumnSize(j)-1;
3199  if ((ki >= 0) && (A.IndexImag(j, ki) == j))
3200  ajj += T0(0, A.ValueImag(j, ki--));
3201 
3202  X(j) *= omega / ajj;
3203  for (int k = 0; k <= kr; k++)
3204  X(A.IndexReal(j, k)) -= A.ValueReal(j, k) * X(j);
3205 
3206  for (int k = 0; k <= ki; k++)
3207  X(A.IndexImag(j, k)) -= T0(0, A.ValueImag(j, k)) * X(j);
3208  }
3209  }
3210  }
3211 
3212 
3214 
3221  template <class T0, class Prop0, class Allocator0,
3222  class T1, class Storage1, class Allocator1,
3223  class T2, class Storage2, class Allocator2, class T3>
3224  void SorVector(const SeldonTranspose& transM,
3225  const Matrix<T0, Prop0, ColComplexSparse, Allocator0>& A,
3226  Vector<complex<T2>, Storage2, Allocator2>& X,
3227  const Vector<complex<T1>, Storage1, Allocator1>& B,
3228  const T3& omega, int iter, int type_ssor)
3229  {
3230  if (transM.NoTrans())
3231  return SorVector(A, X, B, omega, iter, type_ssor);
3232 
3233  complex<T1> temp, zero; T3 one;
3234  SetComplexZero(zero);
3235  SetComplexOne(one);
3236 
3237  int ma = A.GetM();
3238 
3239 #ifdef SELDON_CHECK_BOUNDS
3240  int na = A.GetN();
3241  if (na != ma)
3242  throw WrongDim("SOR", "Matrix must be squared.");
3243 
3244  if (ma != X.GetM() || ma != B.GetM())
3245  throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
3246 #endif
3247 
3248  long* ptr_real = A.GetRealPtr();
3249  long* ptr_imag = A.GetImagPtr();
3250  int* ind_real = A.GetRealInd();
3251  int* ind_imag = A.GetImagInd();
3252 
3253  typename ClassComplexType<T0>::Treal* data_real = A.GetRealData();
3254  typename ClassComplexType<T0>::Treal* data_imag = A.GetImagData();
3255 
3256  T0 ajj;
3257 
3258  if (type_ssor % 2 == 0)
3259  for (int i = 0; i < iter; i++)
3260  for (int j = 0; j < ma; j++)
3261  {
3262  temp = zero;
3263  ajj = zero;
3264  for (long k = ptr_real[j]; k < ptr_real[j+1]; k++)
3265  {
3266  if (ind_real[k] != j)
3267  temp += data_real[k] * X(ind_real[k]);
3268  else
3269  ajj = T0(data_real[k], 0);
3270  }
3271 
3272  for (long k = ptr_imag[j]; k < ptr_imag[j+1]; k++)
3273  {
3274  if (ind_imag[k] != j)
3275  temp += complex<T1>(0, data_imag[k]) * X(ind_imag[k]);
3276  else
3277  ajj += complex<T1>(0, data_imag[k]);
3278  }
3279 
3280 #ifdef SELDON_CHECK_BOUNDS
3281  if (ajj == zero)
3282  throw WrongArgument("SOR", "Matrix must contain"
3283  " a non-null diagonal");
3284 #endif
3285 
3286  X(j) = (one-omega) * X(j) + omega * (B(j) - temp) / ajj;
3287  }
3288 
3289  if (type_ssor % 3 == 0)
3290  for (int i = 0; i < iter; i++)
3291  for (int j = ma-1; j >= 0; j--)
3292  {
3293  temp = zero;
3294  ajj = zero;
3295  for (long k = ptr_real[j]; k < ptr_real[j+1]; k++)
3296  {
3297  if (ind_real[k] != j)
3298  temp += data_real[k] * X(ind_real[k]);
3299  else
3300  ajj = T0(data_real[k], 0);
3301  }
3302 
3303  for (long k = ptr_imag[j]; k < ptr_imag[j+1]; k++)
3304  {
3305  if (ind_imag[k] != j)
3306  temp += complex<T1>(0, data_imag[k]) * X(ind_imag[k]);
3307  else
3308  ajj += complex<T1>(0, data_imag[k]);
3309  }
3310 
3311  X(j) = (one-omega) * X(j) + omega * (B(j) - temp) / ajj;
3312  }
3313  }
3314 
3315 
3317 
3324  template <class T0, class Prop0, class Allocator0,
3325  class T1, class Storage1, class Allocator1,
3326  class T2, class Storage2, class Allocator2, class T3>
3327  void SorVector(const SeldonTranspose& transM,
3328  const Matrix<T0, Prop0, ArrayColComplexSparse, Allocator0>& A,
3329  Vector<complex<T2>, Storage2, Allocator2>& X,
3330  const Vector<complex<T1>, Storage1, Allocator1>& B,
3331  const T3& omega, int iter, int type_ssor)
3332  {
3333  if (transM.NoTrans())
3334  return SorVector(A, X, B, omega, iter, type_ssor);
3335 
3336  complex<T1> temp, zero; T3 one;
3337  SetComplexZero(zero);
3338  SetComplexOne(one);
3339 
3340  int ma = A.GetM();
3341 
3342 #ifdef SELDON_CHECK_BOUNDS
3343  int na = A.GetN();
3344  if (na != ma)
3345  throw WrongDim("SOR", "Matrix must be squared.");
3346 
3347  if (ma != X.GetM() || ma != B.GetM())
3348  throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
3349 #endif
3350 
3351  T0 ajj;
3352 
3353  if (type_ssor%2 == 0)
3354  for (int i = 0; i < iter; i++)
3355  for (int j = 0; j < ma; j++)
3356  {
3357  temp = zero;
3358  ajj = zero;
3359  for (int k = 0; k < A.GetRealColumnSize(j); k++)
3360  {
3361  if (A.IndexReal(j, k) != j)
3362  temp += A.ValueReal(j,k) * X(A.IndexReal(j,k));
3363  else
3364  ajj += A.ValueReal(j, k);
3365  }
3366 
3367  for (int k = 0; k < A.GetImagColumnSize(j); k++)
3368  {
3369  if (A.IndexImag(j, k) != j)
3370  temp += T0(0, A.ValueImag(j,k))
3371  * X(A.IndexImag(j, k));
3372  else
3373  ajj += T0(0, A.ValueImag(j,k));
3374  }
3375 
3376 #ifdef SELDON_CHECK_BOUNDS
3377  if (ajj == zero)
3378  throw WrongArgument("SOR", "Matrix must contain"
3379  " a non-null diagonal");
3380 #endif
3381 
3382  X(j) = (one-omega) * X(j) + omega * (B(j) - temp) / ajj;
3383  }
3384 
3385  if (type_ssor % 3 == 0)
3386  for (int i = 0; i < iter; i++)
3387  for (int j = ma-1; j >= 0; j--)
3388  {
3389  temp = zero;
3390  ajj = zero;
3391  for (int k = 0; k < A.GetRealColumnSize(j); k++)
3392  {
3393  if (A.IndexReal(j, k) != j)
3394  temp += A.ValueReal(j,k) * X(A.IndexReal(j,k));
3395  else
3396  ajj += A.ValueReal(j, k);
3397  }
3398 
3399  for (int k = 0; k < A.GetImagColumnSize(j); k++)
3400  {
3401  if (A.IndexImag(j, k) != j)
3402  temp += T0(0, A.ValueImag(j,k))
3403  * X(A.IndexImag(j, k));
3404  else
3405  ajj += T0(0, A.ValueImag(j,k));
3406  }
3407 
3408  X(j) = (one-omega) * X(j) + omega * (B(j) - temp) / ajj;
3409  }
3410  }
3411 
3412 
3414 
3421  template <class T0, class Prop0, class Allocator0,
3422  class T1, class Storage1, class Allocator1,
3423  class T2, class Storage2, class Allocator2, class T3>
3424  void SorVector(const SeldonTranspose& transM,
3425  const Matrix<T0, Prop0, RowComplexSparse, Allocator0>& A,
3426  Vector<complex<T2>, Storage2, Allocator2>& X,
3427  const Vector<complex<T1>, Storage1, Allocator1>& B,
3428  const T3& omega,int iter, int type_ssor)
3429  {
3430  if (transM.NoTrans())
3431  return SorVector(A, X, B, omega, iter, type_ssor);
3432 
3433  complex<T1> zero; T3 one;
3434  SetComplexZero(zero);
3435  SetComplexOne(one);
3436 
3437  int ma = A.GetM();
3438 
3439 #ifdef SELDON_CHECK_BOUNDS
3440  int na = A.GetN();
3441  if (na != ma)
3442  throw WrongDim("SOR", "Matrix must be squared.");
3443 
3444  if (ma != X.GetM() || ma != B.GetM())
3445  throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
3446 #endif
3447 
3448  long* ptr_real = A.GetRealPtr();
3449  int* ind_real = A.GetRealInd();
3450  typename ClassComplexType<T0>::Treal* data_real = A.GetRealData();
3451 
3452  long* ptr_imag = A.GetImagPtr();
3453  int* ind_imag = A.GetImagInd();
3454  typename ClassComplexType<T0>::Treal* data_imag = A.GetImagData();
3455 
3456  // Let us consider the following splitting : A = D - L - U
3457  // D diagonal of A
3458  // L lower part of A
3459  // U upper part of A
3460 
3461  // Forward sweep
3462  // (D/omega - L) X^{n+1/2} = (U + (1-omega)/omega D) X^n + B
3463  T0 ajj;
3464  T3 coef = (one - omega) / omega;
3465  if (type_ssor % 2 == 0)
3466  for (int i = 0; i < iter; i++)
3467  {
3468  // First we compute X = (U + (1-omega)/omega D) X + B
3469  for (int j = 0; j < ma; j++)
3470  {
3471  long kr = ptr_real[j];
3472  while ( (kr < ptr_real[j+1]) && (ind_real[kr] < j))
3473  {
3474  X(ind_real[kr]) -= data_real[kr]*X(j);
3475  kr++;
3476  }
3477 
3478  if ( (kr < ptr_real[j+1]) && (ind_real[kr] == j))
3479  ajj = T0(data_real[kr], 0);
3480 
3481  long ki = ptr_imag[j];
3482  while ( (ki < ptr_imag[j+1]) && (ind_imag[ki] < j))
3483  {
3484  X(ind_imag[ki]) -= T0(0, data_imag[ki])*X(j);
3485  ki++;
3486  }
3487 
3488  if ( (ki < ptr_imag[j+1]) && (ind_imag[ki] == j))
3489  ajj += T0(0, data_imag[ki]);
3490 
3491 #ifdef SELDON_CHECK_BOUNDS
3492  if (ajj == zero)
3493  throw WrongArgument("SOR", "Matrix must contain"
3494  " a non-null diagonal");
3495 #endif
3496 
3497  X(j) = B(j) + coef * ajj * X(j);
3498  }
3499 
3500  // Then we solve (D/omega - L) X = X
3501  for (int j = 0; j < ma; j++)
3502  {
3503  long kr = ptr_real[j];
3504  while ( (kr < ptr_real[j+1]) && (ind_real[kr] < j))
3505  kr++;
3506 
3507  if ( (kr < ptr_real[j+1]) && (ind_real[kr] == j))
3508  {
3509  ajj = T0(data_real[kr], 0);
3510  kr++;
3511  }
3512 
3513  long ki = ptr_imag[j];
3514  while ( (ki < ptr_imag[j+1]) && (ind_imag[ki] < j))
3515  ki++;
3516 
3517  if ( (ki < ptr_imag[j+1]) && (ind_imag[ki] == j))
3518  {
3519  ajj += T0(0, data_imag[ki]);
3520  ki++;
3521  }
3522 
3523  X(j) *= omega/ajj;
3524  while (kr < ptr_real[j+1])
3525  {
3526  X(ind_real[kr]) -= data_real[kr]*X(j);
3527  kr++;
3528  }
3529 
3530  while (ki < ptr_imag[j+1])
3531  {
3532  X(ind_imag[ki]) -= T0(0, data_imag[ki])*X(j);
3533  ki++;
3534  }
3535  }
3536  }
3537 
3538  // Backward sweep.
3539  // (D/omega - U) X^{n+1} = (L + (1-omega)/omega D) X^{n+1/2} + B
3540  if (type_ssor % 3 == 0)
3541  for (int i = 0; i < iter; i++)
3542  {
3543  // First we compute X = (L + (1-omega)/omega D) X + B.
3544  for (int j = ma-1; j >= 0; j--)
3545  {
3546  long kr = ptr_real[j+1]-1;
3547  while ( (kr >= ptr_real[j]) && (ind_real[kr] > j) )
3548  {
3549  X(ind_real[kr]) -= data_real[kr]*X(j);
3550  kr--;
3551  }
3552 
3553  if ( (kr >= ptr_real[j]) && (ind_real[kr] == j))
3554  ajj = T0(data_real[kr], 0);
3555 
3556  long ki = ptr_imag[j+1]-1;
3557  while ( (ki >= ptr_imag[j]) && (ind_imag[ki] > j) )
3558  {
3559  X(ind_imag[ki]) -= T0(0, data_imag[ki])*X(j);
3560  ki--;
3561  }
3562 
3563  if ( (ki >= ptr_imag[j]) && (ind_imag[ki] == j))
3564  ajj += T0(0, data_imag[ki]);
3565 
3566 #ifdef SELDON_CHECK_BOUNDS
3567  if (ajj == zero)
3568  throw WrongArgument("SOR", "Matrix must contain"
3569  " a non-null diagonal");
3570 #endif
3571 
3572  X(j) = B(j) + coef * ajj * X(j);
3573  }
3574 
3575  // sThen we solve (D/omega - U) X = X.
3576  for (int j = ma-1; j >= 0; j--)
3577  {
3578  long kr = ptr_real[j+1]-1;
3579  while ( (kr >= ptr_real[j]) && (ind_real[kr] > j) )
3580  kr--;
3581 
3582  if ( (kr >= ptr_real[j]) && (ind_real[kr] == j))
3583  {
3584  ajj = T0(data_real[kr], 0);
3585  kr--;
3586  }
3587 
3588  long ki = ptr_imag[j+1]-1;
3589  while ( (ki >= ptr_imag[j]) && (ind_imag[ki] > j) )
3590  ki--;
3591 
3592  if ( (ki >= ptr_imag[j]) && (ind_imag[ki] == j))
3593  {
3594  ajj += T0(0, data_imag[ki]);
3595  ki--;
3596  }
3597 
3598  X(j) *= omega/ajj;
3599  while (kr >= ptr_real[j])
3600  {
3601  X(ind_real[kr]) -= data_real[kr]*X(j);
3602  kr--;
3603  }
3604 
3605  while (ki >= ptr_imag[j])
3606  {
3607  X(ind_imag[ki]) -= T0(0, data_imag[ki])*X(j);
3608  ki--;
3609  }
3610  }
3611  }
3612  }
3613 
3614 
3616 
3623  template <class T0, class Prop0, class Allocator0,
3624  class T1, class Storage1, class Allocator1,
3625  class T2, class Storage2, class Allocator2, class T3>
3626  void SorVector(const SeldonTranspose& transM,
3627  const Matrix<T0, Prop0, ArrayRowComplexSparse, Allocator0>& A,
3628  Vector<complex<T2>, Storage2, Allocator2>& X,
3629  const Vector<complex<T1>, Storage1, Allocator1>& B,
3630  const T3& omega,int iter, int type_ssor)
3631  {
3632  if (transM.NoTrans())
3633  return SorVector(A, X, B, omega, iter, type_ssor);
3634 
3635  complex<T1> zero; T3 one;
3636  SetComplexZero(zero);
3637  SetComplexOne(one);
3638 
3639  int ma = A.GetM();
3640 
3641 #ifdef SELDON_CHECK_BOUNDS
3642  int na = A.GetN();
3643  if (na != ma)
3644  throw WrongDim("SOR", "Matrix must be squared.");
3645 
3646  if (ma != X.GetM() || ma != B.GetM())
3647  throw WrongDim("SOR", "Matrix and vector dimensions are incompatible.");
3648 #endif
3649 
3650  // Let us consider the following splitting : A = D - L - U
3651  // D diagonal of A
3652  // L lower part of A
3653  // U upper part of A
3654 
3655  // Forward sweep
3656  // (D/omega - L) X^{n+1/2} = (U + (1-omega)/omega D) X^n + B
3657  T0 ajj;
3658  T3 coef = (one - omega) / omega;
3659  if (type_ssor % 2 == 0)
3660  for (int i = 0; i < iter; i++)
3661  {
3662  // First we compute X = (U + (1-omega)/omega D) X + B
3663  for (int j = 0; j < ma; j++)
3664  {
3665  int kr = 0;
3666  while ( (kr < A.GetRealRowSize(j)) && (A.IndexReal(j, kr) < j))
3667  {
3668  X(A.IndexReal(j, kr)) -= A.ValueReal(j, kr) * X(j);
3669  kr++;
3670  }
3671 
3672  if ( (kr < A.GetRealRowSize(j)) && (A.IndexReal(j, kr) == j))
3673  ajj = T0(A.ValueReal(j, kr), 0);
3674 
3675  int ki = 0;
3676  while ( (ki < A.GetImagRowSize(j)) && (A.IndexImag(j, ki) < j))
3677  {
3678  X(A.IndexImag(j, ki)) -= T0(0, A.ValueImag(j, ki))*X(j);
3679  ki++;
3680  }
3681 
3682  if ( (ki < A.GetImagRowSize(j)) && (A.IndexImag(j, ki) == j) )
3683  ajj += T0(0, A.ValueImag(j, ki));
3684 
3685 #ifdef SELDON_CHECK_BOUNDS
3686  if (ajj == zero)
3687  throw WrongArgument("SOR", "Matrix must contain"
3688  " a non-null diagonal");
3689 #endif
3690 
3691  X(j) = B(j) + coef * ajj * X(j);
3692  }
3693 
3694  // Then we solve (D/omega - L) X = X
3695  for (int j = 0; j < ma; j++)
3696  {
3697  int kr = 0;
3698  while ( (kr < A.GetRealRowSize(j)) && (A.IndexReal(j, kr) < j))
3699  kr++;
3700 
3701  if ( (kr < A.GetRealRowSize(j)) && (A.IndexReal(j, kr) == j) )
3702  {
3703  ajj = T0(A.ValueReal(j, kr), 0);
3704  kr++;
3705  }
3706 
3707  int ki = 0;
3708  while ( (ki < A.GetImagRowSize(j)) && (A.IndexImag(j, ki) < j))
3709  ki++;
3710 
3711  if ( (ki < A.GetImagRowSize(j)) && (A.IndexImag(j, ki) == j))
3712  {
3713  ajj += T0(0, A.ValueImag(j, ki));
3714  ki++;
3715  }
3716 
3717  X(j) *= omega/ajj;
3718  while (kr < A.GetRealRowSize(j))
3719  {
3720  X(A.IndexReal(j, kr)) -= A.ValueReal(j, kr)*X(j);
3721  kr++;
3722  }
3723 
3724  while (ki < A.GetImagRowSize(j))
3725  {
3726  X(A.IndexImag(j, ki)) -= T0(0, A.ValueImag(j, ki))*X(j);
3727  ki++;
3728  }
3729  }
3730  }
3731 
3732  // Backward sweep.
3733  // (D/omega - U) X^{n+1} = (L + (1-omega)/omega D) X^{n+1/2} + B
3734  if (type_ssor % 3 == 0)
3735  for (int i = 0; i < iter; i++)
3736  {
3737  // First we compute X = (L + (1-omega)/omega D) X + B.
3738  for (int j = ma-1; j >= 0; j--)
3739  {
3740  int kr = A.GetRealRowSize(j)-1;
3741  while ( (kr >= 0) && (A.IndexReal(j, kr) > j) )
3742  {
3743  X(A.IndexReal(j, kr)) -= A.ValueReal(j, kr)*X(j);
3744  kr--;
3745  }
3746 
3747  if ( (kr >= 0) && (A.IndexReal(j, kr) == j))
3748  ajj = T0(A.ValueReal(j, kr), 0);
3749 
3750  int ki = A.GetImagRowSize(j)-1;
3751  while ( (ki >= 0) && (A.IndexImag(j, ki) > j) )
3752  {
3753  X(A.IndexImag(j, ki)) -= T0(0, A.ValueImag(j, ki))*X(j);
3754  ki--;
3755  }
3756 
3757  if ( (ki >= 0) && (A.IndexImag(j, ki) == j))
3758  ajj += T0(0, A.ValueImag(j, ki));
3759 
3760 #ifdef SELDON_CHECK_BOUNDS
3761  if (ajj == zero)
3762  throw WrongArgument("SOR", "Matrix must contain"
3763  " a non-null diagonal");
3764 #endif
3765 
3766  X(j) = B(j) + coef * ajj * X(j);
3767  }
3768 
3769  // Then we solve (D/omega - U) X = X.
3770  for (int j = ma-1; j >= 0; j--)
3771  {
3772  int kr = A.GetRealRowSize(j)-1;
3773  while ( (kr >= 0) && (A.IndexReal(j, kr) > j) )
3774  kr--;
3775 
3776  if ( (kr >= 0) && (A.IndexReal(j, kr) == j))
3777  {
3778  ajj = T0(A.ValueReal(j, kr), 0);
3779  kr--;
3780  }
3781 
3782  int ki = A.GetImagRowSize(j)-1;
3783  while ( (ki >= 0) && (A.IndexImag(j, ki) > j) )
3784  ki--;
3785 
3786  if ( (ki >= 0) && (A.IndexImag(j, ki) == j))
3787  {
3788  ajj += T0(0, A.ValueImag(j, ki));
3789  ki--;
3790  }
3791 
3792  X(j) *= omega/ajj;
3793  while (kr >= 0)
3794  {
3795  X(A.IndexReal(j, kr)) -= A.ValueReal(j, kr)*X(j);
3796  kr--;
3797  }
3798 
3799  while (ki >= 0)
3800  {
3801  X(A.IndexImag(j, ki)) -= T0(0, A.ValueImag(j, ki))*X(j);
3802  ki--;
3803  }
3804  }
3805  }
3806  }
3807 
3808 
3810 
3817  template <class T0, class Prop0, class Allocator0,
3818  class T1, class Storage1, class Allocator1,
3819  class T2, class Storage2, class Allocator2, class T3>
3820  void SorVector(const SeldonTranspose& transM,
3821  const Matrix<T0, Prop0, RowSymComplexSparse, Allocator0>& A,
3822  Vector<complex<T2>, Storage2, Allocator2>& X,
3823  const Vector<complex<T1>, Storage1, Allocator1>& B,
3824  const T3& omega, int iter, int type_ssor)
3825  {
3826  SorVector(A, X, B, omega, iter, type_ssor);
3827  }
3828 
3829 
3831 
3838  template <class T0, class Prop0, class Allocator0,
3839  class T1, class Storage1, class Allocator1,
3840  class T2, class Storage2, class Allocator2, class T3>
3841  void SorVector(const SeldonTranspose& transM,
3842  const Matrix<T0, Prop0, ArrayRowSymComplexSparse, Allocator0>& A,
3843  Vector<complex<T2>, Storage2, Allocator2>& X,
3844  const Vector<complex<T1>, Storage1, Allocator1>& B,
3845  const T3& omega, int iter, int type_ssor)
3846  {
3847  SorVector(A, X, B, omega, iter, type_ssor);
3848  }
3849 
3850 
3852 
3859  template <class T0, class Prop0, class Allocator0,
3860  class T1, class Storage1, class Allocator1,
3861  class T2, class Storage2, class Allocator2, class T3>
3862  void SorVector(const SeldonTranspose& transM,
3863  const Matrix<T0, Prop0, ColSymComplexSparse, Allocator0>& A,
3864  Vector<complex<T2>, Storage2, Allocator2>& X,
3865  const Vector<complex<T1>, Storage1, Allocator1>& B,
3866  const T3& omega, int iter, int type_ssor)
3867  {
3868  SorVector(A, X, B, omega, iter, type_ssor);
3869  }
3870 
3871 
3873 
3880  template <class T0, class Prop0, class Allocator0,
3881  class T1, class Storage1, class Allocator1,
3882  class T2, class Storage2, class Allocator2, class T3>
3883  void SorVector(const SeldonTranspose& transM,
3884  const Matrix<T0, Prop0, ArrayColSymComplexSparse, Allocator0>& A,
3885  Vector<complex<T2>, Storage2, Allocator2>& X,
3886  const Vector<complex<T1>, Storage1, Allocator1>& B,
3887  const T3& omega, int iter, int type_ssor)
3888  {
3889  SorVector(A, X, B, omega, iter, type_ssor);
3890  }
3891 
3892 
3893  // SOR //
3895 
3896 }
3897 
3898 #define SELDON_FILE_FUNCTIONS_MATVECT_COMPLEX_CXX
3899 #endif
3900 
Seldon::SeldonTranspose
Definition: MatrixFlag.hxx:32
Seldon::Vector
Definition: SeldonHeader.hxx:207
Seldon::Matrix
Definition: SeldonHeader.hxx:226
Seldon::CheckDim
void CheckDim(const Matrix< T0, Prop0, Storage0, Allocator0 > &A, const Matrix< T1, Prop1, Storage1, Allocator1 > &B, const Matrix< T2, Prop2, Storage2, Allocator2 > &C, string function)
Checks the compatibility of the dimensions.
Definition: Functions_Matrix.cxx:2132
Seldon::WrongArgument
Definition: Errors.hxx:76
Seldon::SorVector
void SorVector(const Matrix< T0, Prop0, Storage0, Allocator0 > &M, Vector< T2, Storage2, Allocator2 > &Y, const Vector< T1, Storage1, Allocator1 > &X, const T3 &omega, int iter, int type_ssor)
Solve M Y = X with S.O.R. method.
Definition: Functions_MatVect.cxx:1638
Seldon::WrongDim
Definition: Errors.hxx:102
Seldon
Seldon namespace.
Definition: Array.cxx:24
Seldon::SetComplexReal
void SetComplexReal(int n, std::complex< T > &number)
Sets a complex number to (n, 0).
Definition: CommonInline.cxx:234