Functions_BaseInline.cxx
1 #ifndef SELDON_FILE_FUNCTIONS_BASE_INLINE_CXX
2 
3 // in this file, we put functions such as Add, Mlt, MltAdd
4 // with a reduced number of templates in order
5 // to forbid undesired instantiations of generic functions
6 
7 /*
8  Functions defined in this file:
9 
10  M X -> Y
11  Mlt(M, X, Y)
12 
13  alpha M X -> Y
14  Mlt(alpha, M, X, Y)
15 
16  M X -> Y
17  M^T X -> Y
18  Mlt(Trans, M, X, Y)
19 
20 */
21 
22 namespace Seldon
23 {
24 
25  /************************
26  * Functions in Array3D *
27  ************************/
28 
29  template<class T, class Allocator>
30  inline void Mlt(const T& alpha, Array3D<T, Allocator>& A)
31  {
32  MltScalar(alpha, A);
33  }
34 
35  template<class T, class Allocator>
36  inline void Mlt(const T& alpha, Array3D<complex<T>, Allocator>& A)
37  {
38  MltScalar(alpha, A);
39  }
40 
41  /*********************************
42  * Functions in Functions_Vector *
43  *********************************/
44 
45 #ifdef SELDON_WITH_REDUCED_TEMPLATE
46  template<class T, class Storage, class Allocator>
47  inline void Mlt(const T& alpha, Vector<T, Storage, Allocator>& A)
48  {
49  MltScalar(alpha, A);
50  }
51 #else
52  template<class T, class T0, class Storage, class Allocator>
53  inline void Mlt(const T& alpha, Vector<T0, Storage, Allocator>& A)
54  {
55  MltScalar(alpha, A);
56  }
57 #endif
58 
59  template<class T, class Storage, class Allocator>
60  inline void Mlt(const T& alpha, Vector<complex<T>, Storage, Allocator>& A)
61  {
62  MltScalar(alpha, A);
63  }
64 
65  template<class T, class Storage1, class Allocator1,
66  class Storage2, class Allocator2>
67  inline void Add(const T& alpha, const Vector<T, Storage1, Allocator1>& X,
68  Vector<T, Storage2, Allocator2>& Y)
69  {
70  AddVector(alpha, X, Y);
71  }
72 
73  template<class T, class Storage1, class Allocator1,
74  class Storage2, class Allocator2>
75  inline void Add(const T& alpha, const Vector<complex<T>, Storage1, Allocator1>& X,
76  Vector<complex<T>, Storage2, Allocator2>& Y)
77  {
78  AddVector(alpha, X, Y);
79  }
80 
81  template<class T, class Storage1, class Allocator1,
82  class Storage2, class Allocator2>
83  inline void Add(const T& alpha, const Vector<T, Storage1, Allocator1>& X,
84  const T& beta, Vector<T, Storage2, Allocator2>& Y)
85  {
86  AddVector(alpha, X, beta, Y);
87  }
88 
89  template<class T, class Storage1, class Allocator1,
90  class Storage2, class Allocator2>
91  inline void Copy(const Vector<T, Storage1, Allocator1>& X,
92  Vector<T, Storage2, Allocator2>& Y)
93  {
94  CopyVector(X, Y);
95  }
96 
97  template<class T, class Storage1, class Allocator1,
98  class Storage2, class Allocator2>
99  inline T DotProd(const Vector<T, Storage1, Allocator1>& X,
100  const Vector<T, Storage2, Allocator2>& Y)
101  {
102  return DotProdVector(X, Y);
103  }
104 
105  template<class T, class Storage1, class Allocator1,
106  class Storage2, class Allocator2>
107  inline T DotProdConj(const Vector<T, Storage1, Allocator1>& X,
108  const Vector<T, Storage2, Allocator2>& Y)
109  {
110  return DotProdVector(X, Y);
111  }
112 
113  template<class T, class Storage1, class Allocator1,
114  class Storage2, class Allocator2>
115  inline complex<T> DotProdConj(const Vector<complex<T>, Storage1, Allocator1>& X,
116  const Vector<complex<T>, Storage2, Allocator2>& Y)
117  {
118  return DotProdConjVector(X, Y);
119  }
120 
121  template<class T, class Storage1, class Allocator1,
122  class Storage2, class Allocator2>
123  inline complex<T> DotProdConj(const Vector<complex<T>, Storage1, Allocator1>& X,
124  const Vector<T, Storage2, Allocator2>& Y)
125  {
126  return DotProdConjVector(X, Y);
127  }
128 
129  template<class T, class Storage1, class Allocator1,
130  class Storage2, class Allocator2>
131  inline complex<T> DotProdConj(const Vector<T, Storage1, Allocator1>& X,
132  const Vector<complex<T>, Storage2, Allocator2>& Y)
133  {
134  return DotProdVector(X, Y);
135  }
136 
138  template<class T, class E>
139  inline typename ClassComplexType<T>::Treal
140  Norm1(const VectorExpression<T, E>& X)
141  {
142  typename ClassComplexType<T>::Treal value(0);
143 
144  for (long i = 0; i < X.GetSize(); i++)
145  value += ComplexAbs(X(i));
146 
147  return value;
148  }
149 
151  template<class T, class E>
152  typename ClassComplexType<T>::Treal
153  Norm2(const VectorExpression<T, E>& X)
154  {
155  typename ClassComplexType<T>::Treal value(0);
156 
157  for (long i = 0; i < X.GetSize(); i++)
158  value += absSquare(X(i));
159 
160  return sqrt(value);
161  }
162 
163 
164  /**********************************
165  * Functions in Functions_MatVect *
166  **********************************/
167 
168 #ifdef SELDON_WITH_REDUCED_TEMPLATE
169  template <class T, class Prop0, class Storage0, class Allocator0,
170  class Storage1, class Allocator1,
171  class Storage2, class Allocator2>
172  inline void Mlt(const Matrix<T, Prop0, Storage0, Allocator0>& M,
173  const Vector<T, Storage1, Allocator1>& X,
174  Vector<T, Storage2, Allocator2>& Y)
175  {
176  MltVector(M, X, Y);
177  }
178 #else
179  template <class T0, class Prop0, class Storage0, class Allocator0,
180  class T1, class Storage1, class Allocator1,
181  class T2, class Storage2, class Allocator2>
182  inline void Mlt(const Matrix<T0, Prop0, Storage0, Allocator0>& M,
183  const Vector<T1, Storage1, Allocator1>& X,
184  Vector<T2, Storage2, Allocator2>& Y)
185  {
186  MltVector(M, X, Y);
187  }
188 #endif
189 
190  template <class T, class Prop0, class Storage0, class Allocator0,
191  class Storage1, class Allocator1,
192  class Storage2, class Allocator2>
193  inline void Mlt(const Matrix<complex<T>, Prop0, Storage0, Allocator0>& M,
194  const Vector<T, Storage1, Allocator1>& X,
195  Vector<T, Storage2, Allocator2>& Y)
196  {
197  throw WrongArgument("Mlt", "Incompatible matrix-vector product");
198  }
199 
200  template <class T, class Prop0, class Storage0, class Allocator0,
201  class Storage1, class Allocator1,
202  class Storage2, class Allocator2>
203  inline void Mlt(const Matrix<T, Prop0, Storage0, Allocator0>& M,
204  const Vector<complex<T>, Storage1, Allocator1>& X,
205  Vector<complex<T>, Storage2, Allocator2>& Y)
206  {
207  MltVector(M, X, Y);
208  }
209 
210 
211  template<class T, class Prop1, class Storage1, class Allocator1,
212  class Storage2, class Allocator2,
213  class Storage3, class Allocator3>
214  void Mlt(const SeldonTranspose& Trans,
215  const Matrix<T, Prop1, Storage1, Allocator1>& M,
216  const Vector<T, Storage2, Allocator2>& X,
217  Vector<T, Storage3, Allocator3>& Y)
218  {
219  MltVector(Trans, M, X, Y);
220  }
221 
222 
223  template<class T, class Prop1, class Storage1, class Allocator1,
224  class Storage2, class Allocator2,
225  class Storage3, class Allocator3>
226  inline void Mlt(const SeldonTranspose& Trans,
227  const Matrix<T, Prop1, Storage1, Allocator1>& M,
228  const Vector<complex<T>, Storage2, Allocator2>& X,
229  Vector<complex<T>, Storage3, Allocator3>& Y)
230  {
231  MltVector(Trans, M, X, Y);
232  }
233 
234  template<class T, class Prop1, class Storage1, class Allocator1,
235  class Storage2, class Allocator2,
236  class Storage3, class Allocator3>
237  inline void Mlt(const SeldonTranspose& Trans,
238  const Matrix<complex<T>, Prop1, Storage1, Allocator1>& M,
239  const Vector<T, Storage2, Allocator2>& X,
240  Vector<T, Storage3, Allocator3>& Y)
241  {
242  throw WrongArgument("Mlt", "Incompatible matrix-vector product");
243  }
244 
245 #ifdef SELDON_WITH_REDUCED_TEMPLATE
246  template <class T, class Prop1, class Storage1, class Allocator1,
247  class Storage2, class Allocator2,
248  class Storage3, class Allocator3>
249  inline void Mlt(const T& alpha,
250  const Matrix<T, Prop1, Storage1, Allocator1>& M,
251  const Vector<T, Storage2, Allocator2>& X,
252  Vector<T, Storage3, Allocator3>& Y)
253  {
254  MltVector(M, X, Y);
255  Mlt(alpha, Y);
256  }
257 #else
258  template <class T0,
259  class T1, class Prop1, class Storage1, class Allocator1,
260  class T2, class Storage2, class Allocator2,
261  class T3, class Storage3, class Allocator3>
262  inline void Mlt(const T0& alpha,
263  const Matrix<T1, Prop1, Storage1, Allocator1>& M,
264  const Vector<T2, Storage2, Allocator2>& X,
265  Vector<T3, Storage3, Allocator3>& Y)
266  {
267  MltVector(M, X, Y);
268  Mlt(alpha, Y);
269  }
270 #endif
271 
272  template <class T, class Prop1, class Storage1, class Allocator1,
273  class Storage2, class Allocator2,
274  class Storage3, class Allocator3>
275  inline void Mlt(const complex<T>& alpha,
276  const Matrix<T, Prop1, Storage1, Allocator1>& M,
277  const Vector<complex<T>, Storage2, Allocator2>& X,
278  Vector<complex<T>, Storage3, Allocator3>& Y)
279  {
280  MltVector(M, X, Y);
281  Mlt(alpha, Y);
282  }
283 
284  template <class T, class Prop1, class Storage1, class Allocator1,
285  class Storage2, class Allocator2,
286  class Storage3, class Allocator3>
287  inline void Mlt(const T& alpha,
288  const Matrix<complex<T>, Prop1, Storage1, Allocator1>& M,
289  const Vector<T, Storage2, Allocator2>& X,
290  Vector<T, Storage3, Allocator3>& Y)
291  {
292  throw WrongArgument("Mlt", "Incompatible matrix-vector product");
293  }
294 
295 #ifdef SELDON_WITH_REDUCED_TEMPLATE
296  template<class T,
297  class Prop1, class Storage1, class Allocator1,
298  class Storage2, class Allocator2,
299  class Storage4, class Allocator4>
300  inline void MltAdd(const T& alpha,
301  const Matrix<T, Prop1, Storage1, Allocator1>& M,
302  const Vector<T, Storage2, Allocator2>& X,
303  const T& beta,
304  Vector<T, Storage4, Allocator4>& Y)
305  {
306  MltAddVector(alpha, M, X, beta, Y);
307  }
308 #else
309  template<class T,
310  class T1, class Prop1, class Storage1, class Allocator1,
311  class T2, class Storage2, class Allocator2,
312  class T3,
313  class T4, class Storage4, class Allocator4>
314  inline void MltAdd(const T& alpha,
315  const Matrix<T1, Prop1, Storage1, Allocator1>& M,
316  const Vector<T2, Storage2, Allocator2>& X,
317  const T3& beta,
318  Vector<T4, Storage4, Allocator4>& Y)
319  {
320  MltAddVector(alpha, M, X, beta, Y);
321  }
322 #endif
323 
324  template<class T,
325  class Prop1, class Storage1, class Allocator1,
326  class Storage2, class Allocator2,
327  class Storage4, class Allocator4>
328  inline void MltAdd(const T& alpha,
329  const Matrix<complex<T>, Prop1, Storage1, Allocator1>& M,
330  const Vector<complex<T>, Storage2, Allocator2>& X,
331  const T& beta,
332  Vector<complex<T>, Storage4, Allocator4>& Y)
333  {
334  MltAddVector(alpha, M, X, beta, Y);
335  }
336 
337  template<class T,
338  class Prop1, class Storage1, class Allocator1,
339  class Storage2, class Allocator2,
340  class Storage4, class Allocator4>
341  inline void MltAdd(const T& alpha,
342  const Matrix<T, Prop1, Storage1, Allocator1>& M,
343  const Vector<complex<T>, Storage2, Allocator2>& X,
344  const T& beta,
345  Vector<complex<T>, Storage4, Allocator4>& Y)
346  {
347  MltAddVector(alpha, M, X, beta, Y);
348  }
349 
350  template<class T,
351  class Prop1, class Storage1, class Allocator1,
352  class Storage2, class Allocator2,
353  class Storage4, class Allocator4>
354  inline void MltAdd(const T& alpha,
355  const Matrix<complex<T>, Prop1, Storage1, Allocator1>& M,
356  const Vector<T, Storage2, Allocator2>& X,
357  const T& beta,
358  Vector<T, Storage4, Allocator4>& Y)
359  {
360  throw WrongArgument("MltAdd", "Incompatible matrix-vector product");
361  }
362 
363  template<class T,
364  class Prop1, class Storage1, class Allocator1,
365  class Storage2, class Allocator2,
366  class Storage4, class Allocator4>
367  inline void MltAdd(const complex<T>& alpha,
368  const Matrix<T, Prop1, Storage1, Allocator1>& M,
369  const Vector<complex<T>, Storage2, Allocator2>& X,
370  const complex<T>& beta,
371  Vector<complex<T>, Storage4, Allocator4>& Y)
372  {
373  MltAddVector(alpha, M, X, beta, Y);
374  }
375 
376 
377  template<class T,
378  class Prop1, class Storage1, class Allocator1,
379  class Storage2, class Allocator2,
380  class Storage4, class Allocator4>
381  inline void MltAdd(const T& alpha, const SeldonTranspose& Trans,
382  const Matrix<T, Prop1, Storage1, Allocator1>& M,
383  const Vector<T, Storage2, Allocator2>& X,
384  const T& beta,
385  Vector<T, Storage4, Allocator4>& Y)
386  {
387  MltAddVector(alpha, Trans, M, X, beta, Y);
388  }
389 
390 
391  template<class T,
392  class Prop1, class Storage1, class Allocator1,
393  class Storage2, class Allocator2,
394  class Storage4, class Allocator4>
395  inline void MltAdd(const T& alpha, const SeldonTranspose& Trans,
396  const Matrix<complex<T>, Prop1, Storage1, Allocator1>& M,
397  const Vector<complex<T>, Storage2, Allocator2>& X,
398  const T& beta,
399  Vector<complex<T>, Storage4, Allocator4>& Y)
400  {
401  MltAddVector(alpha, Trans, M, X, beta, Y);
402  }
403 
404  template<class T,
405  class Prop1, class Storage1, class Allocator1,
406  class Storage2, class Allocator2,
407  class Storage4, class Allocator4>
408  inline void MltAdd(const T& alpha, const SeldonTranspose& Trans,
409  const Matrix<T, Prop1, Storage1, Allocator1>& M,
410  const Vector<complex<T>, Storage2, Allocator2>& X,
411  const T& beta,
412  Vector<complex<T>, Storage4, Allocator4>& Y)
413  {
414  MltAddVector(alpha, Trans, M, X, beta, Y);
415  }
416 
417  template<class T,
418  class Prop1, class Storage1, class Allocator1,
419  class Storage2, class Allocator2,
420  class Storage4, class Allocator4>
421  inline void MltAdd(const T& alpha, const SeldonTranspose& Trans,
422  const Matrix<complex<T>, Prop1, Storage1, Allocator1>& M,
423  const Vector<T, Storage2, Allocator2>& X,
424  const T& beta,
425  Vector<T, Storage4, Allocator4>& Y)
426  {
427  throw WrongArgument("MltAdd", "Incompatible matrix-vector product");
428  }
429 
430  template<class T,
431  class Prop1, class Storage1, class Allocator1,
432  class Storage2, class Allocator2,
433  class Storage4, class Allocator4>
434  inline void MltAdd(const complex<T>& alpha, const SeldonTranspose& Trans,
435  const Matrix<T, Prop1, Storage1, Allocator1>& M,
436  const Vector<complex<T>, Storage2, Allocator2>& X,
437  const complex<T>& beta,
438  Vector<complex<T>, Storage4, Allocator4>& Y)
439  {
440 
441  MltAddVector(alpha, Trans, M, X, beta, Y);
442  }
443 
444  template <class T, class Prop0, class Storage0, class Allocator0,
445  class Storage1, class Allocator1,
446  class Storage2, class Allocator2>
447  inline void SOR(const Matrix<T, Prop0, Storage0, Allocator0>& M,
448  Vector<T, Storage2, Allocator2>& Y,
449  const Vector<T, Storage1, Allocator1>& X,
450  const typename ClassComplexType<T>::Treal& omega,
451  int iter, int type_ssor)
452  {
453  SorVector(M, Y, X, omega, iter, type_ssor);
454  }
455 
456 
457  template <class T, class Prop0, class Storage0, class Allocator0,
458  class Storage1, class Allocator1,
459  class Storage2, class Allocator2>
460  inline void SOR(const SeldonTranspose& trans,
461  const Matrix<T, Prop0, Storage0, Allocator0>& M,
462  Vector<T, Storage2, Allocator2>& Y,
463  const Vector<T, Storage1, Allocator1>& X,
464  const typename ClassComplexType<T>::Treal& omega, int iter, int type_ssor)
465  {
466  SorVector(trans, M, Y, X, omega, iter, type_ssor);
467  }
468 
469 
470  template <class T, class Prop0, class Storage0, class Allocator0,
471  class Storage1, class Allocator1,
472  class Storage2, class Allocator2>
473  inline void SOR(const Matrix<T, Prop0, Storage0, Allocator0>& M,
474  Vector<complex<T>, Storage2, Allocator2>& Y,
475  const Vector<complex<T>, Storage1, Allocator1>& X,
476  const typename ClassComplexType<T>::Treal& omega,
477  int iter, int type_ssor)
478  {
479  SorVector(M, Y, X, omega, iter, type_ssor);
480  }
481 
482 
483  template <class T, class Prop0, class Storage0, class Allocator0,
484  class Storage1, class Allocator1,
485  class Storage2, class Allocator2>
486  inline void SOR(const SeldonTranspose& trans,
487  const Matrix<T, Prop0, Storage0, Allocator0>& M,
488  Vector<complex<T>, Storage2, Allocator2>& Y,
489  const Vector<complex<T>, Storage1, Allocator1>& X,
490  const typename ClassComplexType<T>::Treal& omega, int iter, int type_ssor)
491  {
492  SorVector(trans, M, Y, X, omega, iter, type_ssor);
493  }
494 
495 
496  template <class T, class Prop0, class Storage0, class Allocator0,
497  class Storage1, class Allocator1,
498  class Storage2, class Allocator2>
499  inline void SOR(const Matrix<complex<T>, Prop0, Storage0, Allocator0>& M,
500  Vector<T, Storage2, Allocator2>& Y,
501  const Vector<T, Storage1, Allocator1>& X,
502  const typename ClassComplexType<T>::Treal& omega,
503  int iter, int type_ssor)
504  {
505  throw WrongArgument("SOR", "incompatible types");
506  }
507 
508 
509  template <class T, class Prop0, class Storage0, class Allocator0,
510  class Storage1, class Allocator1,
511  class Storage2, class Allocator2>
512  inline void SOR(const SeldonTranspose&,
513  const Matrix<complex<T>, Prop0, Storage0, Allocator0>& M,
514  Vector<T, Storage2, Allocator2>& Y,
515  const Vector<T, Storage1, Allocator1>& X,
516  const typename ClassComplexType<T>::Treal& omega, int iter, int type_ssor)
517  {
518  throw WrongArgument("SOR", "incompatible types");
519  }
520 
521 
522  template <class T, class Prop0, class Allocator0,
523  class Allocator1>
524  void
525  inline Solve(const SeldonUplo& Uplo,
526  const SeldonTranspose& TransA,
527  const SeldonDiag& DiagA,
528  const Matrix<T, Prop0, RowSparse, Allocator0>& A,
529  const Vector<T, VectFull, Allocator1>& X,
530  Vector<T, VectFull, Allocator1>& Y)
531  {
532  SolveTriangular(Uplo, TransA, DiagA, A, X, Y);
533  }
534 
535  // SolveLU without pivot
536  template <class T, class Prop0, class Storage0, class Allocator0,
537  class Storage1, class Allocator1>
538  inline void SolveLU(const Matrix<T, Prop0, Storage0, Allocator0>& M,
539  Vector<T, Storage1, Allocator1>& Y)
540  {
541  SolveLuVector(M, Y);
542  }
543 
544  template <class T, class Prop0, class Storage0, class Allocator0,
545  class Storage1, class Allocator1>
546  inline void SolveLU(const Matrix<complex<T>, Prop0, Storage0, Allocator0>& M,
547  Vector<T, Storage1, Allocator1>& Y)
548  {
549  throw WrongArgument("SolveLU", "incompatible types");
550  }
551 
552  template <class T, class Prop0, class Storage0, class Allocator0,
553  class Storage1, class Allocator1>
554  inline void SolveLU(const Matrix<T, Prop0, Storage0, Allocator0>& M,
555  Vector<complex<T>, Storage1, Allocator1>& Y)
556  {
557  SolveLuVector(M, Y);
558  }
559 
560  template <class T, class Prop0, class Storage0, class Allocator0,
561  class Storage1, class Allocator1>
562  inline void SolveLU(const SeldonTranspose& transA,
563  const Matrix<T, Prop0, Storage0, Allocator0>& M,
564  Vector<T, Storage1, Allocator1>& Y)
565  {
566  SolveLuVector(transA, M, Y);
567  }
568 
569  template <class T, class Prop0, class Storage0, class Allocator0,
570  class Storage1, class Allocator1>
571  inline void SolveLU(const SeldonTranspose& transA,
572  const Matrix<complex<T>, Prop0, Storage0, Allocator0>& M,
573  Vector<T, Storage1, Allocator1>& Y)
574  {
575  throw WrongArgument("SolveLU", "incompatible types");
576  }
577 
578  template <class T, class Prop0, class Storage0, class Allocator0,
579  class Storage1, class Allocator1>
580  inline void SolveLU(const SeldonTranspose& transA,
581  const Matrix<T, Prop0, Storage0, Allocator0>& M,
582  Vector<complex<T>, Storage1, Allocator1>& Y)
583  {
584  SolveLuVector(transA, M, Y);
585  }
586 
587  // SolveLU with pivot
588 
589  template<class T0, class Prop0, class Storage0, class Allocator0>
590  inline void SolveLU(const Matrix<T0, Prop0, Storage0, Allocator0>& A,
591  const Vector<int>& pivot, Vector<T0>& x)
592  {
593  SolveLuVector(A, pivot, x);
594  }
595 
596  template<class T0, class Prop0, class Storage0, class Allocator0>
597  inline void SolveLU(const Matrix<complex<T0>, Prop0, Storage0, Allocator0>& A,
598  const Vector<int>& pivot, Vector<T0>& x)
599  {
600  throw WrongArgument("SolveLU", "incompatible types");
601  }
602 
603  template<class T0, class Storage0, class Allocator0>
604  inline void SolveLU(const SeldonTranspose& trans,
605  const Matrix<T0, General, Storage0, Allocator0>& A,
606  const Vector<int>& pivot, Vector<T0>& x)
607  {
608  SolveLuVector(trans, A, pivot, x);
609  }
610 
611  template<class T0, class Storage0, class Allocator0>
612  inline void SolveLU(const SeldonTranspose& trans,
613  const Matrix<complex<T0>, General, Storage0, Allocator0>& A,
614  const Vector<int>& pivot, Vector<T0>& x)
615  {
616  throw WrongArgument("SolveLU", "incompatible types");
617  }
618 
619 
620  template<class T0, class Storage0, class Allocator0>
621  inline void SolveLU(const SeldonTranspose& trans,
622  const Matrix<T0, Symmetric, Storage0, Allocator0>& A,
623  const Vector<int>& pivot, Vector<T0>& x)
624  {
625  SolveLuVector(A, pivot, x);
626  }
627 
628  template<class T0, class Storage0, class Allocator0>
629  inline void SolveLU(const SeldonTranspose& trans,
630  const Matrix<complex<T0>, Symmetric, Storage0, Allocator0>& A,
631  const Vector<int>& pivot, Vector<T0>& x)
632  {
633  throw WrongArgument("SolveLU", "incompatible types");
634  }
635 
636  /*************************************
637  * Functions in Functions_Matrix.cxx *
638  *************************************/
639 
640  template<class T, class Prop, class Storage, class Allocator>
641  inline void Mlt(const T& alpha, Matrix<T, Prop, Storage, Allocator>& A)
642  {
643  MltScalar(alpha, A);
644  }
645 
646  template<class T, class Prop, class Storage, class Allocator>
647  inline void Mlt(const T& alpha, Matrix<complex<T>, Prop, Storage, Allocator>& A)
648  {
649  MltScalar(alpha, A);
650  }
651 
652  template <class T,
653  class Prop1, class Storage1, class Allocator1,
654  class Prop2, class Storage2, class Allocator2,
655  class Prop3, class Storage3, class Allocator3>
656  inline void Mlt(const T& alpha,
657  const Matrix<T, Prop1, Storage1, Allocator1>& A,
658  const Matrix<T, Prop2, Storage2, Allocator2>& B,
659  Matrix<T, Prop3, Storage3, Allocator3>& C)
660  {
661  T zero;
662  SetComplexZero(zero);
663  C.Zero();
664  MltAddMatrix(alpha, A, B, zero, C);
665  }
666 
667  template <class T, class Prop0, class Storage0, class Allocator0,
668  class Prop1, class Storage1, class Allocator1,
669  class Prop2, class Storage2, class Allocator2>
670  inline void Mlt(const Matrix<T, Prop0, Storage0, Allocator0>& A,
671  const Matrix<T, Prop1, Storage1, Allocator1>& B,
672  Matrix<T, Prop2, Storage2, Allocator2>& C)
673  {
674  T one, zero;
675  SetComplexZero(zero);
676  SetComplexOne(one);
677  C.Zero();
678  MltAddMatrix(one, A, B, zero, C);
679  }
680 
681  template <class T,
682  class Prop1, class Storage1, class Allocator1,
683  class Prop2, class Storage2, class Allocator2,
684  class Prop3, class Storage3, class Allocator3>
685  inline void MltAdd(const T& alpha,
686  const Matrix<T, Prop1, Storage1, Allocator1>& A,
687  const Matrix<T, Prop2, Storage2, Allocator2>& B,
688  const T& beta,
689  Matrix<T, Prop3, Storage3, Allocator3>& C)
690  {
691  MltAddMatrix(alpha, A, B, beta, C);
692  }
693 
694  template <class T,
695  class Prop1, class Storage1, class Allocator1,
696  class Prop2, class Storage2, class Allocator2,
697  class Prop3, class Storage3, class Allocator3>
698  inline void MltAdd(const T& alpha, const SeldonTranspose& transA,
699  const Matrix<T, Prop1, Storage1, Allocator1>& A,
700  const SeldonTranspose& transB,
701  const Matrix<T, Prop2, Storage2, Allocator2>& B,
702  const T& beta,
703  Matrix<T, Prop3, Storage3, Allocator3>& C)
704  {
705  MltAddMatrix(alpha, transA, A, transB, B, beta, C);
706  }
707 
708 
709  template <class T,
710  class Prop1, class Storage1, class Allocator1,
711  class Prop2, class Storage2, class Allocator2>
712  inline void Add(const T& alpha,
713  const Matrix<T, Prop1, Storage1, Allocator1>& A,
714  Matrix<T, Prop2, Storage2, Allocator2>& B)
715  {
716  AddMatrix(alpha, A, B);
717  }
718 
719  template <class T,
720  class Prop1, class Storage1, class Allocator1,
721  class Prop2, class Storage2, class Allocator2>
722  inline void Add(const complex<T>& alpha,
723  const Matrix<T, Prop1, Storage1, Allocator1>& A,
724  Matrix<complex<T>, Prop2, Storage2, Allocator2>& B)
725  {
726  AddMatrix(alpha, A, B);
727  }
728 
729 
730  template <class T,
731  class Prop1, class Storage1, class Allocator1,
732  class Prop2, class Storage2, class Allocator2>
733  inline void Add(const T& alpha,
734  const Matrix<complex<T>, Prop1, Storage1, Allocator1>& A,
735  Matrix<complex<T>, Prop2, Storage2, Allocator2>& B)
736  {
737  AddMatrix(alpha, A, B);
738  }
739 
740  template <class T,
741  class Prop1, class Storage1, class Allocator1,
742  class Prop2, class Storage2, class Allocator2>
743  inline void Add(const T& alpha,
744  const Matrix<complex<T>, Prop1, Storage1, Allocator1>& A,
745  Matrix<T, Prop2, Storage2, Allocator2>& B)
746  {
747  throw WrongArgument("Add", "incompatible types");
748  }
749 
750 
751  template<class T, class Prop1, class Storage1, class Allocator1,
752  class Prop2, class Storage2, class Allocator2>
753  inline void Copy(const Matrix<T, Prop1, Storage1, Allocator1>& A,
754  Matrix<T, Prop2, Storage2, Allocator2>& B)
755  {
756  CopyMatrix(A, B);
757  }
758 
759  template<class T, class Prop1, class Storage1, class Allocator1,
760  class Prop2, class Storage2, class Allocator2>
761  inline void Copy(const Matrix<complex<T>, Prop1, Storage1, Allocator1>& A,
762  Matrix<T, Prop2, Storage2, Allocator2>& B)
763  {
764  throw WrongArgument("Copy", "incompatible types");
765  }
766 
767  template<class T, class Prop1, class Storage1, class Allocator1,
768  class Prop2, class Storage2, class Allocator2>
769  inline void Copy(const Matrix<T, Prop1, Storage1, Allocator1>& A,
770  Matrix<complex<T>, Prop2, Storage2, Allocator2>& B)
771  {
772  CopyMatrix(A, B);
773  }
774 
775 
777  template<class T, class Prop, class Storage, class Allocator>
778  inline bool IsSymmetricMatrix(const Matrix<T, Prop, Storage, Allocator>& A)
779  {
780  return false;
781  }
782 
783 
785  template<class T, class Storage, class Allocator>
786  inline bool IsSymmetricMatrix(const Matrix<T, Symmetric, Storage, Allocator>& A)
787  {
788  return true;
789  }
790 
791 
793  template<class T, class Prop, class Storage, class Allocator>
794  inline bool IsComplexMatrix(const Matrix<T, Prop, Storage, Allocator>& A)
795  {
796  return false;
797  }
798 
799 
801  template<class T, class Prop, class Storage, class Allocator>
802  inline bool IsComplexMatrix(const Matrix<complex<T>, Prop, Storage, Allocator>& A)
803  {
804  return true;
805  }
806 
807 }
808 
809 #define SELDON_FILE_FUNCTIONS_BASE_INLINE_CXX
810 #endif
Seldon::VectorExpression::GetSize
long GetSize() const
returns the size of the associated vector
Definition: VectorExpressionInline.cxx:8
Seldon::MltScalar
void MltScalar(const T0 &alpha, Array3D< T, Allocator > &A)
Multiplication of all elements of a 3D array by a scalar.
Definition: Array3D.cxx:539
Seldon::Norm2
ClassComplexType< T >::Treal Norm2(const VectorExpression< T, E > &X)
returns 2-norm of an expression with vectors
Definition: Functions_BaseInline.cxx:153
Seldon::Norm1
ClassComplexType< T >::Treal Norm1(const VectorExpression< T, E > &X)
returns 1-norm of an expression with vectors
Definition: Functions_BaseInline.cxx:140
Seldon::IsComplexMatrix
bool IsComplexMatrix(const Matrix< T, Prop, Storage, Allocator > &A)
returns true if the matrix is complex
Definition: Functions_BaseInline.cxx:794
Seldon::VectorExpression
Expression between vectors.
Definition: VectorExpression.hxx:8
Seldon::IsSymmetricMatrix
bool IsSymmetricMatrix(const Matrix< T, Prop, Storage, Allocator > &A)
returns true if the matrix is symmetric
Definition: Functions_BaseInline.cxx:778
Seldon::DotProdVector
T1 DotProdVector(const Vector< T1, Storage1, Allocator1 > &X, const Vector< T2, Storage2, Allocator2 > &Y)
Scalar product between two vectors.
Definition: Functions_Vector.cxx:399
Seldon::Matrix
Definition: SeldonHeader.hxx:226
Seldon::CopyVector
void CopyVector(const Vector< T1, Storage1, Allocator1 > &X, Vector< T2, Storage2, Allocator2 > &Y)
Y = X.
Definition: Functions_Vector.cxx:293
Seldon::SolveLuVector
void SolveLuVector(const Matrix< T0, Prop0, Storage0, Allocator0 > &M, Vector< T1, Storage1, Allocator1 > &Y)
Solves a linear system whose matrix has been LU-factorized.
Definition: Functions_MatVect.cxx:1782
Seldon::absSquare
T absSquare(const T &x)
returns the square modulus of z
Definition: CommonInline.cxx:340
Seldon::SolveTriangular
void SolveTriangular(const SeldonUplo &Uplo, const SeldonTranspose &TransA, const SeldonDiag &DiagA, const Matrix< T0, Prop0, RowSparse, Allocator0 > &A, const Vector< T1, VectFull, Allocator1 > &X, Vector< T1, VectFull, Allocator1 > &Y)
solves by triangular upper part or lower part of A
Definition: Functions_MatVect.cxx:1871
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
Seldon namespace.
Definition: Array.cxx:24
Seldon::DotProdConjVector
T1 DotProdConjVector(const Vector< T1, Storage1, Allocator1 > &X, const Vector< T2, Storage2, Allocator2 > &Y)
Scalar product between two vectors conj(X).Y .
Definition: Functions_Vector.cxx:461
Seldon::MltAddMatrix
void MltAddMatrix(const T0 &alpha, const Matrix< T1, Prop1, Storage1, Allocator1 > &A, const Matrix< T2, Prop2, Storage2, Allocator2 > &B, const T3 &beta, Matrix< T4, Prop4, Storage4, Allocator4 > &C)
Multiplies two matrices, and adds the result to a third matrix.
Definition: Functions_Matrix.cxx:617
Seldon::AddMatrix
void AddMatrix(const T0 &alpha, const Matrix< T1, Prop1, Storage1, Allocator1 > &A, Matrix< T2, Prop2, Storage2, Allocator2 > &B)
Adds two matrices.
Definition: Functions_Matrix.cxx:1619
Seldon::ComplexAbs
T ComplexAbs(const T &val)
returns absolute value of val
Definition: CommonInline.cxx:309
Seldon::AddVector
void AddVector(const T0 &alpha, const Vector< T1, Storage1, Allocator1 > &X, Vector< T2, Storage2, Allocator2 > &Y)
Adds two vectors Y = Y + alpha X.
Definition: Functions_Vector.cxx:94