Functions_MatrixArray.cxx
1 // Copyright (C) 2003-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_MATRIX_ARRAY_CXX
22 
23 /*
24  Functions defined in this file:
25  (storage ArrayRowSparse, ArrayColSparse, etc)
26 
27  X = A(i, :)
28  GetRow(A, i, X)
29 
30  X = A(:, j)
31  GetCol(A, j, X)
32 
33  A(i, :) = X
34  SetRow(X, i, A)
35 
36  A(:, j) = X
37  SetCol(X, j, A)
38 
39  alpha.M*X + beta.Y -> Y
40  MltAdd(alpha, M, X, beta, Y)
41 
42  alpha.A + B -> B
43  Add(alpha, A, B)
44 
45  alpha.M -> M
46  Mlt(alpha, M)
47 
48  Highest absolute value of A.
49  MaxAbs(A)
50 
51  1-norm of matrix A.
52  Norm1(A)
53 
54  infinity norm of matrix A.
55  NormInf(A)
56 
57  transpose of matrix A
58  Transpose(A)
59 
60  B = transpose(A)
61  Transpose(A, B)
62 
63  conjugate of transpose of matrix A
64  TransposeConj(A)
65 
66  conjugate of matrix A
67  Conjugate(A)
68 
69  alpha.A*B + beta.C -> C
70  MltAdd(alpha, A, B, beta, C)
71 
72 */
73 
74 namespace Seldon
75 {
76 
77 
79  // GetRow, SetRow //
80 
81 
83 
89  template<class T0, class Allocator0,
90  class T1, class Allocator1>
91  void GetRow(const Matrix<T0, General, ArrayRowSparse, Allocator0>& A,
92  int i, Vector<T1, VectSparse, Allocator1>& X)
93  {
94 #ifdef SELDON_CHECK_BOUNDS
95  int m = A.GetM();
96  if (i < 0 || i >= m)
97  throw WrongIndex("GetRow()",
98  string("Index should be in [0, ") + to_str(m - 1)
99  + "], but is equal to " + to_str(i) + ".");
100 #endif
101 
102  int size_row = A.GetRowSize(i);
103  X.Reallocate(size_row);
104  for (int j = 0; j < size_row; j++)
105  {
106  X.Index(j) = A.Index(i, j);
107  X.Value(j) = A.Value(i, j);
108  }
109  }
110 
111 
113 
119  template <class T0, class Allocator0, class T1, class Allocator1>
120  void GetRow(const Matrix<T0, General, ArrayColSparse, Allocator0>& M,
121  int i, Vector<T1, VectSparse, Allocator1>& X)
122  {
123 #ifdef SELDON_CHECK_BOUNDS
124  int m = M.GetM();
125  if (i < 0 || i >= m)
126  throw WrongIndex("GetRow()",
127  string("Index should be in [0, ") + to_str(m - 1)
128  + "], but is equal to " + to_str(i) + ".");
129 #endif
130 
131  list<pair<int, T0> > vec;
132  for (int j = 0; j < M.GetN(); j++)
133  for (int k = 0; k < M.GetColumnSize(j); k++)
134  if (M.Index(j, k) == i)
135  vec.push_back(make_pair(j, M.Value(j, k)));
136 
137  typename list<pair<int, T0> >::iterator it;
138  X.Reallocate(vec.size());
139  int j = 0;
140  for (it = vec.begin(); it != vec.end(); it++)
141  {
142  X.Index(j) = it->first;
143  X.Value(j) = it->second;
144  j++;
145  }
146  }
147 
148 
150 
156  template <class T0, class Allocator0, class T1, class Allocator1>
157  void GetRow(const Matrix<T0, Symmetric, ArrayRowSymSparse, Allocator0>& M,
158  int i, Vector<T1, VectSparse, Allocator1>& X)
159  {
160 #ifdef SELDON_CHECK_BOUNDS
161  int m = M.GetM();
162  if (i < 0 || i >= m)
163  throw WrongIndex("GetRow()",
164  string("Index should be in [0, ") + to_str(m - 1)
165  + "], but is equal to " + to_str(i) + ".");
166 #endif
167 
168  list<pair<int, T0> > vec;
169  // beginning of row i
170  for (int j = 0; j < i; j++)
171  for (int k = 0; k < M.GetRowSize(j); k++)
172  if (M.Index(j, k) == i)
173  vec.push_back(make_pair(j, M.Value(j, k)));
174 
175  // end of row i
176  for (int k = 0; k < M.GetRowSize(i); k++)
177  vec.push_back(make_pair(M.Index(i, k), M.Value(i, k)));
178 
179  typename list<pair<int, T0> >::iterator it;
180  X.Reallocate(vec.size());
181  int j = 0;
182  for (it = vec.begin(); it != vec.end(); it++)
183  {
184  X.Index(j) = it->first;
185  X.Value(j) = it->second;
186  j++;
187  }
188  }
189 
190 
192 
198  template <class T0, class Allocator0, class T1, class Allocator1>
199  void GetRow(const Matrix<T0, Symmetric, ArrayColSymSparse, Allocator0>& M,
200  int i, Vector<T1, VectSparse, Allocator1>& X)
201  {
202 #ifdef SELDON_CHECK_BOUNDS
203  int m = M.GetM();
204  if (i < 0 || i >= m)
205  throw WrongIndex("GetRow()",
206  string("Index should be in [0, ") + to_str(m - 1)
207  + "], but is equal to " + to_str(i) + ".");
208 #endif
209 
210  list<pair<int, T0> > vec;
211  // beginning of row i
212  for (int k = 0; k < M.GetColumnSize(i); k++)
213  vec.push_back(make_pair(M.Index(i, k), M.Value(i, k)));
214 
215  // end of row i
216  for (int j = i+1; j < M.GetN(); j++)
217  for (int k = 0; k < M.GetColumnSize(j); k++)
218  if (M.Index(j, k) == i)
219  vec.push_back(make_pair(j, M.Value(j, k)));
220 
221  typename list<pair<int, T0> >::iterator it;
222  X.Reallocate(vec.size());
223  int j = 0;
224  for (it = vec.begin(); it != vec.end(); it++)
225  {
226  X.Index(j) = it->first;
227  X.Value(j) = it->second;
228  j++;
229  }
230  }
231 
232 
234 
240  template <class T0, class Allocator0,
241  class T1, class Allocator1>
242  void SetRow(const Vector<T1, VectSparse, Allocator1>& X,
243  int i, Matrix<T0, General, ArrayRowSparse, Allocator0>& M)
244  {
245  M.ClearRow(i);
246  M.ReallocateRow(i, X.GetM());
247  for (int k = 0; k < X.GetM(); k++)
248  {
249  M.Index(i, k) = X.Index(k);
250  M.Value(i, k) = X.Value(k);
251  }
252  }
253 
254 
256 
262  template <class T0, class Allocator0,
263  class T1, class Allocator1>
264  void SetRow(const Vector<T1, VectSparse, Allocator1>& X,
265  int i, Matrix<T0, General, ArrayColSparse, Allocator0>& M)
266  {
267  T0 val;
268  int p = 0;
269  for (int j = 0; j < M.GetN(); j++)
270  {
271  while ( (p < X.GetM()) && (X.Index(p) < j))
272  p++;
273 
274  bool present_X = false;
275  if ( (p < X.GetM()) && (X.Index(p) == j))
276  {
277  present_X = true;
278  val = X.Value(p);
279  }
280 
281  int size_col = M.GetColumnSize(j);
282  bool present_val = false;
283  int k = 0;
284  while ( (k < size_col) && (M.Index(j, k) < i))
285  k++;
286 
287  if ( (k < size_col) && (M.Index(j, k) == i))
288  {
289  if (!present_X)
290  {
291  // reducing size of column
292  if (size_col > 1)
293  {
294  if (k == size_col-1)
295  M.ResizeColumn(j, size_col-1);
296  else
297  {
298  int last_col = M.Index(j, size_col-1);
299  val = M.Value(j, size_col-1);
300  M.ResizeColumn(j, size_col-1);
301  for (int q = k; q < size_col-2; q++)
302  {
303  M.Index(j, q) = M.Index(j, q+1);
304  M.Value(j, q) = M.Value(j, q+1);
305  }
306 
307  M.Index(j, size_col-2) = last_col;
308  M.Value(j, size_col-2) = val;
309  }
310  }
311  else
312  M.ClearColumn(j);
313  }
314  else
315  M.Value(j, k) = val;
316 
317  present_val = true;
318  }
319 
320  if (!present_val && present_X)
321  {
322  // increasing size of column
323  M.ResizeColumn(j, size_col+1);
324  for (int q = size_col; q > k; q--)
325  {
326  M.Index(j, q) = M.Index(j, q-1);
327  M.Value(j, q) = M.Value(j, q-1);
328  }
329 
330  M.Index(j, k) = i;
331  M.Value(j, k) = val;
332  }
333  }
334  }
335 
336 
338 
344  template <class T0, class Allocator0,
345  class T1, class Allocator1>
346  void SetRow(const Vector<T1, VectSparse, Allocator1>& X,
347  int i, Matrix<T0, Symmetric, ArrayRowSymSparse, Allocator0>& M)
348  {
349  T0 val;
350  SetComplexZero(val);
351  int p = 0;
352  for (int j = 0; j < i; j++)
353  {
354  while ( (p < X.GetM()) && (X.Index(p) < j))
355  p++;
356 
357  bool present_X = false;
358  if ( (p < X.GetM()) && (X.Index(p) == j))
359  {
360  present_X = true;
361  val = X.Value(p);
362  p++;
363  }
364 
365  int size_row = M.GetRowSize(j);
366  bool present_val = false;
367  int k = 0;
368  while ( (k < size_row) && (M.Index(j, k) < i))
369  k++;
370 
371  if ( (k < size_row) && (M.Index(j, k) == i))
372  {
373  if (!present_X)
374  {
375  // reducing size of row
376  if (size_row > 1)
377  {
378  if (k == size_row-1)
379  M.ResizeRow(j, size_row-1);
380  else
381  {
382  int last_col = M.Index(j, size_row-1);
383  val = M.Value(j, size_row-1);
384  M.ResizeRow(j, size_row-1);
385  for (int q = k; q < size_row-2; q++)
386  {
387  M.Index(j, q) = M.Index(j, q+1);
388  M.Value(j, q) = M.Value(j, q+1);
389  }
390 
391  M.Index(j, size_row-2) = last_col;
392  M.Value(j, size_row-2) = val;
393  }
394  }
395  else
396  M.ClearRow(j);
397  }
398  else
399  M.Value(j, k) = val;
400 
401  present_val = true;
402  }
403 
404  if (!present_val && present_X)
405  {
406  // increasing size of row
407  M.ResizeRow(j, size_row+1);
408  for (int q = size_row; q > k; q--)
409  {
410  M.Index(j, q) = M.Index(j, q-1);
411  M.Value(j, q) = M.Value(j, q-1);
412  }
413 
414  M.Index(j, k) = i;
415  M.Value(j, k) = val;
416  }
417  }
418 
419  // and changing row i
420  M.ClearRow(i);
421  if (p < X.GetM())
422  {
423  M.ReallocateRow(i, X.GetM() - p);
424  int k = 0;
425  while (p < X.GetM())
426  {
427  M.Index(i, k) = X.Index(p);
428  M.Value(i, k) = X.Value(p);
429  k++;
430  p++;
431  }
432  }
433  }
434 
435 
437 
443  template <class T0, class Allocator0,
444  class T1, class Allocator1>
445  void SetRow(const Vector<T1, VectSparse, Allocator1>& X,
446  int i, Matrix<T0, Symmetric, ArrayColSymSparse, Allocator0>& M)
447  {
448  T1 val;
449  int p = 0;
450  // setting column i
451  M.ClearColumn(i);
452  while ( (p < X.GetM()) && (X.Index(p) <= i))
453  p++;
454 
455  if (p > 0)
456  {
457  M.ReallocateColumn(i, p);
458  for (int k = 0; k < p; k++)
459  {
460  M.Index(i, k) = X.Index(k);
461  M.Value(i, k) = X.Value(k);
462  }
463  }
464 
465  // then modifying last columns
466  for (int j = i+1; j < M.GetN(); j++)
467  {
468  while ( (p < X.GetM()) && (X.Index(p) < j))
469  p++;
470 
471  bool present_X = false;
472  if ( (p < X.GetM()) && (X.Index(p) == j))
473  {
474  present_X = true;
475  val = X.Value(p);
476  }
477 
478  int size_col = M.GetColumnSize(j);
479  bool present_val = false;
480  int k = 0;
481  while ( (k < size_col) && (M.Index(j, k) < i))
482  k++;
483 
484  if ( (k < size_col) && (M.Index(j, k) == i))
485  {
486  if (!present_X)
487  {
488  // reducing size of column
489  if (size_col > 1)
490  {
491  if (k == size_col-1)
492  M.ResizeColumn(j, size_col-1);
493  else
494  {
495  int last_col = M.Index(j, size_col-1);
496  val = M.Value(j, size_col-1);
497  M.ResizeColumn(j, size_col-1);
498  for (int q = k; q < size_col-2; q++)
499  {
500  M.Index(j, q) = M.Index(j, q+1);
501  M.Value(j, q) = M.Value(j, q+1);
502  }
503 
504  M.Index(j, size_col-2) = last_col;
505  M.Value(j, size_col-2) = val;
506  }
507  }
508  else
509  M.ClearColumn(j);
510  }
511  else
512  M.Value(j, k) = val;
513 
514  present_val = true;
515  }
516 
517  if (!present_val && present_X)
518  {
519  // increasing size of column
520  M.ResizeColumn(j, size_col+1);
521  for (int q = size_col; q > k; q--)
522  {
523  M.Index(j, q) = M.Index(j, q-1);
524  M.Value(j, q) = M.Value(j, q-1);
525  }
526 
527  M.Index(j, k) = i;
528  M.Value(j, k) = val;
529  }
530 
531  }
532  }
533 
534 
535  // GetRow, SetRow //
537 
538 
540  // GetCol, SetCol //
541 
542 
544 
550  template <class T0, class Allocator0, class T1, class Allocator1>
551  void GetCol(const Matrix<T0, General, ArrayRowSparse, Allocator0>& M,
552  int j, Vector<T1, VectSparse, Allocator1>& X)
553  {
554 #ifdef SELDON_CHECK_BOUNDS
555  int n = M.GetN();
556  if (j < 0 || j >= n)
557  throw WrongIndex("GetCol()",
558  string("Index should be in [0, ") + to_str(n - 1)
559  + "], but is equal to " + to_str(j) + ".");
560 #endif
561 
562  int m = M.GetM();
563 
564  list<pair<int, T0> > vec;
565  for (int i = 0; i < m; i++)
566  for (int k = 0; k < M.GetRowSize(i); k++)
567  if (M.Index(i, k) == j)
568  vec.push_back(make_pair(i, M.Value(i, k)));
569 
570  typename list<pair<int, T0> >::iterator it;
571  X.Reallocate(vec.size());
572  int i = 0;
573  for (it = vec.begin(); it != vec.end(); it++)
574  {
575  X.Index(i) = it->first;
576  X.Value(i) = it->second;
577  i++;
578  }
579  }
580 
581 
583 
589  template<class T0, class Allocator0,
590  class T1, class Allocator1>
591  void GetCol(const Matrix<T0, General, ArrayColSparse, Allocator0>& A,
592  int j, Vector<T1, VectSparse, Allocator1>& X)
593  {
594 #ifdef SELDON_CHECK_BOUNDS
595  int n = A.GetN();
596  if (j < 0 || j >= n)
597  throw WrongIndex("GetCol()",
598  string("Index should be in [0, ") + to_str(n - 1)
599  + "], but is equal to " + to_str(j) + ".");
600 #endif
601 
602  int size_col = A.GetColumnSize(j);
603  X.Reallocate(size_col);
604  for (int k = 0; k < size_col; k++)
605  {
606  X.Index(k) = A.Index(j, k);
607  X.Value(k) = A.Value(j, k);
608  }
609  }
610 
611 
613 
619  template <class T0, class Allocator0, class T1, class Allocator1>
620  void GetCol(const Matrix<T0, Symmetric, ArrayRowSymSparse, Allocator0>& M,
621  int j, Vector<T1, VectSparse, Allocator1>& X)
622  {
623  // symmetric matrix row = col
624  GetRow(M, j, X);
625  }
626 
627 
629 
635  template <class T0, class Allocator0, class T1, class Allocator1>
636  void GetCol(const Matrix<T0, Symmetric, ArrayColSymSparse, Allocator0>& M,
637  int j, Vector<T1, VectSparse, Allocator1>& X)
638  {
639  // symmetric matrix row = col
640  GetRow(M, j, X);
641  }
642 
643 
645 
651  template <class T0, class Allocator0,
652  class T1, class Allocator1>
653  void SetCol(const Vector<T1, VectSparse, Allocator1>& X,
654  int j, Matrix<T0, General, ArrayRowSparse, Allocator0>& M)
655  {
656  T0 val;
657  SetComplexZero(val);
658  int p = 0;
659  for (int i = 0; i < M.GetM(); i++)
660  {
661  while ( (p < X.GetM()) && (X.Index(p) < i))
662  p++;
663 
664  bool present_X = false;
665  if ( (p < X.GetM()) && (X.Index(p) == i))
666  {
667  present_X = true;
668  val = X.Value(p);
669  }
670 
671  int size_row = M.GetRowSize(i);
672  bool present_val = false;
673  int k = 0;
674  while ( (k < size_row) && (M.Index(i, k) < j))
675  k++;
676 
677 
678  if ( (k < size_row) && (M.Index(i, k) == j))
679  {
680  if (!present_X)
681  {
682  // reducing size of row
683  if (size_row > 1)
684  {
685  if (k == size_row-1)
686  M.ResizeRow(j, size_row-1);
687  else
688  {
689  int last_row = M.Index(i, size_row-1);
690  val = M.Value(i, size_row-1);
691  M.ResizeRow(i, size_row-1);
692  for (int q = k; q < size_row-2; q++)
693  {
694  M.Index(i, q) = M.Index(i, q+1);
695  M.Value(i, q) = M.Value(i, q+1);
696  }
697 
698  M.Index(i, size_row-2) = last_row;
699  M.Value(i, size_row-2) = val;
700  }
701  }
702  else
703  M.ClearRow(i);
704  }
705  else
706  M.Value(i, k) = val;
707 
708  present_val = true;
709  }
710 
711  if (!present_val && present_X)
712  {
713  // increasing size of row
714  M.ResizeRow(i, size_row+1);
715  for (int q = size_row; q > k; q--)
716  {
717  M.Index(i, q) = M.Index(i, q-1);
718  M.Value(i, q) = M.Value(i, q-1);
719  }
720 
721  M.Index(i, k) = j;
722  M.Value(i, k) = val;
723  }
724  }
725  }
726 
727 
729 
735  template <class T0, class Allocator0,
736  class T1, class Allocator1>
737  void SetCol(const Vector<T1, VectSparse, Allocator1>& X,
738  int j, Matrix<T0, General, ArrayColSparse, Allocator0>& M)
739  {
740  M.ClearColumn(j);
741  M.ReallocateColumn(j, X.GetM());
742  for (int k = 0; k < X.GetM(); k++)
743  {
744  M.Index(j, k) = X.Index(k);
745  M.Value(j, k) = X.Value(k);
746  }
747  }
748 
749 
751 
757  template <class T0, class Allocator0,
758  class T1, class Allocator1>
759  void SetCol(const Vector<T1, VectSparse, Allocator1>& X,
760  int j, Matrix<T0, Symmetric, ArrayRowSymSparse, Allocator0>& M)
761  {
762  // symmetric matrix, row = column
763  SetRow(X, j, M);
764  }
765 
766 
768 
774  template <class T0, class Allocator0,
775  class T1, class Allocator1>
776  void SetCol(const Vector<T1, VectSparse, Allocator1>& X,
777  int j, Matrix<T0, Symmetric, ArrayColSymSparse, Allocator0>& M)
778  {
779  // symmetric matrix, row = column
780  SetRow(X, j, M);
781  }
782 
783 
784  // GetCol, SetCol //
786 
787 
789  // Mlt //
790 
791  /*** ArrayRowSymSparse ***/
792 
793 
794  template<class T1, class T2, class T3,
795  class Allocator1, class Allocator2, class Allocator3>
796  void MltVector(const Matrix<T1, Symmetric,
797  ArrayRowSymSparse, Allocator1>& A,
798  const Vector<T2, VectFull, Allocator2>& B,
799  Vector<T3, VectFull, Allocator3>& C)
800  {
801  T1 val;
802  int m = A.GetM(), n, p;
803  C.Fill(0);
804  for (int i = 0 ; i < m ; i++)
805  {
806  n = A.GetRowSize(i);
807  for (int k = 0; k < n ; k++)
808  {
809  p = A.Index(i, k);
810  val = A.Value(i, k);
811 
812  if (p == i)
813  C(i) += val * B(i);
814  else
815  {
816  C(i) += val * B(p);
817  C(p) += val * B(i);
818  }
819  }
820  }
821  }
822 
823 
824  template<class T1, class T2, class T3,
825  class Allocator1, class Allocator2, class Allocator3>
826  void MltVector(const SeldonTranspose& Trans, const Matrix<T1, Symmetric,
827  ArrayRowSymSparse, Allocator1>& A,
828  const Vector<T2, VectFull, Allocator2>& B,
829  Vector<T3, VectFull, Allocator3>& C)
830  {
831  if (!Trans.ConjTrans())
832  {
833  MltVector(A, B, C);
834  return;
835  }
836 
837  C.Fill(0);
838 
839  int m = A.GetM(), n, p;
840  T1 val;
841  for (int i = 0 ; i < m ; i++)
842  {
843  n = A.GetRowSize(i);
844  for (int k = 0; k < n ; k++)
845  {
846  p = A.Index(i, k);
847  val = conjugate(A.Value(i, k));
848 
849  if (p == i)
850  C(i) += val * B(i);
851  else
852  {
853  C(i) += val * B(p);
854  C(p) += val * B(i);
855  }
856  }
857  }
858  }
859 
860 
861  /*** ArrayColSymSparse ***/
862 
863 
864  template<class T1, class T2, class T3,
865  class Allocator1, class Allocator2, class Allocator3>
866  void MltVector(const Matrix<T1, Symmetric,
867  ArrayColSymSparse, Allocator1>& A,
868  const Vector<T2, VectFull, Allocator2>& B,
869  Vector<T3, VectFull, Allocator3>& C)
870  {
871  C.Fill(0);
872  int m = A.GetM(), n, p;
873  T1 val;
874  for (int i = 0 ; i < m ; i++)
875  {
876  n = A.GetColumnSize(i);
877  for (int k = 0; k < n ; k++)
878  {
879  p = A.Index(i, k);
880  val = A.Value(i, k);
881 
882  if (p == i)
883  C(i) += val * B(i);
884  else
885  {
886  C(i) += val * B(p);
887  C(p) += val * B(i);
888  }
889  }
890  }
891  }
892 
893 
894  template<class T1, class T2, class T3,
895  class Allocator1, class Allocator2, class Allocator3>
896  void MltVector(const SeldonTranspose& Trans, const Matrix<T1, Symmetric,
897  ArrayColSymSparse, Allocator1>& A,
898  const Vector<T2, VectFull, Allocator2>& B,
899  Vector<T3, VectFull, Allocator3>& C)
900  {
901  if (!Trans.ConjTrans())
902  {
903  MltVector(A, B, C);
904  return;
905  }
906 
907  C.Fill(0);
908 
909  int m = A.GetM(), n, p;
910  T3 val;
911  for (int i = 0 ; i < m ; i++)
912  {
913  n = A.GetColumnSize(i);
914  for (int k = 0; k < n ; k++)
915  {
916  p = A.Index(i, k);
917  val = conjugate(A.Value(i, k));
918 
919  if (p == i)
920  C(i) += val * B(i);
921  else
922  {
923  C(i) += val * B(p);
924  C(p) += val * B(i);
925  }
926  }
927  }
928  }
929 
930 
931  /*** ArrayRowSparse ***/
932 
933 
934  template<class T1, class T2, class T3,
935  class Allocator1, class Allocator2, class Allocator3>
936  void MltVector(const Matrix<T1, General, ArrayRowSparse, Allocator1>& A,
937  const Vector<T2, VectFull, Allocator2>& B,
938  Vector<T3, VectFull, Allocator3>& C)
939  {
940  T3 zero, temp; SetComplexZero(zero);
941  int m = A.GetM(), n;
942  for (int i = 0 ; i < m ; i++)
943  {
944  n = A.GetRowSize(i);
945  temp = zero;
946  for (int k = 0; k < n ; k++)
947  temp += A.Value(i, k)*B(A.Index(i, k));
948 
949  C(i) = temp;
950  }
951  }
952 
953 
954  template<class T1, class T2, class T3,
955  class Allocator1, class Allocator2, class Allocator3>
956  void MltVector(const SeldonTranspose& Trans,
957  const Matrix<T1, General, ArrayRowSparse, Allocator1>& A,
958  const Vector<T2, VectFull, Allocator2>& B,
959  Vector<T3, VectFull, Allocator3>& C)
960  {
961  if (Trans.NoTrans())
962  {
963  MltVector(A, B, C);
964  return;
965  }
966 
967  C.Fill(0);
968  int m = A.GetM(), n;
969 
970  if (Trans.Trans())
971  {
972  for (int i = 0 ; i < m ; i++)
973  {
974  n = A.GetRowSize(i);
975  for (int k = 0; k < n ; k++)
976  C(A.Index(i, k)) += A.Value(i, k)*B(i);
977  }
978  }
979  else
980  {
981  for (int i = 0 ; i < m ; i++)
982  {
983  n = A.GetRowSize(i);
984  for (int k = 0; k < n ; k++)
985  C(A.Index(i, k)) += conjugate(A.Value(i, k))*B(i);
986  }
987  }
988  }
989 
990 
991  /*** ArrayColSparse ***/
992 
993 
994  template<class T1, class T2, class T3,
995  class Allocator1, class Allocator2, class Allocator3>
996  void MltVector(const Matrix<T1, General, ArrayColSparse, Allocator1>& A,
997  const Vector<T2, VectFull, Allocator2>& B,
998  Vector<T3, VectFull, Allocator3>& C)
999  {
1000  C.Fill(0);
1001  for (int i = 0 ; i < A.GetN(); i++)
1002  for (int k = 0; k < A.GetColumnSize(i); k++)
1003  C(A.Index(i, k)) += A.Value(i, k) * B(i);
1004  }
1005 
1006 
1007  template<class T1, class T2, class T3,
1008  class Allocator1, class Allocator2, class Allocator3>
1009  void MltVector(const SeldonTranspose& Trans,
1010  const Matrix<T1, General, ArrayColSparse, Allocator1>& A,
1011  const Vector<T2, VectFull, Allocator2>& B,
1012  Vector<T3, VectFull, Allocator3>& C)
1013  {
1014  if (Trans.NoTrans())
1015  {
1016  MltVector(A, B, C);
1017  return;
1018  }
1019 
1020  T3 zero, temp;
1021  SetComplexZero(zero);
1022 
1023  if (Trans.Trans())
1024  {
1025  for (int i = 0 ; i < A.GetN(); i++)
1026  {
1027  temp = zero;
1028  for (int k = 0; k < A.GetColumnSize(i); k++)
1029  temp += A.Value(i, k) * B(A.Index(i, k));
1030 
1031  C(i) = temp;
1032  }
1033  }
1034  else
1035  {
1036  for (int i = 0 ; i < A.GetN(); i++)
1037  {
1038  temp = zero;
1039  for (int k = 0; k < A.GetColumnSize(i); k++)
1040  temp += conjugate(A.Value(i, k)) * B(A.Index(i, k));
1041 
1042  C(i) = temp;
1043  }
1044  }
1045  }
1046 
1047 
1048  // Mlt //
1050 
1051 
1053  // MltAdd //
1054 
1055 
1056  /*** ArrayRowSymSparse ***/
1057 
1058 
1059  template<class T0, class T1, class T2, class T3, class T4,
1060  class Allocator1, class Allocator2, class Allocator3>
1061  void MltAddVector(const T0& alpha, const Matrix<T1, Symmetric,
1062  ArrayRowSymSparse, Allocator1>& A,
1063  const Vector<T2, VectFull, Allocator2>& B,
1064  const T4& beta,
1065  Vector<T3, VectFull, Allocator3>& C)
1066  {
1067  T4 zero; T0 one;
1068  SetComplexZero(zero);
1069  SetComplexOne(one);
1070 
1071  if (beta == zero)
1072  C.Fill(0);
1073  else
1074  Mlt(beta, C);
1075 
1076  int m = A.GetM(), n, p;
1077  T1 val;
1078  if (alpha == one)
1079  {
1080  for (int i = 0 ; i < m ; i++)
1081  {
1082  n = A.GetRowSize(i);
1083  for (int k = 0; k < n ; k++)
1084  {
1085  p = A.Index(i, k);
1086  val = A.Value(i, k);
1087 
1088  if (p == i)
1089  C(i) += val * B(i);
1090  else
1091  {
1092  C(i) += val * B(p);
1093  C(p) += val * B(i);
1094  }
1095  }
1096  }
1097  }
1098  else
1099  {
1100  // alpha different from 1
1101  for (int i = 0 ; i < m ; i++)
1102  {
1103  n = A.GetRowSize(i);
1104  for (int k = 0; k < n ; k++)
1105  {
1106  p = A.Index(i, k);
1107  val = A.Value(i, k);
1108 
1109  if (p==i)
1110  C(i) += alpha * val * B(i);
1111  else
1112  {
1113  C(i) += alpha * val * B(p);
1114  C(p) += alpha * val * B(i);
1115  }
1116  }
1117  }
1118  }
1119  }
1120 
1121 
1122  template<class T0, class T1, class T2, class T3, class T4,
1123  class Allocator1, class Allocator2, class Allocator3>
1124  void MltAddVector(const T0& alpha,
1125  const SeldonTranspose& Trans, const Matrix<T1, Symmetric,
1126  ArrayRowSymSparse, Allocator1>& A,
1127  const Vector<T2, VectFull, Allocator2>& B,
1128  const T4& beta,
1129  Vector<T3, VectFull, Allocator3>& C)
1130  {
1131  if (!Trans.ConjTrans())
1132  {
1133  MltAddVector(alpha, A, B, beta, C);
1134  return;
1135  }
1136 
1137  T4 zero; T0 one;
1138  SetComplexZero(zero);
1139  SetComplexOne(one);
1140 
1141  if (beta == zero)
1142  C.Fill(0);
1143  else
1144  Mlt(beta, C);
1145 
1146  int m = A.GetM(), n, p;
1147  T1 val;
1148  if (alpha == one)
1149  {
1150  for (int i = 0 ; i < m ; i++)
1151  {
1152  n = A.GetRowSize(i);
1153  for (int k = 0; k < n ; k++)
1154  {
1155  p = A.Index(i, k);
1156  val = conjugate(A.Value(i, k));
1157 
1158  if (p == i)
1159  C(i) += val * B(i);
1160  else
1161  {
1162  C(i) += val * B(p);
1163  C(p) += val * B(i);
1164  }
1165  }
1166  }
1167  }
1168  else
1169  {
1170  // alpha different from 1
1171  for (int i = 0 ; i < m ; i++)
1172  {
1173  n = A.GetRowSize(i);
1174  for (int k = 0; k < n ; k++)
1175  {
1176  p = A.Index(i, k);
1177  val = conjugate(A.Value(i, k));
1178 
1179  if (p==i)
1180  C(i) += alpha * val * B(i);
1181  else
1182  {
1183  C(i) += alpha * val * B(p);
1184  C(p) += alpha * val * B(i);
1185  }
1186  }
1187  }
1188  }
1189  }
1190 
1191 
1192  /*** ArrayColSymSparse ***/
1193 
1194 
1195  template<class T0, class T1, class T2, class T3, class T4,
1196  class Allocator1, class Allocator2, class Allocator3>
1197  void MltAddVector(const T0& alpha, const Matrix<T1, Symmetric,
1198  ArrayColSymSparse, Allocator1>& A,
1199  const Vector<T2, VectFull, Allocator2>& B,
1200  const T4& beta,
1201  Vector<T3, VectFull, Allocator3>& C)
1202  {
1203  T4 zero; T0 one;
1204  SetComplexZero(zero);
1205  SetComplexOne(one);
1206 
1207  if (beta == zero)
1208  C.Fill(0);
1209  else
1210  Mlt(beta, C);
1211 
1212  int m = A.GetM(), n, p;
1213  T1 val;
1214  if (alpha == one)
1215  {
1216  for (int i = 0 ; i < m ; i++)
1217  {
1218  n = A.GetColumnSize(i);
1219  for (int k = 0; k < n ; k++)
1220  {
1221  p = A.Index(i, k);
1222  val = A.Value(i, k);
1223 
1224  if (p == i)
1225  C(i) += val * B(i);
1226  else
1227  {
1228  C(i) += val * B(p);
1229  C(p) += val * B(i);
1230  }
1231  }
1232  }
1233  }
1234  else
1235  {
1236  // alpha different from 1
1237  for (int i = 0 ; i < m ; i++)
1238  {
1239  n = A.GetColumnSize(i);
1240  for (int k = 0; k < n ; k++)
1241  {
1242  p = A.Index(i, k);
1243  val = A.Value(i, k);
1244 
1245  if (p==i)
1246  C(i) += alpha * val * B(i);
1247  else
1248  {
1249  C(i) += alpha * val * B(p);
1250  C(p) += alpha * val * B(i);
1251  }
1252  }
1253  }
1254  }
1255  }
1256 
1257 
1258  template<class T0, class T1, class T2, class T3, class T4,
1259  class Allocator1, class Allocator2, class Allocator3>
1260  void MltAddVector(const T0& alpha,
1261  const SeldonTranspose& Trans, const Matrix<T1, Symmetric,
1262  ArrayColSymSparse, Allocator1>& A,
1263  const Vector<T2, VectFull, Allocator2>& B,
1264  const T4& beta,
1265  Vector<T3, VectFull, Allocator3>& C)
1266  {
1267  if (!Trans.ConjTrans())
1268  {
1269  MltAddVector(alpha, A, B, beta, C);
1270  return;
1271  }
1272 
1273  T4 zero; T0 one;
1274  SetComplexZero(zero);
1275  SetComplexOne(one);
1276 
1277  if (beta == zero)
1278  C.Fill(0);
1279  else
1280  Mlt(beta, C);
1281 
1282  int m = A.GetM(), n, p;
1283  T3 val;
1284  if (alpha == one)
1285  {
1286  for (int i = 0 ; i < m ; i++)
1287  {
1288  n = A.GetColumnSize(i);
1289  for (int k = 0; k < n ; k++)
1290  {
1291  p = A.Index(i, k);
1292  val = conjugate(A.Value(i, k));
1293 
1294  if (p == i)
1295  C(i) += val * B(i);
1296  else
1297  {
1298  C(i) += val * B(p);
1299  C(p) += val * B(i);
1300  }
1301  }
1302  }
1303  }
1304  else
1305  {
1306  // alpha different from 1
1307  for (int i = 0 ; i < m ; i++)
1308  {
1309  n = A.GetColumnSize(i);
1310  for (int k = 0; k < n ; k++)
1311  {
1312  p = A.Index(i, k);
1313  val = alpha * conjugate(A.Value(i, k));
1314 
1315  if (p==i)
1316  C(i) += val * B(i);
1317  else
1318  {
1319  C(i) += val * B(p);
1320  C(p) += val * B(i);
1321  }
1322  }
1323  }
1324  }
1325  }
1326 
1327 
1328  /*** ArrayRowSparse ***/
1329 
1330 
1331  template<class T0, class T1, class T2, class T3, class T4,
1332  class Allocator1, class Allocator2, class Allocator3>
1333  void MltAddVector(const T0& alpha,
1334  const Matrix<T1, General, ArrayRowSparse, Allocator1>& A,
1335  const Vector<T2, VectFull, Allocator2>& B,
1336  const T4& beta,
1337  Vector<T3, VectFull, Allocator3>& C)
1338  {
1339  T4 zero; SetComplexZero(zero);
1340  T0 one; SetComplexOne(one);
1341 
1342  if (beta == zero)
1343  C.Fill(0);
1344  else
1345  Mlt(beta, C);
1346 
1347  int m = A.GetM(), n, p;
1348  T1 val;
1349  if (alpha == one)
1350  {
1351  for (int i = 0 ; i < m ; i++)
1352  {
1353  n = A.GetRowSize(i);
1354  for (int k = 0; k < n ; k++)
1355  {
1356  p = A.Index(i, k);
1357  val = A.Value(i, k);
1358  C(i) += val * B(p);
1359  }
1360  }
1361  }
1362  else // alpha != 1.
1363  {
1364  for (int i = 0 ; i < m ; i++)
1365  {
1366  n = A.GetRowSize(i);
1367  for (int k = 0; k < n ; k++)
1368  {
1369  p = A.Index(i, k);
1370  val = A.Value(i, k);
1371  C(i) += alpha * val * B(p);
1372  }
1373  }
1374  }
1375  }
1376 
1377 
1378  template<class T0, class T1, class T2, class T3, class T4,
1379  class Allocator1, class Allocator2, class Allocator3>
1380  void MltAddVector(const T0& alpha,
1381  const SeldonTranspose& Trans,
1382  const Matrix<T1, General, ArrayRowSparse, Allocator1>& A,
1383  const Vector<T2, VectFull, Allocator2>& B,
1384  const T4& beta,
1385  Vector<T3, VectFull, Allocator3>& C)
1386  {
1387  if (Trans.NoTrans())
1388  {
1389  MltAddVector(alpha, A, B, beta, C);
1390  return;
1391  }
1392 
1393  T4 zero; SetComplexZero(zero);
1394  T0 one; SetComplexOne(one);
1395 
1396  if (beta == zero)
1397  C.Fill(0);
1398  else
1399  Mlt(beta, C);
1400 
1401  int m = A.GetM(), n, p;
1402  T1 val;
1403 
1404  if (Trans.Trans())
1405  {
1406  if (alpha == one)
1407  {
1408  for (int i = 0 ; i < m ; i++)
1409  {
1410  n = A.GetRowSize(i);
1411  for (int k = 0; k < n ; k++)
1412  {
1413  p = A.Index(i, k);
1414  val = A.Value(i, k);
1415  C(p) += val * B(i);
1416  }
1417  }
1418  }
1419  else // alpha != 1.
1420  {
1421  for (int i = 0 ; i < m ; i++)
1422  {
1423  n = A.GetRowSize(i);
1424  for (int k = 0; k < n ; k++)
1425  {
1426  p = A.Index(i, k);
1427  val = A.Value(i, k);
1428  C(p) += alpha * val * B(i);
1429  }
1430  }
1431  }
1432  }
1433  else
1434  {
1435  if (alpha == one)
1436  {
1437  for (int i = 0 ; i < m ; i++)
1438  {
1439  n = A.GetRowSize(i);
1440  for (int k = 0; k < n ; k++)
1441  {
1442  p = A.Index(i, k);
1443  val = conjugate(A.Value(i, k));
1444  C(p) += val * B(i);
1445  }
1446  }
1447  }
1448  else // alpha != 1.
1449  {
1450  for (int i = 0 ; i < m ; i++)
1451  {
1452  n = A.GetRowSize(i);
1453  for (int k = 0; k < n ; k++)
1454  {
1455  p = A.Index(i, k);
1456  val = conjugate(A.Value(i, k));
1457  C(p) += alpha * val * B(i);
1458  }
1459  }
1460  }
1461  }
1462  }
1463 
1464 
1465  /*** ArrayColSparse ***/
1466 
1467 
1468  template<class T0, class T1, class T2, class T3, class T4,
1469  class Allocator1, class Allocator2, class Allocator3>
1470  void MltAddVector(const T0& alpha,
1471  const Matrix<T1, General, ArrayColSparse, Allocator1>& A,
1472  const Vector<T2, VectFull, Allocator2>& B,
1473  const T4& beta,
1474  Vector<T3, VectFull, Allocator3>& C)
1475  {
1476  T4 zero; T0 one;
1477  SetComplexZero(zero);
1478  SetComplexOne(one);
1479 
1480  if (beta == zero)
1481  C.Fill(0);
1482  else
1483  Mlt(beta, C);
1484 
1485  if (alpha == one)
1486  {
1487  for (int i = 0 ; i < A.GetN(); i++)
1488  for (int k = 0; k < A.GetColumnSize(i); k++)
1489  C(A.Index(i, k)) += A.Value(i, k) * B(i);
1490  }
1491  else // alpha != 1.
1492  {
1493  for (int i = 0 ; i < A.GetN(); i++)
1494  for (int k = 0; k < A.GetColumnSize(i); k++)
1495  C(A.Index(i, k)) += alpha * A.Value(i, k) * B(i);
1496  }
1497  }
1498 
1499 
1500  template<class T0, class T1, class T2, class T3, class T4,
1501  class Allocator1, class Allocator2, class Allocator3>
1502  void MltAddVector(const T0& alpha,
1503  const SeldonTranspose& Trans,
1504  const Matrix<T1, General, ArrayColSparse, Allocator1>& A,
1505  const Vector<T2, VectFull, Allocator2>& B,
1506  const T4& beta,
1507  Vector<T3, VectFull, Allocator3>& C)
1508  {
1509  if (Trans.NoTrans())
1510  {
1511  MltAddVector(alpha, A, B, beta, C);
1512  return;
1513  }
1514 
1515  T4 zero; T0 one;
1516  SetComplexZero(zero);
1517  SetComplexOne(one);
1518 
1519  if (beta == zero)
1520  C.Fill(0);
1521  else
1522  Mlt(beta, C);
1523 
1524  if (Trans.Trans())
1525  {
1526  if (alpha == one)
1527  {
1528  for (int i = 0 ; i < A.GetN(); i++)
1529  for (int k = 0; k < A.GetColumnSize(i); k++)
1530  C(i) += A.Value(i, k) * B(A.Index(i, k));
1531  }
1532  else // alpha != 1.
1533  {
1534  for (int i = 0 ; i < A.GetN(); i++)
1535  for (int k = 0; k < A.GetColumnSize(i); k++)
1536  C(i) += alpha * A.Value(i, k) * B(A.Index(i, k));
1537  }
1538  }
1539  else
1540  {
1541  if (alpha == one)
1542  {
1543  for (int i = 0 ; i < A.GetN(); i++)
1544  for (int k = 0; k < A.GetColumnSize(i); k++)
1545  C(i) += conjugate(A.Value(i, k)) * B(A.Index(i, k));
1546  }
1547  else // alpha != 1.
1548  {
1549  for (int i = 0 ; i < A.GetN(); i++)
1550  for (int k = 0; k < A.GetColumnSize(i); k++)
1551  C(i) += alpha * conjugate(A.Value(i, k)) * B(A.Index(i, k));
1552  }
1553  }
1554  }
1555 
1556 
1557  // MltAdd //
1559 
1560 
1561 
1563  // Add //
1564 
1565 
1566  template<class T0, class T1, class T2, class Allocator1, class Allocator2>
1567  void AddMatrix(const T0& alpha,
1568  const Matrix<T1, General, ArrayRowSparse, Allocator1>& A,
1569  Matrix<T2, General, ArrayRowSparse, Allocator2>& B)
1570  {
1571  int m = B.GetM(), n;
1572  Vector<T2, VectFull, Allocator2> value(B.GetN());
1573  for (int i = 0 ; i < m ; i++)
1574  {
1575  n = A.GetRowSize(i);
1576  for (int j = 0; j < n; j++)
1577  value(j) = alpha*A.Value(i, j);
1578 
1579  B.AddInteractionRow(i, n, A.GetIndex(i), value.GetData());
1580  }
1581  }
1582 
1583 
1584  template<class T0, class T1, class T2, class Allocator1, class Allocator2>
1585  void AddMatrix(const T0& alpha,
1586  const Matrix<T1, General, ArrayColSparse, Allocator1>& A,
1587  Matrix<T2, General, ArrayColSparse, Allocator2>& B)
1588  {
1589  int m = B.GetN(), n;
1590  Vector<T2, VectFull, Allocator2> value;
1591  for (int i = 0 ; i < m ; i++)
1592  {
1593  n = A.GetColumnSize(i);
1594  value.Reallocate(n);
1595  for (int j = 0; j < n; j++)
1596  value(j) = A.Value(i, j);
1597 
1598  Mlt(alpha, value);
1599  B.AddInteractionColumn(i, n, A.GetIndex(i), value.GetData());
1600  }
1601  }
1602 
1603 
1604  template<class T0, class T1, class T2, class Allocator1, class Allocator2>
1605  void AddMatrix(const T0& alpha,
1606  const Matrix<T1, Symmetric, ArrayRowSymSparse, Allocator1>& A,
1607  Matrix<T2, Symmetric, ArrayRowSymSparse, Allocator2>& B)
1608  {
1609  int m = B.GetM(),n;
1610  Vector<T2, VectFull, Allocator2> value(B.GetN());
1611  for (int i = 0 ; i < m ; i++)
1612  {
1613  n = A.GetRowSize(i);
1614  for (int j = 0; j < n; j++)
1615  value(j) = alpha*A.Value(i, j);
1616 
1617  B.AddInteractionRow(i, n, A.GetIndex(i), value.GetData());
1618  }
1619  }
1620 
1621 
1622  template<class T0, class T1, class T2, class Allocator1, class Allocator2>
1623  void AddMatrix(const T0& alpha,
1624  const Matrix<T1, Symmetric, ArrayColSymSparse, Allocator1>& A,
1625  Matrix<T2, Symmetric, ArrayColSymSparse, Allocator2>& B)
1626  {
1627  int m = B.GetN(), n;
1628  Vector<T2, VectFull, Allocator2> value;
1629  for (int i = 0 ; i < m ; i++)
1630  {
1631  n = A.GetColumnSize(i);
1632  value.Reallocate(n);
1633  for (int j = 0; j < n; j++)
1634  value(j) = A.Value(i, j);
1635 
1636  Mlt(alpha, value);
1637  B.AddInteractionColumn(i, n, A.GetIndex(i), value.GetData());
1638  }
1639  }
1640 
1641 
1642  // C = C + alpha complex(A,B)
1643  template<class T0, class T1, class T2, class T3, class Allocator1,
1644  class Allocator2, class Allocator3>
1645  void AddMatrix(const T0& alpha,
1646  const Matrix<T1, General, ArrayRowSparse, Allocator1>& A,
1647  const Matrix<T2, General, ArrayRowSparse, Allocator2>& B,
1648  Matrix<complex<T3>, General, ArrayRowSparse, Allocator3>& C)
1649  {
1650  int m = B.GetM(),n1,n2,size_row;;
1651  Vector<complex<T3>, VectFull, Allocator3> val_row;
1652  IVect ind_row;
1653  for (int i = 0 ; i < m ; i++)
1654  {
1655  n1 = A.GetRowSize(i);
1656  n2 = B.GetRowSize(i);
1657  size_row = n1 + n2;
1658  val_row.Reallocate(size_row);
1659  ind_row.Reallocate(size_row);
1660  for (int j = 0 ; j < n1 ; j++)
1661  {
1662  ind_row(j) = A.Index(i, j);
1663  val_row(j) = alpha*complex<T3>(A.Value(i, j), 0);
1664  }
1665 
1666  for (int j = 0 ; j < n2 ; j++)
1667  {
1668  ind_row(j+n1) = B.Index(i, j);
1669  val_row(j+n1) = alpha * complex<T3>(0, B.Value(i, j));
1670  }
1671 
1672  C.AddInteractionRow(i, size_row, ind_row, val_row);
1673  }
1674  }
1675 
1676 
1677  // C = C + alpha complex(A,B)
1678  template<class T0, class T1, class T2, class T3,
1679  class Allocator1, class Allocator2, class Allocator3>
1680  void AddMatrix(const T0& alpha,
1681  const Matrix<T1, Symmetric, ArrayRowSymSparse, Allocator1>& A,
1682  const Matrix<T2, Symmetric, ArrayRowSymSparse, Allocator2>& B,
1683  Matrix<complex<T3>, Symmetric, ArrayRowSymSparse, Allocator3>& C)
1684  {
1685  int m = B.GetM(), n1, n2, size_row;
1686  Vector<complex<T3>, VectFull, Allocator3> val_row;
1687  IVect ind_row;
1688  for (int i = 0 ; i < m ; i++)
1689  {
1690  n1 = A.GetRowSize(i);
1691  n2 = B.GetRowSize(i);
1692  size_row = n1 + n2;
1693  val_row.Reallocate(size_row);
1694  ind_row.Reallocate(size_row);
1695  for (int j = 0 ; j < n1 ; j++)
1696  {
1697  ind_row(j) = A.Index(i, j);
1698  val_row(j) = alpha * complex<T3>(A.Value(i, j), 0);
1699  }
1700 
1701  for (int j = 0 ; j < n2 ; j++)
1702  {
1703  ind_row(j+n1) = B.Index(i, j);
1704  val_row(j+n1) = alpha * complex<T3>(0, B.Value(i, j));
1705  }
1706 
1707  C.AddInteractionRow(i, size_row, ind_row, val_row);
1708  }
1709  }
1710 
1711 
1712  template<class T0, class T1, class T2, class Allocator1, class Allocator2>
1713  void AddMatrix(const T0& alpha,
1714  const Matrix<T1, Symmetric, ArrayRowSymSparse, Allocator1>& A,
1715  Matrix<T2, General, ArrayRowSparse, Allocator2>& B)
1716  {
1717  Vector<T2, VectFull, Allocator2> value;
1718  for (int i = 0; i < B.GetM(); i++)
1719  {
1720  int n = A.GetRowSize(i);
1721  value.Reallocate(n);
1722  for (int j = 0; j < n; j++)
1723  {
1724  value(j) = A.Value(i, j);
1725  if (A.Index(i, j) != i)
1726  B.AddInteraction(A.Index(i, j), i, alpha*value(j));
1727  }
1728 
1729  Mlt(alpha, value);
1730  B.AddInteractionRow(i, n, A.GetIndex(i), value.GetData());
1731  }
1732  }
1733 
1734 
1735  // Add //
1737 
1738 
1739 
1741  // Mlt //
1742 
1743 
1744  template<class T0, class T, class Allocator>
1745  void MltScalar(const T0& alpha, Matrix<T, General, ArrayRowSparse, Allocator>& A)
1746  {
1747  for (int i = 0; i < A.GetM(); i++)
1748  for (int j = 0; j < A.GetRowSize(i); j++)
1749  A.Value(i,j) *= alpha;
1750  }
1751 
1752 
1753  template<class T0, class T, class Allocator>
1754  void MltScalar(const T0& alpha, Matrix<T, General, ArrayColSparse, Allocator>& A)
1755  {
1756  for (int i = 0; i < A.GetN(); i++)
1757  for (int j = 0; j < A.GetColumnSize(i); j++)
1758  A.Value(i,j) *= alpha;
1759  }
1760 
1761 
1762  template<class T0, class T, class Allocator>
1763  void MltScalar(const T0& alpha,
1764  Matrix<T, Symmetric, ArrayRowSymSparse, Allocator>& A)
1765  {
1766  for (int i = 0; i < A.GetM(); i++)
1767  for (int j = 0; j < A.GetRowSize(i); j++)
1768  A.Value(i,j) *= alpha;
1769  }
1770 
1771 
1772  template<class T0, class T, class Allocator>
1773  void MltScalar(const T0& alpha,
1774  Matrix<T, Symmetric, ArrayColSymSparse, Allocator>& A)
1775  {
1776  for (int i = 0; i < A.GetN(); i++)
1777  for (int j = 0; j < A.GetColumnSize(i); j++)
1778  A.Value(i,j) *= alpha;
1779  }
1780 
1781 
1782  // Matrix-matrix product (sparse matrix against dense matrix)
1783  template<class T0, class T1, class Prop1, class Allocator1, class T4,
1784  class T2, class Prop2, class Storage2, class Allocator2,
1785  class T3, class Prop3, class Storage3, class Allocator3>
1786  void MltAddMatrix(const T0& alpha,
1787  const Matrix<T1, Prop1, ArrayRowSparse, Allocator1>& A,
1788  const Matrix<T2, Prop2, Storage2, Allocator2>& B,
1789  const T4& beta,
1790  Matrix<T3, Prop3, Storage3, Allocator3>& C)
1791  {
1792  if (Storage2::Sparse || Storage3::Sparse)
1793  throw WrongArgument("Mlt", "Function intended for product "
1794  " between a sparse matrix and a dense matrix");
1795 
1796 #ifdef SELDON_CHECK_DIMENSIONS
1797  CheckDim(A, B, "MltAdd(alpha, A, B, beta, C)");
1798 #endif
1799 
1800  T4 zero; SetComplexZero(zero);
1801  if (beta == zero)
1802  C.Fill(0);
1803  else
1804  Mlt(beta, C);
1805 
1806  int m = A.GetM();
1807  int n = B.GetN();
1808  T3 val;
1809  for (int i = 0; i < m; i++)
1810  for (int j = 0; j < n; j++)
1811  {
1812  SetComplexZero(val);
1813  for (int ind = 0; ind < A.GetRowSize(i); ind++)
1814  {
1815  int k = A.Index(i, ind);
1816  val += A.Value(i, ind) * B(k, j);
1817  }
1818 
1819  C(i, j) += alpha*val;
1820  }
1821  }
1822 
1823 
1824  // Matrix-matrix product (sparse matrix against dense matrix)
1825  template<class T0, class T1, class Prop1, class Storage1, class Allocator1,
1826  class T4, class T2, class Prop2, class Allocator2,
1827  class T3, class Prop3, class Storage3, class Allocator3>
1828  void MltAddMatrix(const T0& alpha,
1829  const Matrix<T1, Prop1, Storage1, Allocator1>& A,
1830  const Matrix<T2, Prop2, ArrayRowSparse, Allocator2>& B,
1831  const T4& beta,
1832  Matrix<T3, Prop3, Storage3, Allocator3>& C)
1833  {
1834  if (Storage1::Sparse || Storage3::Sparse)
1835  throw WrongArgument("Mlt", "Function intended for product "
1836  " between a sparse matrix and a dense matrix");
1837 
1838 #ifdef SELDON_CHECK_DIMENSIONS
1839  CheckDim(A, B, "MltAdd(alpha, A, B, beta, C)");
1840 #endif
1841 
1842  T4 zero; SetComplexZero(zero);
1843  if (beta == zero)
1844  C.Fill(0);
1845  else
1846  Mlt(beta, C);
1847 
1848  int m = A.GetM();
1849  T2 val;
1850  for (int j = 0; j < B.GetM(); j++)
1851  {
1852  for (int ind = 0; ind < B.GetRowSize(j); ind++)
1853  {
1854  int k = B.Index(j, ind);
1855  val = B.Value(j, ind);
1856  for (int i = 0; i < m; i++)
1857  C(i, k) += alpha*val*A(i, j);
1858  }
1859  }
1860  }
1861 
1862 
1863  // Mlt //
1865 
1866 
1868  // Norms //
1869 
1870 
1872 
1876  template <class T, class Prop, class Allocator>
1877  typename ClassComplexType<T>::Treal
1878  MaxAbs(const Matrix<T, Prop, ArrayRowSparse, Allocator>& A)
1879  {
1880  typename ClassComplexType<T>::Treal res(0);
1881  for (int i = 0; i < A.GetM(); i++)
1882  for (int j = 0; j < A.GetRowSize(i); j++)
1883  res = max(res, ComplexAbs(A.Value(i, j)));
1884 
1885  return res;
1886  }
1887 
1888 
1890 
1894  template <class T, class Prop, class Allocator>
1895  typename ClassComplexType<T>::Treal
1896  Norm1(const Matrix<T, Prop, ArrayRowSparse, Allocator>& A)
1897  {
1898  typedef typename ClassComplexType<T>::Treal Treal;
1899  Vector<Treal> sum(A.GetN());
1900  sum.Fill(Treal(0));
1901  for (int i = 0; i < A.GetM(); i++)
1902  for (int j = 0; j < A.GetRowSize(i); j++)
1903  sum(A.Index(i, j)) += ComplexAbs( A.Value(i, j));
1904 
1905  return sum.GetNormInf();
1906  }
1907 
1908 
1910 
1914  template <class T, class Prop, class Allocator>
1915  typename ClassComplexType<T>::Treal
1916  NormInf(const Matrix<T, Prop, ArrayRowSparse, Allocator>& A)
1917  {
1918  typedef typename ClassComplexType<T>::Treal Treal;
1919  Treal res(0), sum;
1920  for (int i = 0; i < A.GetM(); i++)
1921  {
1922  sum = Treal(0);
1923  for (int j = 0; j < A.GetRowSize(i); j++)
1924  sum += ComplexAbs(A.Value(i, j));
1925 
1926  res = max(res, sum);
1927  }
1928 
1929  return res;
1930  }
1931 
1932 
1934 
1938  template <class T, class Prop, class Allocator>
1939  typename ClassComplexType<T>::Treal
1940  MaxAbs(const Matrix<T, Prop, ArrayColSparse, Allocator>& A)
1941  {
1942  typename ClassComplexType<T>::Treal res(0);
1943  for (int i = 0; i < A.GetN(); i++)
1944  for (int j = 0; j < A.GetColumnSize(i); j++)
1945  res = max(res, ComplexAbs(A.Value(i, j)));
1946 
1947  return res;
1948  }
1949 
1950 
1952 
1956  template <class T, class Prop, class Allocator>
1957  typename ClassComplexType<T>::Treal
1958  Norm1(const Matrix<T, Prop, ArrayColSparse, Allocator>& A)
1959  {
1960  typedef typename ClassComplexType<T>::Treal Treal;
1961  Treal res(0), sum;
1962  for (int i = 0; i < A.GetN(); i++)
1963  {
1964  sum = Treal(0);
1965  for (int j = 0; j < A.GetColumnSize(i); j++)
1966  sum += ComplexAbs(A.Value(i, j));
1967 
1968  res = max(res, sum);
1969  }
1970 
1971  return res;
1972  }
1973 
1974 
1976 
1980  template <class T, class Prop, class Allocator>
1981  typename ClassComplexType<T>::Treal
1982  NormInf(const Matrix<T, Prop, ArrayColSparse, Allocator>& A)
1983  {
1984  typedef typename ClassComplexType<T>::Treal Treal;
1985  Vector<Treal> sum(A.GetM());
1986  sum.Fill(Treal(0));
1987  for (int i = 0; i < A.GetN(); i++)
1988  for (int j = 0; j < A.GetColumnSize(i); j++)
1989  sum(A.Index(i, j)) += ComplexAbs(A.Value(i, j));
1990 
1991  return sum.GetNormInf();
1992  }
1993 
1994 
1996 
2000  template <class T, class Prop, class Allocator>
2001  typename ClassComplexType<T>::Treal
2002  MaxAbs(const Matrix<T, Prop, ArrayRowSymSparse, Allocator>& A)
2003  {
2004  typename ClassComplexType<T>::Treal res(0);
2005  for (int i = 0; i < A.GetM(); i++)
2006  for (int j = 0; j < A.GetRowSize(i); j++)
2007  res = max(res, ComplexAbs(A.Value(i, j)));
2008 
2009  return res;
2010  }
2011 
2012 
2014 
2018  template <class T, class Prop, class Allocator>
2019  typename ClassComplexType<T>::Treal
2020  Norm1(const Matrix<T, Prop, ArrayRowSymSparse, Allocator>& A)
2021  {
2022  typedef typename ClassComplexType<T>::Treal Treal;
2023  Vector<Treal> sum(A.GetN());
2024  sum.Fill(Treal(0));
2025  for (int i = 0; i < A.GetM(); i++)
2026  for (int j = 0; j < A.GetRowSize(i); j++)
2027  {
2028  sum(A.Index(i, j)) += ComplexAbs( A.Value(i, j));
2029  if (A.Index(i, j) != i)
2030  sum(i) += ComplexAbs(A.Value(i, j));
2031  }
2032 
2033  return sum.GetNormInf();
2034  }
2035 
2036 
2038 
2042  template <class T, class Prop, class Allocator>
2043  typename ClassComplexType<T>::Treal
2044  NormInf(const Matrix<T, Prop, ArrayRowSymSparse, Allocator>& A)
2045  {
2046  return Norm1(A);
2047  }
2048 
2049 
2051 
2055  template <class T, class Prop, class Allocator>
2056  typename ClassComplexType<T>::Treal
2057  MaxAbs(const Matrix<T, Prop, ArrayColSymSparse, Allocator>& A)
2058  {
2059  typename ClassComplexType<T>::Treal res(0);
2060  for (int i = 0; i < A.GetM(); i++)
2061  for (int j = 0; j < A.GetColumnSize(i); j++)
2062  res = max(res, ComplexAbs(A.Value(i, j)));
2063 
2064  return res;
2065  }
2066 
2067 
2069 
2073  template <class T, class Prop, class Allocator>
2074  typename ClassComplexType<T>::Treal
2075  Norm1(const Matrix<T, Prop, ArrayColSymSparse, Allocator>& A)
2076  {
2077  typedef typename ClassComplexType<T>::Treal Treal;
2078  Vector<Treal> sum(A.GetN());
2079  sum.Fill(Treal(0));
2080  for (int i = 0; i < A.GetM(); i++)
2081  for (int j = 0; j < A.GetColumnSize(i); j++)
2082  {
2083  sum(A.Index(i, j)) += ComplexAbs( A.Value(i, j));
2084  if (A.Index(i, j) != i)
2085  sum(i) += ComplexAbs(A.Value(i, j));
2086  }
2087 
2088  return sum.GetNormInf();
2089  }
2090 
2091 
2093 
2097  template <class T, class Prop, class Allocator>
2098  typename ClassComplexType<T>::Treal
2099  NormInf(const Matrix<T, Prop, ArrayColSymSparse, Allocator>& A)
2100  {
2101  return Norm1(A);
2102  }
2103 
2104 
2106  template<class T, class Prop, class Allocator>
2107  typename ClassComplexType<T>::Treal
2108  NormFro(const Matrix<T, Prop, ArrayRowSparse, Allocator>& A)
2109  {
2110  typename ClassComplexType<T>::Treal res(0);
2111  for (int i = 0; i < A.GetM(); i++)
2112  for (int j = 0; j < A.GetRowSize(i); j++)
2113  res += absSquare(A.Value(i, j));
2114 
2115  res = sqrt(res);
2116  return res;
2117  }
2118 
2119 
2121  template<class T, class Prop, class Allocator>
2122  typename ClassComplexType<T>::Treal
2123  NormFro(const Matrix<T, Prop, ArrayRowSymSparse, Allocator>& A)
2124  {
2125  typename ClassComplexType<T>::Treal res(0);
2126  for (int i = 0; i < A.GetM(); i++)
2127  for (int j = 0; j < A.GetRowSize(i); j++)
2128  {
2129  if (i == A.Index(i, j))
2130  res += absSquare(A.Value(i, j));
2131  else
2132  res += 2.0*absSquare(A.Value(i, j));
2133  }
2134 
2135  res = sqrt(res);
2136  return res;
2137  }
2138 
2139 
2140  // Norms //
2142 
2143 
2145  // Transpose //
2146 
2147 
2149  template<class T, class Allocator>
2150  void Transpose(const Matrix<T, General, ArrayRowSparse, Allocator>& A,
2151  Matrix<T, General, ArrayRowSparse, Allocator>& B)
2152  {
2153  B.Clear();
2154 
2155  int m = A.GetM();
2156  int n = A.GetN();
2157  Vector<int> ptr_T(n);
2158 
2159  B.Reallocate(n, m);
2160 
2161  // For each column j, computes number of its non-zeroes and stores it in
2162  // ptr_T[j].
2163  ptr_T.Zero();
2164  for (int i = 0; i < m; i++)
2165  for (int j = 0; j < A.GetRowSize(i); j++)
2166  ptr_T(A.Index(i, j))++;
2167 
2168  for (int i = 0; i < n; i++)
2169  B.ReallocateRow(i, ptr_T(i));
2170 
2171  // filling matrix B
2172  ptr_T.Zero();
2173  for (int i = 0; i < m; i++)
2174  for (int jp = 0; jp < A.GetRowSize(i); jp++)
2175  {
2176  int j = A.Index(i, jp);
2177  int k = ptr_T(j);
2178  ++ptr_T(j);
2179  B.Value(j, k) = A.Value(i, jp);
2180  B.Index(j, k) = i;
2181  }
2182 
2183  B.Assemble();
2184  }
2185 
2186 
2188  template<class T, class Allocator>
2189  void Transpose(Matrix<T, General, ArrayRowSparse, Allocator>& A)
2190  {
2191  Matrix<T, General, ArrayRowSparse, Allocator> Acopy(A);
2192  Transpose(Acopy, A);
2193  }
2194 
2195 
2197  template<class T, class Allocator>
2198  void Transpose(const Matrix<T, General, ArrayColSparse, Allocator>& A,
2199  Matrix<T, General, ArrayColSparse, Allocator>& B)
2200  {
2201  B.Clear();
2202 
2203  int m = A.GetM();
2204  int n = A.GetN();
2205  Vector<int> ptr_T(m);
2206 
2207  B.Reallocate(n, m);
2208 
2209  // For each row j, computes number of its non-zeroes and stores it in
2210  // ptr_T[j].
2211  ptr_T.Zero();
2212  for (int i = 0; i < n; i++)
2213  for (int j = 0; j < A.GetColumnSize(i); j++)
2214  ptr_T(A.Index(i, j))++;
2215 
2216  for (int i = 0; i < m; i++)
2217  B.ReallocateColumn(i, ptr_T(i));
2218 
2219  // filling matrix B
2220  ptr_T.Zero();
2221  for (int i = 0; i < n; i++)
2222  for (int jp = 0; jp < A.GetColumnSize(i); jp++)
2223  {
2224  int j = A.Index(i, jp);
2225  int k = ptr_T(j);
2226  ++ptr_T(j);
2227  B.Value(j, k) = A.Value(i, jp);
2228  B.Index(j, k) = i;
2229  }
2230 
2231  B.Assemble();
2232  }
2233 
2234 
2236  template<class T, class Allocator>
2237  void Transpose(Matrix<T, General, ArrayColSparse, Allocator>& A)
2238  {
2239  Matrix<T, General, ArrayColSparse, Allocator> Acopy(A);
2240  Transpose(Acopy, A);
2241  }
2242 
2243 
2245  template<class T, class Allocator>
2246  void Conjugate(Matrix<T, General, ArrayRowSparse, Allocator>& A)
2247  {
2248  for (int i = 0; i < A.GetM(); i++)
2249  for (int j = 0; j < A.GetRowSize(i); j++)
2250  A.Value(i, j) = conjugate(A.Value(i, j));
2251  }
2252 
2253 
2255  template<class T, class Allocator>
2256  void Conjugate(Matrix<T, General, ArrayColSparse, Allocator>& A)
2257  {
2258  for (int i = 0; i < A.GetN(); i++)
2259  for (int j = 0; j < A.GetColumnSize(i); j++)
2260  A.Value(i, j) = conjugate(A.Value(i, j));
2261  }
2262 
2263 
2265  template<class T, class Allocator>
2266  void Conjugate(Matrix<T, Symmetric, ArrayRowSymSparse, Allocator>& A)
2267  {
2268  for (int i = 0; i < A.GetM(); i++)
2269  for (int j = 0; j < A.GetRowSize(i); j++)
2270  A.Value(i, j) = conjugate(A.Value(i, j));
2271  }
2272 
2273 
2275  template<class T, class Allocator>
2276  void Conjugate(Matrix<T, Symmetric, ArrayColSymSparse, Allocator>& A)
2277  {
2278  for (int i = 0; i < A.GetN(); i++)
2279  for (int j = 0; j < A.GetColumnSize(i); j++)
2280  A.Value(i, j) = conjugate(A.Value(i, j));
2281  }
2282 
2283 
2284  // Transpose //
2286 
2287 
2289  // MltAdd (matrix-matrix product) //
2290 
2291 
2292  template<class T1, class Prop1, class Allocator1,
2293  class T2, class Prop2, class Allocator2,
2294  class T4, class Prop4, class Allocator4>
2295  void MltMatrix(const Matrix<T1, Prop1, ArrayRowSparse, Allocator1>& A,
2296  const Matrix<T2, Prop2, ArrayRowSparse, Allocator2>& B,
2297  Matrix<T4, Prop4, ArrayRowSparse, Allocator4>& C)
2298  {
2299  T4 one, zero;
2300  SetComplexOne(one);
2301  SetComplexZero(zero);
2302  MltAdd(one, A, B, zero, C);
2303  }
2304 
2305 
2306  // C = beta * C + alpha * A * B
2307  template<class T0, class T1, class Prop1, class Allocator1,
2308  class T2, class Prop2, class Allocator2, class T3,
2309  class T4, class Prop4, class Allocator4>
2310  void MltAddMatrix(const T0& alpha,
2311  const Matrix<T1, Prop1, ArrayRowSparse, Allocator1>& A,
2312  const Matrix<T2, Prop2, ArrayRowSparse, Allocator2>& B,
2313  const T3& beta,
2314  Matrix<T4, Prop4, ArrayRowSparse, Allocator4>& C)
2315  {
2316  int m = A.GetM();
2317  int n = B.GetN();
2318 
2319  T3 zero;
2320  SetComplexZero(zero);
2321  if (beta == zero)
2322  {
2323  C.Clear();
2324  C.Reallocate(m, n);
2325  }
2326  else
2327  Mlt(beta, C);
2328 
2329 #ifdef SELDON_CHECK_DIMENSIONS
2330  CheckDim(A, B, C, "MltAdd(alpha, A, B, beta, C)");
2331 #endif
2332 
2333  Vector<int> Index(n), IndCol(n);
2334  Vector<T4, VectFull, Allocator4> Value(n);
2335  Index.Fill(-1);
2336  int col, ind; T4 vloc;
2337  // loop over rows of matrix A
2338  for (int i = 0; i < m; i++)
2339  {
2340  int nnz = 0;
2341  for (int j = 0; j < A.GetRowSize(i); j++)
2342  {
2343  col = A.Index(i, j);
2344  // for each non-zero entry of the row i of A
2345  // loop over non-zeros entries of the corresponding row of B
2346  for (int k = 0; k < B.GetRowSize(col); k++)
2347  {
2348  ind = B.Index(col, k);
2349  vloc = alpha*A.Value(i, j)*B.Value(col, k);
2350  if (Index(ind) >= 0)
2351  {
2352  // already existing entry, we add it
2353  Value(Index(ind)) += vloc;
2354  }
2355  else
2356  {
2357  // new non-zero entry
2358  Value(nnz) = vloc;
2359  IndCol(nnz) = ind;
2360  Index(ind) = nnz++;
2361  }
2362  }
2363  }
2364 
2365  // sorting entries
2366  Sort(nnz, IndCol, Value);
2367 
2368  // adding interactions to matrix C
2369  C.AddInteractionRow(i, nnz, IndCol, Value, true);
2370 
2371  // resetting Index
2372  for (int j = 0; j < nnz; j++)
2373  Index(IndCol(j)) = -1;
2374  }
2375  }
2376 
2377 
2378  // C = beta * C + alpha * A * B
2379  template<class T0, class T1, class Prop1, class Allocator1,
2380  class T2, class Prop2, class Allocator2, class T3,
2381  class T4, class Prop4, class Allocator4>
2382  void MltAddMatrix(const T0& alpha,
2383  const SeldonTranspose& TransA,
2384  const Matrix<T1, Prop1, ArrayRowSparse, Allocator1>& A,
2385  const SeldonTranspose& TransB,
2386  const Matrix<T2, Prop2, ArrayRowSparse, Allocator2>& B,
2387  const T3& beta,
2388  Matrix<T4, Prop4, ArrayRowSparse, Allocator4>& C)
2389  {
2390  int m = A.GetM();
2391  int n = B.GetN();
2392  if (!TransA.NoTrans())
2393  m = A.GetN();
2394 
2395  if (!TransB.NoTrans())
2396  n = B.GetM();
2397 
2398  T3 zero, one;
2399  SetComplexZero(zero);
2400  SetComplexOne(one);
2401  if (beta == zero)
2402  {
2403  C.Clear();
2404  C.Reallocate(m, n);
2405  }
2406  else
2407  {
2408  if ((!TransA.NoTrans()) || (!TransB.NoTrans()))
2409  Mlt(beta, C);
2410  }
2411 
2412 #ifdef SELDON_CHECK_DIMENSIONS
2413  CheckDim(TransA, A, TransB, B, C, "MltAdd(alpha, A, B, beta, C)");
2414 #endif
2415 
2416  int col, ind; T4 vloc;
2417 
2418  if (TransA.NoTrans())
2419  {
2420  if (TransB.NoTrans())
2421  MltAdd(alpha, A, B, beta, C);
2422  else if (TransB.Trans())
2423  {
2424  Vector<int> Index(n), IndCol(n);
2425  Vector<T4, VectFull, Allocator4> Value(n);
2426  Index.Fill(-1);
2427 
2428  // loop over rows of matrix A
2429  for (int i = 0; i < A.GetM(); i++)
2430  {
2431  int nnz = 0;
2432  for (int j = 0; j < A.GetRowSize(i); j++)
2433  {
2434  col = A.Index(i, j);
2435  // for each non-zero entry of the row i of A
2436  // loop over non-zeros entries of the corresponding row of B
2437  for (int ib = 0; ib < B.GetM(); ib++)
2438  for (int k = 0; k < B.GetRowSize(ib); k++)
2439  if (B.Index(ib, k) == col)
2440  {
2441  ind = ib;
2442  vloc = alpha*A.Value(i, j)*B.Value(ib, k);
2443  if (Index(ind) >= 0)
2444  {
2445  // already existing entry, we add it
2446  Value(Index(ind)) += vloc;
2447  }
2448  else
2449  {
2450  // new non-zero entry
2451  Value(nnz) = vloc;
2452  IndCol(nnz) = ind;
2453  Index(ind) = nnz++;
2454  }
2455  }
2456  }
2457 
2458  // sorting entries
2459  Sort(nnz, IndCol, Value);
2460 
2461  // adding interactions to matrix C
2462  C.AddInteractionRow(i, nnz, IndCol, Value, true);
2463 
2464  // resetting Index
2465  for (int j = 0; j < nnz; j++)
2466  Index(IndCol(j)) = -1;
2467  }
2468  }
2469  else
2470  {
2471  Vector<int> Index(n), IndCol(n);
2472  Vector<T4, VectFull, Allocator4> Value(n);
2473  Index.Fill(-1);
2474 
2475  // loop over rows of matrix A
2476  for (int i = 0; i < A.GetM(); i++)
2477  {
2478  int nnz = 0;
2479  for (int j = 0; j < A.GetRowSize(i); j++)
2480  {
2481  col = A.Index(i, j);
2482  // for each non-zero entry of the row i of A
2483  // loop over non-zeros entries of the corresponding row of B
2484  for (int ib = 0; ib < B.GetM(); ib++)
2485  for (int k = 0; k < B.GetRowSize(ib); k++)
2486  if (B.Index(ib, k) == col)
2487  {
2488  ind = ib;
2489  vloc = alpha*A.Value(i, j)*conjugate(B.Value(ib, k));
2490  if (Index(ind) >= 0)
2491  {
2492  // already existing entry, we add it
2493  Value(Index(ind)) += vloc;
2494  }
2495  else
2496  {
2497  // new non-zero entry
2498  Value(nnz) = vloc;
2499  IndCol(nnz) = ind;
2500  Index(ind) = nnz++;
2501  }
2502  }
2503  }
2504 
2505  // sorting entries
2506  Sort(nnz, IndCol, Value);
2507 
2508  // adding interactions to matrix C
2509  C.AddInteractionRow(i, nnz, IndCol, Value, true);
2510 
2511  // resetting Index
2512  for (int j = 0; j < nnz; j++)
2513  Index(IndCol(j)) = -1;
2514  }
2515  }
2516  }
2517  else if (TransA.Trans())
2518  {
2519  if (TransB.NoTrans())
2520  {
2521  // C_{ind, col} = C_{ind, col} + \sum_i A_{i, ind} B_{i, col}
2522  for (int i = 0; i < B.GetM(); i++)
2523  {
2524  for (int j = 0; j < B.GetRowSize(i); j++)
2525  {
2526  col = B.Index(i, j);
2527  for (int k = 0; k < A.GetRowSize(i); k++)
2528  {
2529  ind = A.Index(i, k);
2530  vloc = alpha*B.Value(i, j)*A.Value(i, k);
2531  C.AddInteraction(ind, col, vloc);
2532  }
2533  }
2534  }
2535  }
2536  else if (TransB.Trans())
2537  {
2538  Vector<int> Index(m), IndCol(m);
2539  Vector<T4, VectFull, Allocator4> Value(m);
2540  Index.Fill(-1);
2541 
2542  // loop over rows of matrix B
2543  for (int i = 0; i < B.GetM(); i++)
2544  {
2545  int nnz = 0;
2546  for (int j = 0; j < B.GetRowSize(i); j++)
2547  {
2548  col = B.Index(i, j);
2549  // for each non-zero entry of the row i of B
2550  // loop over non-zeros entries of the corresponding row of A
2551  for (int k = 0; k < A.GetRowSize(col); k++)
2552  {
2553  ind = A.Index(col, k);
2554  vloc = alpha*B.Value(i, j)*A.Value(col, k);
2555  if (Index(ind) >= 0)
2556  {
2557  // already existing entry, we add it
2558  Value(Index(ind)) += vloc;
2559  }
2560  else
2561  {
2562  // new non-zero entry
2563  Value(nnz) = vloc;
2564  IndCol(nnz) = ind;
2565  Index(ind) = nnz++;
2566  }
2567  }
2568  }
2569 
2570  // sorting entries
2571  Sort(nnz, IndCol, Value);
2572 
2573  // adding interactions to matrix C
2574  C.AddInteractionColumn(i, nnz, IndCol, Value);
2575 
2576  // resetting Index
2577  for (int j = 0; j < nnz; j++)
2578  Index(IndCol(j)) = -1;
2579  }
2580  }
2581  else
2582  {
2583  Vector<int> Index(m), IndCol(m);
2584  Vector<T4, VectFull, Allocator4> Value(m);
2585  Index.Fill(-1);
2586 
2587  // loop over rows of matrix B
2588  for (int i = 0; i < B.GetM(); i++)
2589  {
2590  int nnz = 0;
2591  for (int j = 0; j < B.GetRowSize(i); j++)
2592  {
2593  col = B.Index(i, j);
2594  // for each non-zero entry of the row i of B
2595  // loop over non-zeros entries of the corresponding row of A
2596  for (int k = 0; k < A.GetRowSize(col); k++)
2597  {
2598  ind = A.Index(col, k);
2599  vloc = alpha*conjugate(B.Value(i, j))*A.Value(col, k);
2600  if (Index(ind) >= 0)
2601  {
2602  // already existing entry, we add it
2603  Value(Index(ind)) += vloc;
2604  }
2605  else
2606  {
2607  // new non-zero entry
2608  Value(nnz) = vloc;
2609  IndCol(nnz) = ind;
2610  Index(ind) = nnz++;
2611  }
2612  }
2613  }
2614 
2615  // sorting entries
2616  Sort(nnz, IndCol, Value);
2617 
2618  // adding interactions to matrix C
2619  C.AddInteractionColumn(i, nnz, IndCol, Value);
2620 
2621  // resetting Index
2622  for (int j = 0; j < nnz; j++)
2623  Index(IndCol(j)) = -1;
2624  }
2625  }
2626  }
2627  else if (TransA.ConjTrans())
2628  {
2629  if (TransB.NoTrans())
2630  {
2631  // C_{ind, col} = C_{ind, col} + \sum_i A_{i, ind} B_{i, col}
2632  for (int i = 0; i < B.GetM(); i++)
2633  {
2634  for (int j = 0; j < B.GetRowSize(i); j++)
2635  {
2636  col = B.Index(i, j);
2637  for (int k = 0; k < A.GetRowSize(i); k++)
2638  {
2639  ind = A.Index(i, k);
2640  vloc = alpha*B.Value(i, j)*conjugate(A.Value(i, k));
2641  C.AddInteraction(ind, col, vloc);
2642  }
2643  }
2644  }
2645  }
2646  else if (TransB.Trans())
2647  {
2648  Vector<int> Index(m), IndCol(m);
2649  Vector<T4, VectFull, Allocator4> Value(m);
2650  Index.Fill(-1);
2651 
2652  // loop over rows of matrix B
2653  for (int i = 0; i < B.GetM(); i++)
2654  {
2655  int nnz = 0;
2656  for (int j = 0; j < B.GetRowSize(i); j++)
2657  {
2658  col = B.Index(i, j);
2659  // for each non-zero entry of the row i of B
2660  // loop over non-zeros entries of the corresponding row of A
2661  for (int k = 0; k < A.GetRowSize(col); k++)
2662  {
2663  ind = A.Index(col, k);
2664  vloc = alpha*B.Value(i, j)*conjugate(A.Value(col, k));
2665  if (Index(ind) >= 0)
2666  {
2667  // already existing entry, we add it
2668  Value(Index(ind)) += vloc;
2669  }
2670  else
2671  {
2672  // new non-zero entry
2673  Value(nnz) = vloc;
2674  IndCol(nnz) = ind;
2675  Index(ind) = nnz++;
2676  }
2677  }
2678  }
2679 
2680  // sorting entries
2681  Sort(nnz, IndCol, Value);
2682 
2683  // adding interactions to matrix C
2684  C.AddInteractionColumn(i, nnz, IndCol, Value);
2685 
2686  // resetting Index
2687  for (int j = 0; j < nnz; j++)
2688  Index(IndCol(j)) = -1;
2689  }
2690  }
2691  else
2692  {
2693  Vector<int> Index(m), IndCol(m);
2694  Vector<T4, VectFull, Allocator4> Value(m);
2695  Index.Fill(-1);
2696 
2697  // loop over rows of matrix B
2698  for (int i = 0; i < B.GetM(); i++)
2699  {
2700  int nnz = 0;
2701  for (int j = 0; j < B.GetRowSize(i); j++)
2702  {
2703  col = B.Index(i, j);
2704  // for each non-zero entry of the row i of B
2705  // loop over non-zeros entries of the corresponding row of A
2706  for (int k = 0; k < A.GetRowSize(col); k++)
2707  {
2708  ind = A.Index(col, k);
2709  vloc = alpha*conjugate(B.Value(i, j))*conjugate(A.Value(col, k));
2710  if (Index(ind) >= 0)
2711  {
2712  // already existing entry, we add it
2713  Value(Index(ind)) += vloc;
2714  }
2715  else
2716  {
2717  // new non-zero entry
2718  Value(nnz) = vloc;
2719  IndCol(nnz) = ind;
2720  Index(ind) = nnz++;
2721  }
2722  }
2723  }
2724 
2725  // sorting entries
2726  Sort(nnz, IndCol, Value);
2727 
2728  // adding interactions to matrix C
2729  C.AddInteractionColumn(i, nnz, IndCol, Value);
2730 
2731  // resetting Index
2732  for (int j = 0; j < nnz; j++)
2733  Index(IndCol(j)) = -1;
2734  }
2735  }
2736  }
2737  }
2738 
2739 
2740  // MltAdd (matrix-matrix product) //
2742 
2743 
2745  template<class T, class Prop, class Storage, class Allocator, class T0>
2746  void RemoveSmallEntry(Matrix<T, Prop, Storage, Allocator>& A,
2747  const T0& epsilon)
2748  {
2749  A.RemoveSmallEntry(epsilon);
2750  }
2751 
2752 
2754  template<class T, class Prop, class Allocator, class T0>
2755  void RemoveSmallEntry(Matrix<T, Prop, RowSparse, Allocator>& A,
2756  const T0& epsilon)
2757  {
2758  // TO BE DONE
2759  }
2760 
2761 
2763  template<class T, class Prop, class Allocator, class T0>
2764  void RemoveSmallEntry(Matrix<T, Prop, RowSymSparse, Allocator>& A,
2765  const T0& epsilon)
2766  {
2767  // TO BE DONE
2768  }
2769 
2770 
2772 
2776  template<class T1, class Prop, class Allocator>
2777  void EraseCol(const IVect& col_number,
2778  Matrix<T1, Prop, ArrayRowSymSparse, Allocator>& A)
2779  {
2780  int m = col_number.GetM();
2781  // index array to know fastly if it is a column to erase
2782  IVect index(A.GetM()); index.Fill(-1);
2783  for (int i = 0; i < m; i++)
2784  index(col_number(i)) = i;
2785 
2786  // first, we remove rows
2787  for (int i = 0; i < A.GetM(); i++)
2788  {
2789  if (index(i) != -1)
2790  A.ClearRow(i);
2791  }
2792 
2793  // then columns
2794  for (int i = 0; i < A.GetM(); i++)
2795  {
2796  bool something_to_remove = false;
2797  for (int j = 0; j < A.GetRowSize(i); j++)
2798  if (index(A.Index(i,j)) != -1)
2799  something_to_remove = true;
2800 
2801  if (something_to_remove)
2802  {
2803  int nb = 0;
2804  for (int j = 0; j < A.GetRowSize(i); j++)
2805  if (index(A.Index(i,j)) == -1)
2806  {
2807  A.Index(i, nb) = A.Index(i, j);
2808  A.Value(i, nb) = A.Value(i, j);
2809  nb++;
2810  }
2811 
2812  A.ResizeRow(i, nb);
2813  }
2814  }
2815  }
2816 
2817 
2819 
2823  template<class T1, class Prop, class Allocator>
2824  void EraseCol(const IVect& col_number,
2825  Matrix<T1, Prop, ArrayRowSparse, Allocator>& A)
2826  {
2827  int m = col_number.GetM();
2828  // index array to know fastly if it is a column to erase
2829  IVect index(A.GetM()); index.Fill(-1);
2830  for (int i = 0; i < m; i++)
2831  index(col_number(i)) = i;
2832 
2833  for (int i = 0; i < A.GetM(); i++)
2834  {
2835  bool something_to_remove = false;
2836  for (int j = 0; j < A.GetRowSize(i); j++)
2837  if (index(A.Index(i,j)) != -1)
2838  something_to_remove = true;
2839 
2840  if (something_to_remove)
2841  {
2842  int nb = 0;
2843  for (int j = 0; j < A.GetRowSize(i); j++)
2844  if (index(A.Index(i,j)) == -1)
2845  {
2846  A.Index(i, nb) = A.Index(i, j);
2847  A.Value(i, nb) = A.Value(i, j);
2848  nb++;
2849  }
2850 
2851  A.ResizeRow(i, nb);
2852  }
2853  }
2854  }
2855 
2856 
2858 
2862  template<class T1, class Prop, class Allocator>
2863  void EraseCol(const IVect& col_number,
2864  Matrix<T1, Prop, RowSparse, Allocator>& A)
2865  {
2866  int m = A.GetM(), n = A.GetN();
2867  long nnz = A.GetIndSize();
2868  long* ptr = A.GetPtr();
2869  int* ind = A.GetInd();
2870  T1* data = A.GetData();
2871  Vector<bool> ColToKeep(n);
2872  ColToKeep.Fill(true);
2873  for (int i = 0; i < col_number.GetM(); i++)
2874  ColToKeep(col_number(i)) = false;
2875 
2876  for (long i = 0; i < A.GetIndSize(); i++)
2877  if (!ColToKeep(ind[i]))
2878  nnz--;
2879 
2880  if (nnz == A.GetIndSize())
2881  return;
2882 
2883  Vector<long> Ptr(m+1); Vector<int> Ind(nnz);
2884  Vector<T1, VectFull, Allocator> Val(nnz);
2885  Ptr(0) = 0;
2886  for (int i = 0; i < m; i++)
2887  {
2888  long jA = Ptr(i); int size_row = 0;
2889  for (long j = ptr[i]; j < ptr[i+1]; j++)
2890  if (ColToKeep(ind[j]))
2891  {
2892  Ind(jA) = ind[j];
2893  Val(jA) = data[j];
2894  size_row++; jA++;
2895  }
2896 
2897  Ptr(i+1) = Ptr(i) + size_row;
2898  }
2899 
2900  A.SetData(m, n, Val, Ptr, Ind);
2901  }
2902 
2903 
2905 
2909  template<class T1, class Prop, class Allocator>
2910  void EraseCol(const IVect& col_number,
2911  Matrix<T1, Prop, RowSymSparse, Allocator>& A)
2912  {
2913  int m = A.GetM(), n = A.GetN();
2914  long nnz = A.GetIndSize();
2915  long* ptr = A.GetPtr();
2916  int* ind = A.GetInd();
2917  T1* data = A.GetData();
2918  Vector<bool> ColToKeep(n);
2919  ColToKeep.Fill(true);
2920  for (int i = 0; i < col_number.GetM(); i++)
2921  ColToKeep(col_number(i)) = false;
2922 
2923  for (int i = 0; i < m; i++)
2924  {
2925  if (!ColToKeep(i))
2926  nnz -= ptr[i+1] - ptr[i];
2927  else
2928  {
2929  for (long j = ptr[i]; j < ptr[i+1]; j++)
2930  if (!ColToKeep(ind[j]))
2931  nnz--;
2932  }
2933  }
2934 
2935  if (nnz == A.GetIndSize())
2936  return;
2937 
2938  Vector<long> Ptr(m+1); Vector<int> Ind(nnz);
2939  Vector<T1, VectFull, Allocator> Val(nnz);
2940  Ptr(0) = 0;
2941  for (int i = 0; i < m; i++)
2942  {
2943  long jA = Ptr(i); int size_row = 0;
2944  if (ColToKeep(i))
2945  for (long j = ptr[i]; j < ptr[i+1]; j++)
2946  if (ColToKeep(ind[j]))
2947  {
2948  Ind(jA) = ind[j];
2949  Val(jA) = data[j];
2950  size_row++; jA++;
2951  }
2952 
2953  Ptr(i+1) = Ptr(i) + size_row;
2954  }
2955 
2956  A.SetData(m, n, Val, Ptr, Ind);
2957  }
2958 
2959 
2961 
2965  template<class T1, class Prop, class Storage, class Allocator>
2966  void EraseCol(const IVect& col_number,
2967  Matrix<T1, Prop, Storage, Allocator>& A)
2968  {
2969  cout << "Not implemented for any matrix" << endl;
2970  abort();
2971  }
2972 
2973 
2975 
2979  template<class T1, class Prop, class Storage, class Allocator>
2980  void EraseRow(const IVect& col_number,
2981  Matrix<T1, Prop, Storage, Allocator>& A)
2982  {
2983  cout << "Not implemented for any matrix" << endl;
2984  abort();
2985  }
2986 
2987 
2989 
2993  template<class T1, class Prop, class Allocator>
2994  void EraseRow(const IVect& col_number,
2995  Matrix<T1, Prop, ArrayRowSymSparse, Allocator>& A)
2996  {
2997  EraseCol(col_number, A);
2998  }
2999 
3000 
3002 
3006  template<class T1, class Prop, class Allocator>
3007  void EraseRow(const IVect& col_number,
3008  Matrix<T1, Prop, ArrayRowSparse, Allocator>& A)
3009  {
3010  for (int i = 0; i < col_number.GetM(); i++)
3011  A.ClearRow(col_number(i));
3012  }
3013 
3014 
3016 
3020  template<class T1, class Prop, class Allocator>
3021  void EraseRow(const IVect& col_number,
3022  Matrix<T1, Prop, RowSparse, Allocator>& A)
3023  {
3024  int m = A.GetM(), n = A.GetN();
3025  long nnz = A.GetIndSize();
3026  long* ptr = A.GetPtr();
3027  int* ind = A.GetInd();
3028  T1* data = A.GetData();
3029  Vector<bool> RowToKeep(m);
3030  RowToKeep.Fill(true);
3031  for (int i = 0; i < col_number.GetM(); i++)
3032  RowToKeep(col_number(i)) = false;
3033 
3034  for (int i = 0; i < m; i++)
3035  if (!RowToKeep(i))
3036  nnz -= ptr[i+1] - ptr[i];
3037 
3038  Vector<long> Ptr(m+1); Vector<int> Ind(nnz);
3039  Vector<T1, VectFull, Allocator> Val(nnz);
3040  Ptr(0) = 0;
3041  for (int i = 0; i < m; i++)
3042  {
3043  if (RowToKeep(i))
3044  {
3045  int size_row = ptr[i+1] - ptr[i];
3046  for (int j = 0; j < size_row; j++)
3047  {
3048  Ind(Ptr(i) + j) = ind[ptr[i] + j];
3049  Val(Ptr(i) + j) = data[ptr[i] + j];
3050  }
3051 
3052  Ptr(i+1) = Ptr(i) + size_row;
3053  }
3054  else
3055  Ptr(i+1) = Ptr(i);
3056  }
3057 
3058  A.SetData(m, n, Val, Ptr, Ind);
3059  }
3060 
3062 
3066  template<class T1, class Prop, class Allocator>
3067  void EraseRow(const IVect& col_number,
3068  Matrix<T1, Prop, RowSymSparse, Allocator>& A)
3069  {
3070  EraseCol(col_number, A);
3071  }
3072 
3073 
3075 
3080  template<class T, class Complexe, class Allocator>
3081  void GetRowSum(Vector<T>& diagonal_scale_left,
3082  const Matrix<Complexe, General,
3083  ArrayRowSparse, Allocator> & mat_direct)
3084  {
3085  int n = mat_direct.GetM();
3086  diagonal_scale_left.Reallocate(n);
3087  diagonal_scale_left.Fill(0);
3088  for (int i = 0; i < n; i++)
3089  for (int j = 0; j < mat_direct.GetRowSize(i); j++)
3090  diagonal_scale_left(i) += abs(mat_direct.Value(i,j));
3091 
3092  }
3093 
3094 
3096 
3101  template<class T, class Complexe, class Allocator>
3102  void GetRowSum(Vector<T>& diagonal_scale_left,
3103  const Matrix<Complexe, Symmetric,
3104  ArrayRowSymSparse, Allocator> & mat_direct)
3105  {
3106  int n = mat_direct.GetM();
3107  diagonal_scale_left.Reallocate(n);
3108  diagonal_scale_left.Fill(0);
3109  for (int i = 0; i < n; i++)
3110  for (int j = 0; j < mat_direct.GetRowSize(i); j++)
3111  {
3112  diagonal_scale_left(i) += abs(mat_direct.Value(i,j));
3113  if (i != mat_direct.Index(i,j))
3114  diagonal_scale_left(mat_direct.Index(i,j))
3115  += abs(mat_direct.Value(i,j));
3116  }
3117  }
3118 
3119 
3121 
3126  template<class T, class Complexe, class Allocator>
3127  void GetRowSum(Vector<T>& diagonal_scale_left,
3128  const Matrix<Complexe, General,
3129  RowSparse, Allocator> & mat_direct)
3130  {
3131  int n = mat_direct.GetM();
3132  diagonal_scale_left.Reallocate(n);
3133  diagonal_scale_left.Fill(0);
3134  long* ptr = mat_direct.GetPtr();
3135  Complexe* data = mat_direct.GetData();
3136  for (int i = 0; i < n; i++)
3137  for (long j = ptr[i]; j < ptr[i+1]; j++)
3138  diagonal_scale_left(i) += abs(data[j]);
3139 
3140  }
3141 
3142 
3144 
3149  template<class T, class Complexe, class Allocator>
3150  void GetRowSum(Vector<T>& diagonal_scale_left,
3151  const Matrix<Complexe, Symmetric,
3152  RowSymSparse, Allocator> & mat_direct)
3153  {
3154  int n = mat_direct.GetM();
3155  diagonal_scale_left.Reallocate(n);
3156  diagonal_scale_left.Fill(0);
3157  long* ptr = mat_direct.GetPtr();
3158  int* ind = mat_direct.GetInd();
3159  Complexe* data = mat_direct.GetData();
3160  for (int i = 0; i < n; i++)
3161  for (long j = ptr[i]; j < ptr[i+1]; j++)
3162  {
3163  diagonal_scale_left(i) += abs(data[j]);
3164  if (i != ind[j])
3165  diagonal_scale_left(ind[j]) += abs(data[j]);
3166  }
3167  }
3168 
3169 
3171 
3176  template<class T, class Complexe, class Allocator>
3177  void GetColSum(Vector<T>& diagonal_scale,
3178  const Matrix<Complexe, General,
3179  ArrayRowSparse, Allocator> & mat_direct)
3180  {
3181  int n = mat_direct.GetM();
3182  diagonal_scale.Reallocate(mat_direct.GetN());
3183  diagonal_scale.Fill(0);
3184  for (int i = 0; i < n; i++)
3185  for (int j = 0; j < mat_direct.GetRowSize(i); j++)
3186  diagonal_scale(mat_direct.Index(i, j)) += abs(mat_direct.Value(i, j));
3187 
3188  }
3189 
3190 
3192 
3197  template<class T, class Complexe, class Allocator>
3198  void GetColSum(Vector<T>& diagonal_scale,
3199  const Matrix<Complexe, Symmetric,
3200  ArrayRowSymSparse, Allocator> & mat_direct)
3201  {
3202  GetRowSum(diagonal_scale, mat_direct);
3203  }
3204 
3205 
3207 
3212  template<class T, class Complexe, class Allocator>
3213  void GetColSum(Vector<T>& diagonal_scale,
3214  const Matrix<Complexe, General,
3215  RowSparse, Allocator> & mat_direct)
3216  {
3217  int n = mat_direct.GetM();
3218  diagonal_scale.Reallocate(mat_direct.GetN());
3219  diagonal_scale.Fill(0);
3220  long* ptr = mat_direct.GetPtr();
3221  int* ind = mat_direct.GetInd();
3222  Complexe* data = mat_direct.GetData();
3223  for (int i = 0; i < n; i++)
3224  for (long j = ptr[i]; j < ptr[i+1]; j++)
3225  diagonal_scale(ind[j]) += abs(data[j]);
3226 
3227  }
3228 
3229 
3231 
3236  template<class T, class Complexe, class Allocator>
3237  void GetColSum(Vector<T>& diagonal_scale,
3238  const Matrix<Complexe, Symmetric,
3239  RowSymSparse, Allocator> & mat_direct)
3240  {
3241  GetRowSum(diagonal_scale, mat_direct);
3242  }
3243 
3244 
3247 
3254  template<class T, class Complexe, class Allocator>
3255  void GetRowColSum(Vector<T>& sum_row,
3256  Vector<T>& sum_col,
3257  const Matrix<Complexe, General,
3258  ArrayRowSparse, Allocator> & A)
3259  {
3260  int n = A.GetM();
3261  sum_row.Reallocate(n);
3262  sum_col.Reallocate(A.GetN());
3263  sum_row.Fill(0);
3264  sum_col.Fill(0);
3265  for (int i = 0; i < n; i++)
3266  for (int j = 0; j < A.GetRowSize(i); j++)
3267  {
3268  sum_row(i) += abs(A.Value(i,j));
3269  sum_col(A.Index(i, j)) += abs(A.Value(i,j));
3270  }
3271  }
3272 
3273 
3276 
3283  template<class T, class Complexe, class Allocator>
3284  void GetRowColSum(Vector<T>& sum_row,
3285  Vector<T>& sum_col,
3286  const Matrix<Complexe, Symmetric,
3287  ArrayRowSymSparse, Allocator> & A)
3288  {
3289  GetRowSum(sum_row, A);
3290  sum_col = sum_row;
3291  }
3292 
3293 
3296 
3303  template<class T, class Complexe, class Allocator>
3304  void GetRowColSum(Vector<T>& sum_row,
3305  Vector<T>& sum_col,
3306  const Matrix<Complexe, General,
3307  RowSparse, Allocator> & A)
3308  {
3309  int n = A.GetM();
3310  sum_row.Reallocate(n);
3311  sum_col.Reallocate(A.GetN());
3312  sum_row.Fill(0);
3313  sum_col.Fill(0);
3314  long* ptr = A.GetPtr();
3315  int* ind = A.GetInd();
3316  Complexe* data = A.GetData();
3317  for (int i = 0; i < n; i++)
3318  for (long j = ptr[i]; j < ptr[i+1]; j++)
3319  {
3320  sum_row(i) += abs(data[j]);
3321  sum_col(ind[j]) += abs(data[j]);
3322  }
3323  }
3324 
3325 
3328 
3335  template<class T, class Complexe, class Allocator>
3336  void GetRowColSum(Vector<T>& sum_row,
3337  Vector<T>& sum_col,
3338  const Matrix<Complexe, Symmetric,
3339  RowSymSparse, Allocator> & A)
3340  {
3341  GetRowSum(sum_row, A);
3342  sum_col = sum_row;
3343  }
3344 
3345 
3347  template<class T0, class Prop0, class Allocator0,
3348  class T1, class Allocator1>
3349  void CopySubMatrix(const Matrix<T0, Prop0, ArrayRowSparse, Allocator0>& A,
3350  const IVect& row, const IVect& col,
3351  Vector<int>& RowNum,
3352  Vector<int>& ColNum,
3353  Vector<T1, VectFull, Allocator1>& Value)
3354  {
3355  int m = A.GetM(), n = A.GetN();
3356  if ((m <= 0) || (n <= 0) || (row.GetM() <= 0) || (col.GetM() <= 0))
3357  {
3358  RowNum.Clear(); ColNum.Clear(); Value.Clear();
3359  return;
3360  }
3361 
3362  Vector<bool> RowKept(m), ColKept(n);
3363  RowKept.Fill(false); ColKept.Fill(false);
3364  for (int i = 0; i < row.GetM(); i++)
3365  RowKept(row(i)) = true;
3366 
3367  for (int i = 0; i < col.GetM(); i++)
3368  ColKept(col(i)) = true;
3369 
3370  long nnz = 0;
3371  for (int i = 0; i < A.GetM(); i++)
3372  if (RowKept(i))
3373  for (int j = 0; j < A.GetRowSize(i); j++)
3374  if (ColKept(A.Index(i, j)))
3375  nnz++;
3376 
3377  RowNum.Reallocate(nnz);
3378  ColNum.Reallocate(nnz);
3379  Value.Reallocate(nnz);
3380  nnz = 0;
3381  for (int i = 0; i < A.GetM(); i++)
3382  if (RowKept(i))
3383  for (int j = 0; j < A.GetRowSize(i); j++)
3384  if (ColKept(A.Index(i, j)))
3385  {
3386  RowNum(nnz) = i;
3387  ColNum(nnz) = A.Index(i, j);
3388  Value(nnz) = A.Value(i, j);
3389  nnz++;
3390  }
3391  }
3392 
3393 
3395  template<class T0, class Prop0, class Allocator0,
3396  class T1, class Allocator1>
3397  void CopySubMatrix(const Matrix<T0, Prop0,
3398  ArrayRowSymSparse, Allocator0>& A,
3399  const IVect& row, const IVect& col,
3400  Vector<int>& RowNum,
3401  Vector<int>& ColNum,
3402  Vector<T1, VectFull, Allocator1>& Value)
3403  {
3404  int m = A.GetM(), n = A.GetN();
3405  if ((m <= 0) || (n <= 0) || (row.GetM() <= 0) || (col.GetM() <= 0))
3406  {
3407  RowNum.Clear(); ColNum.Clear(); Value.Clear();
3408  return;
3409  }
3410 
3411  Vector<bool> RowKept(m), ColKept(n);
3412  RowKept.Fill(false); ColKept.Fill(false);
3413  for (int i = 0; i < row.GetM(); i++)
3414  RowKept(row(i)) = true;
3415 
3416  for (int i = 0; i < col.GetM(); i++)
3417  ColKept(col(i)) = true;
3418 
3419  long nnz = 0;
3420  for (int i = 0; i < A.GetM(); i++)
3421  for (int j = 0; j < A.GetRowSize(i); j++)
3422  {
3423  if (ColKept(A.Index(i, j)) && RowKept(i))
3424  nnz++;
3425 
3426  if (A.Index(i, j) != i)
3427  if (RowKept(A.Index(i, j)) && ColKept(i))
3428  nnz++;
3429  }
3430 
3431  RowNum.Reallocate(nnz);
3432  ColNum.Reallocate(nnz);
3433  Value.Reallocate(nnz);
3434  nnz = 0;
3435  for (int i = 0; i < A.GetM(); i++)
3436  for (int j = 0; j < A.GetRowSize(i); j++)
3437  {
3438  if (ColKept(A.Index(i, j)) && RowKept(i))
3439  {
3440  RowNum(nnz) = i;
3441  ColNum(nnz) = A.Index(i, j);
3442  Value(nnz) = A.Value(i, j);
3443  nnz++;
3444  }
3445 
3446  if (A.Index(i, j) != i)
3447  if (RowKept(A.Index(i, j)) && ColKept(i))
3448  {
3449  RowNum(nnz) = A.Index(i, j);
3450  ColNum(nnz) = i;
3451  Value(nnz) = A.Value(i, j);
3452  nnz++;
3453  }
3454  }
3455  }
3456 
3457 
3459  template<class T0, class Prop0, class Allocator0,
3460  class T1, class Allocator1>
3461  void CopySubMatrix(const Matrix<T0, Prop0, RowSparse, Allocator0>& A,
3462  const IVect& row, const IVect& col,
3463  Vector<int>& RowNum,
3464  Vector<int>& ColNum,
3465  Vector<T1, VectFull, Allocator1>& Value)
3466  {
3467  int m = A.GetM(), n = A.GetN();
3468  if ((m <= 0) || (n <= 0) || (row.GetM() <= 0) || (col.GetM() <= 0))
3469  {
3470  RowNum.Clear(); ColNum.Clear(); Value.Clear();
3471  return;
3472  }
3473 
3474  long* ptr = A.GetPtr();
3475  int* ind = A.GetInd();
3476  T0* data = A.GetData();
3477  Vector<bool> RowKept(m), ColKept(n);
3478  RowKept.Fill(false); ColKept.Fill(false);
3479  for (int i = 0; i < row.GetM(); i++)
3480  RowKept(row(i)) = true;
3481 
3482  for (int i = 0; i < col.GetM(); i++)
3483  ColKept(col(i)) = true;
3484 
3485  long nnz = 0;
3486  for (int i = 0; i < m; i++)
3487  if (RowKept(i))
3488  for (long j = ptr[i]; j < ptr[i+1]; j++)
3489  if (ColKept(ind[j]))
3490  nnz++;
3491 
3492  RowNum.Reallocate(nnz);
3493  ColNum.Reallocate(nnz);
3494  Value.Reallocate(nnz);
3495  nnz = 0;
3496  for (int i = 0; i < m; i++)
3497  if (RowKept(i))
3498  for (long j = ptr[i]; j < ptr[i+1]; j++)
3499  if (ColKept(ind[j]))
3500  {
3501  RowNum(nnz) = i;
3502  ColNum(nnz) = ind[j];
3503  Value(nnz) = data[j];
3504  nnz++;
3505  }
3506  }
3507 
3508 
3510  template<class T0, class Prop0, class Allocator0,
3511  class T1, class Allocator1>
3512  void CopySubMatrix(const Matrix<T0, Prop0, RowSymSparse, Allocator0>& A,
3513  const IVect& row, const IVect& col,
3514  Vector<int>& RowNum,
3515  Vector<int>& ColNum,
3516  Vector<T1, VectFull, Allocator1>& Value)
3517  {
3518  int m = A.GetM(), n = A.GetN();
3519  if ((m <= 0) || (n <= 0) || (row.GetM() <= 0) || (col.GetM() <= 0))
3520  {
3521  RowNum.Clear(); ColNum.Clear(); Value.Clear();
3522  return;
3523  }
3524 
3525  long* ptr = A.GetPtr();
3526  int* ind = A.GetInd();
3527  T0* data = A.GetData();
3528  Vector<bool> RowKept(m), ColKept(n);
3529  RowKept.Fill(false); ColKept.Fill(false);
3530  for (int i = 0; i < row.GetM(); i++)
3531  RowKept(row(i)) = true;
3532 
3533  for (int i = 0; i < col.GetM(); i++)
3534  ColKept(col(i)) = true;
3535 
3536  long nnz = 0;
3537  for (int i = 0; i < m; i++)
3538  for (long j = ptr[i]; j < ptr[i+1]; j++)
3539  {
3540  if (ColKept(ind[j]) && RowKept(i))
3541  nnz++;
3542 
3543  if (ind[j] != i)
3544  if (RowKept(ind[j]) && ColKept(i))
3545  nnz++;
3546  }
3547 
3548  RowNum.Reallocate(nnz);
3549  ColNum.Reallocate(nnz);
3550  Value.Reallocate(nnz);
3551  nnz = 0;
3552  for (int i = 0; i < m; i++)
3553  for (long j = ptr[i]; j < ptr[i+1]; j++)
3554  {
3555  if (ColKept(ind[j]) && RowKept(i))
3556  {
3557  RowNum(nnz) = i;
3558  ColNum(nnz) = ind[j];
3559  Value(nnz) = data[j];
3560  nnz++;
3561  }
3562 
3563  if (ind[j] != i)
3564  if (RowKept(ind[j]) && ColKept(i))
3565  {
3566  RowNum(nnz) = ind[j];
3567  ColNum(nnz) = i;
3568  Value(nnz) = data[j];
3569  nnz++;
3570  }
3571  }
3572  }
3573 
3574 
3576  template<class T0, class Prop0, class Storage0, class Allocator0,
3577  class T1, class Prop1, class Storage1, class Allocator1>
3578  void CopySubMatrix(const Matrix<T0, Prop0, Storage0, Allocator0>& A,
3579  const IVect& row, const IVect& col,
3580  Matrix<T1, Prop1, Storage1, Allocator1>& B)
3581  {
3582  Vector<int> RowNum, ColNum;
3583  typedef typename Matrix<T0, Prop0, Storage0, Allocator0>::entry_type T;
3584  Vector<T> Value;
3585 
3586  // extracts rows/columns in coordinate format
3587  CopySubMatrix(A, row, col, RowNum, ColNum, Value);
3588 
3589  // converts to the sparse matrix B
3590  B.Reallocate(A.GetM(), A.GetN());
3591  ConvertMatrix_from_Coordinates(RowNum, ColNum, Value, B);
3592  }
3593 
3594 } // namespace Seldon
3595 
3596 #define SELDON_FILE_FUNCTIONS_MATRIX_ARRAY_CXX
3597 #endif
Seldon::Matrix< T, Prop, ArrayRowSymSparse, Allocator >::GetRowSize
int GetRowSize(int i) const
Returns the number of non-zero entries of a row.
Definition: Matrix_ArraySparseInline.cxx:1360
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::to_str
std::string to_str(const T &input)
Converts most types to string.
Definition: CommonInline.cxx:137
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::Vector
Definition: SeldonHeader.hxx:207
Seldon::WrongIndex
Definition: Errors.hxx:114
Seldon::Matrix
Definition: SeldonHeader.hxx:226
Seldon::Matrix< T, Prop, ArrayColSparse, Allocator >::GetColumnSize
int GetColumnSize(int i) const
Returns the number of non-zero entries of a column.
Definition: Matrix_ArraySparseInline.cxx:614
Seldon::ArrayRowSymSparse
Definition: Storage.hxx:132
Seldon::Matrix< T, Prop, RowSparse, Allocator >
Row-major sparse-matrix class.
Definition: Matrix_Sparse.hxx:223
Seldon::Conjugate
void Conjugate(Matrix< T, Prop, Storage, Allocator > &A)
A is replaced by its conjugate.
Definition: Functions_Matrix.cxx:2915
Seldon::VirtualMatrix::GetN
int GetN() const
Returns the number of columns.
Definition: Matrix_BaseInline.cxx:80
Seldon::VirtualMatrix::GetM
int GetM() const
Returns the number of rows.
Definition: Matrix_BaseInline.cxx:69
Seldon::Matrix_ArraySparse::Value
const T & Value(int num_row, int i) const
Returns j-th non-zero value of row/column i.
Definition: Matrix_ArraySparseInline.cxx:286
Seldon::MltMatrix
void MltMatrix(const Matrix< T0, Prop0, RowSparse, Allocator0 > &A, const Matrix< T1, Prop1, RowSparse, Allocator1 > &B, Matrix< T2, Prop2, RowSparse, Allocator2 > &C)
Multiplies two row-major sparse matrices in Harwell-Boeing format.
Definition: Functions_Matrix.cxx:207
Seldon::General
Definition: Properties.hxx:26
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::Matrix< T, Prop, ArrayRowSparse, Allocator >
Row-major sparse-matrix class.
Definition: Matrix_ArraySparse.hxx:214
Seldon::Matrix< T, Prop, RowSymSparse, Allocator >
Row-major sparse-matrix class.
Definition: Matrix_SymSparse.hxx:218
Seldon::Matrix_ArraySparse::Index
int Index(int num_row, int i) const
Returns column/row number of j-th non-zero value of row/column i.
Definition: Matrix_ArraySparseInline.cxx:325
Seldon::absSquare
T absSquare(const T &x)
returns the square modulus of z
Definition: CommonInline.cxx:340
Seldon::NormInf
ClassComplexType< T >::Treal NormInf(const Matrix< T, Prop, Storage, Allocator > &A)
Returns the infinity-norm of a matrix.
Definition: Functions_Matrix.cxx:2435
Seldon::Matrix< T, Prop, ArrayColSymSparse, Allocator >
Column-major symmetric sparse-matrix class.
Definition: Matrix_ArraySparse.hxx:254
Seldon::Transpose
void Transpose(Matrix< T, Prop, Storage, Allocator > &A)
Matrix transposition.
Definition: Functions_Matrix.cxx:2699
Seldon::RowSymSparse
Definition: Storage.hxx:114
Seldon::Symmetric
Definition: Properties.hxx:30
Seldon::MaxAbs
ClassComplexType< T >::Treal MaxAbs(const Matrix< T, Prop, Storage, Allocator > &A)
Returns the maximum (in absolute value) of a matrix.
Definition: Functions_Matrix.cxx:2386
Seldon::ArrayRowSparse
Definition: Storage.hxx:124
Seldon::RowSparse
Definition: Storage.hxx:91
Seldon
Seldon namespace.
Definition: Array.cxx:24
Seldon::Matrix< T, Prop, ArrayRowSymSparse, Allocator >
Row-major symmetric sparse-matrix class.
Definition: Matrix_ArraySparse.hxx:311
Seldon::Matrix< T, Prop, ArrayRowSparse, Allocator >::GetRowSize
int GetRowSize(int i) const
Returns the number of non-zero entries of a row.
Definition: Matrix_ArraySparseInline.cxx:751
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::Matrix< T, Prop, ArrayColSparse, Allocator >
Column-major sparse-matrix class.
Definition: Matrix_ArraySparse.hxx:172
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::Matrix< T, Prop, ArrayColSymSparse, Allocator >::GetColumnSize
int GetColumnSize(int i) const
Returns the number of non-zero entries of a column.
Definition: Matrix_ArraySparseInline.cxx:1054
Seldon::ComplexAbs
T ComplexAbs(const T &val)
returns absolute value of val
Definition: CommonInline.cxx:309