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