SuperLU.cxx
1 // Copyright (C) 2003-2009 Marc Duruflé
2 //
3 // This file is part of the linear-algebra library Seldon,
4 // http://seldon.sourceforge.net/.
5 //
6 // Seldon is free software; you can redistribute it and/or modify it under the
7 // terms of the GNU Lesser General Public License as published by the Free
8 // Software Foundation; either version 2.1 of the License, or (at your option)
9 // any later version.
10 //
11 // Seldon is distributed in the hope that it will be useful, but WITHOUT ANY
12 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
13 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
14 // more details.
15 //
16 // You should have received a copy of the GNU Lesser General Public License
17 // along with Seldon. If not, see http://www.gnu.org/licenses/.
18 
19 
20 #ifndef SELDON_FILE_SUPERLU_CXX
21 
22 #include "SuperLU.hxx"
23 
24 namespace Seldon
25 {
26 #ifndef SELDON_WITH_SUPERLU_DIST
27  void SetComplexOne(superlu::doublecomplex& one)
28  {
29  one.r = 1.0;
30  one.i = 0.0;
31  }
32 
33 
34  // The function comes from the Matlab interface to SuperLU. It is part of
35  // SuperLU package. Its copyright is held by University of California
36  // Berkeley, Xerox Palo Alto Research Center and Lawrence Berkeley National
37  // Lab. It is released under a license compatible with the GNU LGPL.
38  template<class T>
39  void LUextract(superlu::SuperMatrix *L, superlu::SuperMatrix *U,
40  T *Lval, int *Lrow, long *Lcol,
41  T *Uval, int *Urow, long *Ucol, int *snnzL,
42  int *snnzU)
43  {
44  int i, j, k;
45  int upper;
46  int fsupc, istart, nsupr;
47  long lastl = 0, lastu = 0;
48  T *SNptr;
49  T one;
50 
51  SetComplexOne(one);
52 
53 #ifdef SELDON_WITH_SUPERLU_MT
54  superlu::SCPformat *Lstore;
55  superlu::NCPformat *Ustore;
56  Lstore = static_cast<superlu::SCPformat*>(L->Store);
57  Ustore = static_cast<superlu::NCPformat*>(U->Store);
58 #else
59  superlu::SCformat *Lstore;
60  superlu::NCformat *Ustore;
61  Lstore = static_cast<superlu::SCformat*>(L->Store);
62  Ustore = static_cast<superlu::NCformat*>(U->Store);
63 #endif
64 
65  Lcol[0] = 0;
66  Ucol[0] = 0;
67 
68  /* for each supernode */
69  for (k = 0; k <= Lstore->nsuper; ++k) {
70 
71  fsupc = L_FST_SUPC(k);
72  istart = L_SUB_START(fsupc);
73  nsupr = L_SUB_START(fsupc+1) - istart;
74  upper = 1;
75 
76  /* for each column in the supernode */
77  for (j = fsupc; j < L_FST_SUPC(k+1); ++j) {
78  SNptr = &(static_cast<T*>(Lstore->nzval))[L_NZ_START(j)];
79 
80  /* Extract U */
81  for (i = U_NZ_START(j); i < U_NZ_START(j+1); ++i) {
82  Uval[lastu] = (static_cast<T*>(Ustore->nzval))[i];
83  Urow[lastu++] = U_SUB(i);
84  }
85  for (i = 0; i < upper; ++i) { /* upper triangle in the supernode */
86  Uval[lastu] = SNptr[i];
87  Urow[lastu++] = L_SUB(istart+i);
88  }
89  Ucol[j+1] = lastu;
90 
91  /* Extract L */
92  Lval[lastl] = one; /* unit diagonal */
93  Lrow[lastl++] = L_SUB(istart + upper - 1);
94  for (i = upper; i < nsupr; ++i) {
95  Lval[lastl] = SNptr[i];
96  Lrow[lastl++] = L_SUB(istart+i);
97  }
98  Lcol[j+1] = lastl;
99 
100  ++upper;
101 
102  } /* for j ... */
103 
104  } /* for k ... */
105 
106  *snnzL = lastl;
107  *snnzU = lastu;
108  }
109 #endif
110 
111 
112  /**********************
113  * MatrixSuperLU_Base *
114  **********************/
115 
116 
118  template<class T>
119  MatrixSuperLU_Base<T>::MatrixSuperLU_Base()
120  {
121  n = 0;
122 
123 #ifndef SELDON_WITH_SUPERLU_DIST
124 
125  permc_spec = superlu::COLAMD;
126  Lstore = NULL;
127  Ustore = NULL;
128 
129 #ifdef SELDON_WITH_SUPERLU_MT
130  nprocs = 1;
131  diag_pivot_thresh = 0.01;
132  usepr = superlu::NO;
133  drop_tol = 0.0;
134 
135 #else
136  superlu::set_default_options(&options);
137 #endif
138 
139 #else
140 
141  // distributed version
142  superlu::set_default_options_dist(&options);
143 
144  options.ParSymbFact = superlu::YES;
145  options.ColPerm = superlu::PARMETIS;
146  options.IterRefine = superlu::NOREFINE;
147 
148 #endif
149 
150  ShowMessages();
151  display_info = false;
152  info_facto = 0;
153  }
154 
155 
157  template<class T>
158  MatrixSuperLU_Base<T>::~MatrixSuperLU_Base()
159  {
160  Clear();
161  }
162 
163 
165  template<class T>
166  void MatrixSuperLU_Base<T>::Clear()
167  {
168  if (n > 0)
169  {
170 #ifndef SELDON_WITH_SUPERLU_DIST
171  // SuperLU objects are cleared
172  superlu::Destroy_SuperNode_Matrix(&L);
173  superlu::Destroy_CompCol_Matrix(&U);
174  if (permc_spec != superlu::MY_PERMC)
175  {
176  perm_r.Clear();
177  perm_c.Clear();
178  }
179 
180  superlu::StatFree(&stat);
181 
182 #else
183 
184  superlu::PStatFree(&stat);
185 
186  // permuted matrix A is cleared
187  superlu::Destroy_CompRowLoc_Matrix_dist(&A);
188 
189  superlu::superlu_gridexit(&grid);
190 #endif
191 
192  n = 0;
193  }
194  }
195 
196 
198  template<class T>
199  void MatrixSuperLU_Base<T>::Init(int size, int_t& panel_size, int_t& relax)
200  {
201  n = size;
202 
203 #ifndef SELDON_WITH_SUPERLU_DIST
204  panel_size = superlu::sp_ienv(1);
205  relax = superlu::sp_ienv(2);
206 
207 #ifdef SELDON_WITH_SUPERLU_MT
208  fact = superlu::EQUILIBRATE;
209  refact = superlu::NO;
210 
211  superlu::StatAlloc(n, nprocs, panel_size, relax, &stat);
212  superlu::StatInit(n, nprocs, &stat);
213 #else
214  superlu::StatInit(&stat);
215 #endif
216 
217 #else
218 
219  // Initialize the statistics variables.
220  superlu::PStatInit(&stat);
221 
222 #endif
223 
224  }
225 
226 
227  template<class T>
228  void MatrixSuperLU_Base<T>::SetNumberOfThreadPerNode(int p)
229  {
230 #ifdef SELDON_WITH_SUPERLU_MT
231  nprocs = p;
232 #endif
233  }
234 
235 
237  template<class T>
238  void MatrixSuperLU_Base<T>::HideMessages()
239  {
240  display_info = false;
241  }
242 
243 
245  template<class T>
246  void MatrixSuperLU_Base<T>::ShowMessages()
247  {
248  display_info = true;
249  }
250 
251 
252  template<class T>
253  bool MatrixSuperLU_Base<T>::UseInteger8() const
254  {
255  if (sizeof(int_t) == 8)
256  return true;
257 
258  return false;
259  }
260 
261 
263  template<class T>
264  int MatrixSuperLU_Base<T>::GetInfoFactorization() const
265  {
266  return info_facto;
267  }
268 
269 
271 
277  template<class T>
278  const Vector<int_t>& MatrixSuperLU_Base<T>::GetRowPermutation() const
279  {
280  return perm_r;
281  }
282 
283 
285 
291  template<class T>
292  const Vector<int_t>& MatrixSuperLU_Base<T>::GetColPermutation() const
293  {
294  return perm_c;
295  }
296 
297 
298  template<class T>
299  void MatrixSuperLU_Base<T>::SelectOrdering(int type)
300  {
301  permc_spec = (colperm_t) type;
302  }
303 
304 
305  template<class T>
306  void MatrixSuperLU_Base<T>::SetPermutation(const IVect& permut)
307  {
308  permc_spec = superlu::MY_PERMC;
309  perm_c.Reallocate(permut.GetM());
310  for (int i = 0; i < permut.GetM(); i++)
311  perm_c(i) = permut(i);
312 
313  perm_r.Reallocate(perm_c.GetM());
314  perm_r.Fill();
315  }
316 
317 
318  /*************************
319  * MatrixSuperLU<double> *
320  *************************/
321 
322  void MatrixSuperLU<double>::Clear()
323  {
324 #ifdef SELDON_WITH_SUPERLU_DIST
325  superlu::dScalePermstructFree(&ScalePermstruct);
326  superlu::dDestroy_LU(n, &grid, &LUstruct);
327  superlu::dLUstructFree(&LUstruct);
328 #endif
329 
330  MatrixSuperLU_Base<double>::Clear();
331  }
332 
333 
334 #ifndef SELDON_WITH_SUPERLU_DIST
335 
345  template<class Prop, class Allocator>
346  void MatrixSuperLU<double>
347  ::GetLU(Matrix<double, Prop, ColSparse, Allocator>& Lmat,
348  Matrix<double, Prop, ColSparse, Allocator>& Umat,
349  bool permuted)
350  {
351 #ifdef SELDON_WITH_SUPERLU_MT
352  Lstore = static_cast<superlu::SCPformat*>(L.Store);
353  Ustore = static_cast<superlu::NCPformat*>(U.Store);
354 #else
355  Lstore = static_cast<superlu::SCformat*>(L.Store);
356  Ustore = static_cast<superlu::NCformat*>(U.Store);
357 #endif
358 
359  int Lnnz = Lstore->nnz;
360  int Unnz = Ustore->nnz;
361 
362  int m = U.nrow;
363  int n = U.ncol;
364 
365  Vector<double, VectFull, Allocator> Lval(Lnnz);
366  Vector<int> Lrow(Lnnz);
367  Vector<long> Lcol(n + 1);
368 
369  Vector<double, VectFull, Allocator> Uval(Unnz);
370  Vector<int> Urow(Unnz);
371  Vector<long> Ucol(n + 1);
372 
373  int Lsnnz;
374  int Usnnz;
375  LUextract(&L, &U, Lval.GetData(), Lrow.GetData(), Lcol.GetData(),
376  Uval.GetData(), Urow.GetData(), Ucol.GetData(), &Lsnnz, &Usnnz);
377 
378  Lmat.SetData(m, n, Lval, Lcol, Lrow);
379  Umat.SetData(m, n, Uval, Ucol, Urow);
380 
381  if (!permuted)
382  {
383  Vector<int> row_perm_orig = perm_r;
384  Vector<int> col_perm_orig = perm_c;
385 
386  Vector<int> row_perm(n);
387  Vector<int> col_perm(n);
388  row_perm.Fill();
389  col_perm.Fill();
390 
391  Sort(row_perm_orig, row_perm);
392  Sort(col_perm_orig, col_perm);
393 
394  ApplyInversePermutation(Lmat, row_perm, col_perm);
395  ApplyInversePermutation(Umat, row_perm, col_perm);
396  }
397  }
398 
399 
401 
412  template<class Prop, class Allocator>
413  void MatrixSuperLU<double>
414  ::GetLU(Matrix<double, Prop, RowSparse, Allocator>& Lmat,
415  Matrix<double, Prop, RowSparse, Allocator>& Umat,
416  bool permuted)
417  {
418  Lmat.Clear();
419  Umat.Clear();
420 
421  Matrix<double, Prop, ColSparse, Allocator> Lmat_col;
422  Matrix<double, Prop, ColSparse, Allocator> Umat_col;
423  GetLU(Lmat_col, Umat_col, permuted);
424 
425  Copy(Lmat_col, Lmat);
426  Lmat_col.Clear();
427  Copy(Umat_col, Umat);
428  Umat_col.Clear();
429  }
430 
431 
433  size_t MatrixSuperLU<double>::GetMemorySize() const
434  {
435  size_t taille = sizeof(int)*(perm_r.GetM()+perm_c.GetM());
436  if (this->n > 0)
437  {
438 #ifdef SELDON_WITH_SUPERLU_MT
439  superlu::superlu_memusage_t mem_usage;
440  int_t panel_size = superlu::sp_ienv(1);
441  superlu::superlu_dQuerySpace(nprocs, const_cast<superlu::SuperMatrix*>(&L),
442  const_cast<superlu::SuperMatrix*>(&U),
443  panel_size, &mem_usage);
444 #else
445  superlu::mem_usage_t mem_usage;
446  superlu::dQuerySpace(const_cast<superlu::SuperMatrix*>(&L),
447  const_cast<superlu::SuperMatrix*>(&U), &mem_usage);
448 #endif
449 
450  taille += mem_usage.total_needed;
451  }
452 
453  return taille;
454  }
455 
456 
458  template<class T0, class Prop, class Storage, class Allocator>
459  void MatrixSuperLU<double>::
460  FactorizeMatrix(Matrix<T0, Prop, Storage, Allocator> & mat,
461  bool keep_matrix)
462  {
463  // clearing previous factorization
464  Clear();
465 
466  // conversion in CSC format
467  General prop;
468  Vector<superlu_int_t> Ptr, IndRow;
469  Vector<double> Val;
470 
471  ConvertToCSC(mat, prop, Ptr, IndRow, Val, false);
472  if (!keep_matrix)
473  mat.Clear();
474 
475  FactorizeCSC(Ptr, IndRow, Val, false);
476  }
477 
478 
480  void MatrixSuperLU<double>
481  ::FactorizeCSC(Vector<superlu_int_t>& Ptr, Vector<superlu_int_t>& IndRow,
482  Vector<double>& Val, bool sym)
483  {
484  // initializing parameters
485  int_t panel_size, relax;
486  int_t lwork = 0;
487  Init(Ptr.GetM()-1, panel_size, relax);
488 
489  superlu::SuperMatrix AA;
490  int_t nnz = IndRow.GetM();
491  superlu::dCreate_CompCol_Matrix(&AA, n, n, nnz, Val.GetData(),
492  reinterpret_cast<int_t*>(IndRow.GetData()),
493  reinterpret_cast<int_t*>(Ptr.GetData()),
494  superlu::SLU_NC,
495  superlu::SLU_D, superlu::SLU_GE);
496 
497  // we get renumbering vectors perm_r and perm_c
498  options.ColPerm = permc_spec;
499  if (permc_spec != superlu::MY_PERMC)
500  {
501  perm_r.Reallocate(n);
502  perm_c.Reallocate(n);
503  perm_r.Fill();
504  perm_c.Fill();
505 
506  superlu::get_perm_c(permc_spec, &AA, perm_c.GetData());
507  }
508 
509 #ifdef SELDON_WITH_SUPERLU_MT
510  superlu::SuperMatrix AC;
511  superlu::trans_t trans = superlu::NOTRANS;
512 
513  superlu::pdgstrf_init(nprocs, fact, trans, refact, panel_size, relax,
514  diag_pivot_thresh, usepr, drop_tol,
515  perm_c.GetData(), perm_r.GetData(),
516  NULL, lwork, &AA, &AC, &options, &stat);
517 
518  int_t info;
519  superlu::pdgstrf(&options, &AC, perm_r.GetData(), &L, &U, &stat, &info);
520 
521  info_facto = info;
522  superlu::pxgstrf_finalize(&options, &AC);
523 
524  if (info_facto == 0 && display_info)
525  {
526  superlu::PrintStat(&stat);
527  }
528 
529 #else
530  superlu::SuperMatrix A;
531  // original matrix AA is permuted to obtain matrix A
532  Vector<int> etree(n);
533  sp_preorder(&options, &AA, perm_c.GetData(), etree.GetData(), &A);
534 
535  // then calling factorisation on permuted matrix
536  superlu::dgstrf(&options, &A, relax, panel_size, etree.GetData(),
537  NULL, lwork, perm_c.GetData(), perm_r.GetData(), &L, &U,
538  &Glu, &stat, &info_facto);
539 
540  if (info_facto == 0 && display_info)
541  {
542  superlu::mem_usage_t mem_usage;
543  Lstore = (superlu::SCformat *) L.Store;
544  Ustore = (superlu::NCformat *) U.Store;
545  cout << "Number of nonzeros in factor L = " << Lstore->nnz << endl;
546  cout << "Number of nonzeros in factor U = " << Ustore->nnz << endl;
547  cout << "Number of nonzeros in L+U = "
548  << Lstore->nnz + Ustore->nnz << endl;
549  superlu::dQuerySpace(&L, &U, &mem_usage);
550  cout << "Memory used for factorization in MB: "
551  << mem_usage.total_needed / (1024. * 1024.) << endl;
552  }
553 
554  superlu::Destroy_CompCol_Permuted(&A);
555 #endif
556 
557  // clearing matrices
558  superlu::Destroy_CompCol_Matrix(&AA);
559 
560  Ptr.Nullify(); IndRow.Nullify(); Val.Nullify();
561  }
562 
563 
565  template<class Allocator2>
566  void MatrixSuperLU<double>::Solve(Vector<double, VectFull, Allocator2>& x)
567  {
568  Solve(SeldonNoTrans, x);
569  }
570 
571 
573  template<class Allocator2>
574  void MatrixSuperLU<double>::Solve(const SeldonTranspose& TransA,
575  Vector<double, VectFull, Allocator2>& x)
576  {
577  Solve(TransA, x.GetData(), 1);
578  }
579 
580 
582  void MatrixSuperLU<double>::Solve(const SeldonTranspose& TransA,
583  double* x_ptr, int nrhs_)
584  {
585  superlu::trans_t trans;
586  if (TransA.NoTrans())
587  trans = superlu::NOTRANS;
588  else
589  trans = superlu::TRANS;
590 
591  int_t nb_rhs = nrhs_, info;
592  // Putting right hand side on SuperLU structure.
593  superlu::dCreate_Dense_Matrix(&B, n, nb_rhs, x_ptr, n,
594  superlu::SLU_DN, superlu::SLU_D, superlu::SLU_GE);
595 
596 #ifdef SELDON_WITH_SUPERLU_MT
597  superlu::dgstrs(trans, &L, &U, perm_r.GetData(),
598  perm_c.GetData(), &B, &stat, &info);
599 #else
600  superlu::dgstrs(trans, &L, &U, perm_c.GetData(),
601  perm_r.GetData(), &B, &stat, &info);
602 #endif
603 
604  superlu::Destroy_SuperMatrix_Store(&B);
605  }
606 
607 
609  template<class Allocator2>
610  void MatrixSuperLU<double>::Solve(Matrix<double, General, ColMajor, Allocator2>& x)
611  {
612  Solve(SeldonNoTrans, x);
613  }
614 
615 
617  template<class Allocator2>
618  void MatrixSuperLU<double>::Solve(const SeldonTranspose& TransA,
619  Matrix<double, General, ColMajor, Allocator2>& x)
620  {
621  Solve(TransA, x.GetData(), x.GetN());
622  }
623 
624 
625 #else
626 
627  /**********************************
628  * Distributed version for double *
629  **********************************/
630 
631 
633  size_t MatrixSuperLU<double>::GetMemorySize() const
634  {
635  size_t size = 0;
636  if (n > 0)
637  {
638  superlu::superlu_dist_mem_usage_t mem_usage;
639  superlu::dQuerySpace_dist(n, const_cast<superlu::dLUstruct_t*>(&LUstruct),
640  const_cast<superlu::gridinfo_t*>(&grid),
641  const_cast<superlu::SuperLUStat_t*>(&stat), &mem_usage);
642 
643  size += mem_usage.total;
644  }
645 
646  return size;
647  }
648 
649 
651  template<class T0, class Prop, class Storage, class Allocator>
652  void MatrixSuperLU<double>::
653  FactorizeMatrix(Matrix<T0, Prop, Storage, Allocator> & mat,
654  bool keep_matrix)
655  {
656  Vector<long> Ptr; Vector<int> IndRow;
657  Vector<double> Val; General prop;
658  ConvertToCSC(mat, prop, Ptr, IndRow, Val, false);
659  if (!keep_matrix)
660  mat.Clear();
661 
662  FactorizeCSC(Ptr, IndRow, Val, IsSymmetricMatrix(mat));
663  }
664 
665 
666  void MatrixSuperLU<double>
667  ::FactorizeCSC(Vector<long>& Ptr, Vector<int>& IndRow,
668  Vector<double>& Val, bool sym)
669  {
670  Vector<int> glob_num(1);
671  glob_num(0) = 0;
672  MPI_Comm comm = MPI_COMM_SELF;
673  FactorizeDistributedMatrix(comm, Ptr, IndRow, Val,
674  glob_num, sym, false);
675  }
676 
677 
679  template<class Allocator2>
680  void MatrixSuperLU<double>::Solve(Vector<double, VectFull, Allocator2>& x)
681  {
682  Solve(SeldonNoTrans, x);
683  }
684 
685 
687  template<class Allocator2>
688  void MatrixSuperLU<double>::Solve(const SeldonTranspose& TransA,
689  Vector<double, VectFull, Allocator2>& x)
690  {
691  Vector<int> glob_num(1);
692  glob_num(0) = 0;
693  MPI_Comm comm = MPI_COMM_SELF;
694  SolveDistributed(comm, TransA, x, glob_num);
695  }
696 
697 
698  void MatrixSuperLU<double>::Solve(const SeldonTranspose& TransA,
699  double* x_ptr, int nrhs_)
700  {
701  Vector<int> glob_num(1);
702  glob_num(0) = 0;
703  MPI_Comm comm = MPI_COMM_SELF;
704  SolveDistributed(comm, TransA, x_ptr, nrhs_, glob_num);
705  }
706 
707 
709  template<class Allocator2>
710  void MatrixSuperLU<double>::Solve(Matrix<double, General, ColMajor, Allocator2>& x)
711  {
712  Solve(SeldonNoTrans, x);
713  }
714 
715 
717  template<class Allocator2>
718  void MatrixSuperLU<double>::Solve(const SeldonTranspose& TransA,
719  Matrix<double, General, ColMajor, Allocator2>& x)
720  {
721  Vector<int> glob_num(1);
722  glob_num(0) = 0;
723  MPI_Comm comm = MPI_COMM_SELF;
724  SolveDistributed(comm, TransA, x, glob_num);
725  }
726 
727 
728  void MatrixSuperLU<double>::
729  FactorizeDistributedMatrix(MPI_Comm& comm_facto, Vector<long>& Ptr,
730  Vector<int>& IndRow, Vector<double>& Val,
731  const Vector<int>& glob_number,
732  bool sym, bool keep_matrix)
733  {
734  Vector<int> PtrInt(Ptr.GetM());
735  for (int i = 0; i < Ptr.GetM(); i++)
736  PtrInt(i) = Ptr(i);
737 
738  FactorizeParallel(comm_facto, PtrInt, IndRow, Val,
739  glob_number, sym, keep_matrix);
740  }
741 
742 
743  void MatrixSuperLU<double>::
744  FactorizeDistributedMatrix(MPI_Comm& comm_facto, Vector<int64_t>& Ptr,
745  Vector<int64_t>& IndRow, Vector<double>& Val,
746  const Vector<int>& glob_number,
747  bool sym, bool keep_matrix)
748  {
749  FactorizeParallel(comm_facto, Ptr, IndRow, Val,
750  glob_number, sym, keep_matrix);
751  }
752 
753 
754  template<class Tint>
755  void MatrixSuperLU<double>::
756  FactorizeParallel(MPI_Comm& comm_facto,
757  Vector<Tint>&, Vector<Tint>&,
758  Vector<double>&,
759  const Vector<int>& glob_num,
760  bool sym, bool keep_matrix)
761  {
762  cout << "Not available for this type of integers " << endl;
763  cout << "Size of Tint = " << sizeof(Tint) << endl;
764  abort();
765  }
766 
767 
768  void MatrixSuperLU<double>::
769  FactorizeParallel(MPI_Comm& comm_facto,
770  Vector<superlu_int_t>& Ptr, Vector<superlu_int_t>& Row,
771  Vector<double>& Val, const Vector<int>& glob_num,
772  bool sym, bool keep_matrix)
773  {
774 #ifdef SELDON_WITH_SUPERLU_DOUBLE
775  // previous factorization is cleared if present
776  Clear();
777 
778  if (!sym)
779  {
780  cout << "Problem with option TRANS of SuperLU" << endl;
781  abort();
782  }
783 
784  // m_loc : local number of rows
785  // N : global number of rows
786  int m_loc_ = Ptr.GetM()-1;
787  int N = m_loc_; nloc = m_loc_;
788  MPI_Allreduce(&m_loc_, &N, 1, MPI_INTEGER, MPI_SUM, comm_facto);
789  int_t m_loc = m_loc_;
790 
791  // structures are initialized with N
792  int_t panel_size, relax;
793  Init(N, panel_size, relax);
794  superlu::dScalePermstructInit(N, N, &ScalePermstruct);
795  superlu::dLUstructInit(n, &LUstruct);
796 
797  // 1-D grid
798  MPI_Comm_size(comm_facto, &nprow);
799  npcol = 1;
800 
801  // initialize the superlu process grid
802  superlu::superlu_gridinit(comm_facto, nprow, npcol, &grid);
803 
804  // fills the superlu structure
805  // global numbers are assumed to be consecutive, we provide the first row number
806  int_t fst_row = glob_num(0);
807  int_t nnz_loc = Row.GetM();
808  //superlu::SuperMatrix A;
809  superlu::
810  dCreate_CompRowLoc_Matrix_dist(&A, n, n, nnz_loc, m_loc, fst_row,
811  Val.GetData(), Row.GetData(), Ptr.GetData(),
812  superlu::SLU_NR_loc, superlu::SLU_D,
813  superlu::SLU_GE);
814 
815  // completes factorization
816  int_t nrhs = 0;
817  options.Trans = superlu::NOTRANS;
818  options.Fact = superlu::DOFACT;
819  superlu::
820  pdgssvx(&options, &A, &ScalePermstruct,
821  NULL, m_loc, nrhs, &grid,
822  &LUstruct, &SOLVEstruct, NULL, &stat, &info_facto);
823 
824  Ptr.Nullify(); Row.Nullify(); Val.Nullify();
825 #else
826  cout << "Recompile Seldon with SELDON_WITH_SUPERLU_DOUBLE" << endl;
827  cout << "Double and complex<double> cannot be handled simultaneously" << endl;
828  abort();
829 #endif
830  }
831 
832 
833  template<class Allocator2>
834  void MatrixSuperLU<double>::
835  SolveDistributed(MPI_Comm& comm_facto,
836  const SeldonTranspose& TransA,
837  Vector<double, VectFull, Allocator2>& x,
838  const IVect& glob_num)
839  {
840  SolveDistributed(comm_facto, TransA, x.GetData(), 1, glob_num);
841  }
842 
843 
844  template<class Allocator2>
845  void MatrixSuperLU<double>::
846  SolveDistributed(MPI_Comm& comm_facto,
847  const SeldonTranspose& TransA,
848  Matrix<double, General, ColMajor, Allocator2>& x,
849  const IVect& glob_num)
850  {
851  SolveDistributed(comm_facto, TransA, x.GetData(), x.GetN(), glob_num);
852  }
853 
854 
855  void MatrixSuperLU<double>::SolveDistributed(MPI_Comm& comm_facto,
856  const SeldonTranspose& TransA,
857  double* x_ptr, int nrhs_,
858  const IVect& glob_num)
859  {
860 #ifdef SELDON_WITH_SUPERLU_DOUBLE
861  options.Fact = superlu::FACTORED;
862  // inverting transpose because we have provided columns instead of rows
863  if (TransA.NoTrans())
864  options.Trans = superlu::TRANS;
865  else
866  options.Trans = superlu::NOTRANS;
867 
868  options.Trans = superlu::NOTRANS;
869  Vector<double> berr(nrhs_);
870  for (int k = 0; k < nrhs_; k++)
871  {
872  // workaround, loop over right hand sides
873  int_t nrhs = 1; int_t info;
874  // int_t nrhs = nrhs_, info;
875  superlu::
876  pdgssvx(&options, &A, &ScalePermstruct, x_ptr,
877  n, nrhs, &grid, &LUstruct, &SOLVEstruct,
878  berr.GetData(), &stat, &info);
879 
880  x_ptr += nloc;
881  }
882 #endif
883  }
884 #endif
885 
886 
887  /********************************
888  * MatrixSuperLU<complexdouble> *
889  ********************************/
890 
891  void MatrixSuperLU<complex<double> >::Clear()
892  {
893 #ifdef SELDON_WITH_SUPERLU_DIST
894  superlu::zScalePermstructFree(&ScalePermstruct);
895  superlu::zDestroy_LU(n, &grid, &LUstruct);
896  superlu::zLUstructFree(&LUstruct);
897 #endif
898 
899  MatrixSuperLU_Base<complex<double> >::Clear();
900  }
901 
902 
903 #ifndef SELDON_WITH_SUPERLU_DIST
904 
914  template<class Prop, class Allocator>
915  void MatrixSuperLU<complex<double> >
916  ::GetLU(Matrix<complex<double>, Prop, ColSparse, Allocator>& Lmat,
917  Matrix<complex<double>, Prop, ColSparse, Allocator>& Umat,
918  bool permuted)
919  {
920 #ifdef SELDON_WITH_SUPERLU_MT
921  Lstore = static_cast<superlu::SCPformat*>(L.Store);
922  Ustore = static_cast<superlu::NCPformat*>(U.Store);
923 #else
924  Lstore = static_cast<superlu::SCformat*>(L.Store);
925  Ustore = static_cast<superlu::NCformat*>(U.Store);
926 #endif
927 
928  int Lnnz = Lstore->nnz;
929  int Unnz = Ustore->nnz;
930 
931  int m = U.nrow;
932  int n = U.ncol;
933 
934  Vector<complex<double>, VectFull, Allocator> Lval(Lnnz);
935  Vector<int> Lrow(Lnnz);
936  Vector<long> Lcol(n + 1);
937 
938  Vector<complex<double>, VectFull, Allocator> Uval(Unnz);
939  Vector<int> Urow(Unnz);
940  Vector<long> Ucol(n + 1);
941 
942  int Lsnnz;
943  int Usnnz;
944  LUextract(&L, &U, reinterpret_cast<superlu::doublecomplex*>(Lval.GetData()),
945  Lrow.GetData(), Lcol.GetData(),
946  reinterpret_cast<superlu::doublecomplex*>(Uval.GetData()),
947  Urow.GetData(), Ucol.GetData(), &Lsnnz, &Usnnz);
948 
949  Lmat.SetData(m, n, Lval, Lcol, Lrow);
950  Umat.SetData(m, n, Uval, Ucol, Urow);
951 
952  if (!permuted)
953  {
954  Vector<int> row_perm_orig = perm_r;
955  Vector<int> col_perm_orig = perm_c;
956 
957  Vector<int> row_perm(n);
958  Vector<int> col_perm(n);
959  row_perm.Fill();
960  col_perm.Fill();
961 
962  Sort(row_perm_orig, row_perm);
963  Sort(col_perm_orig, col_perm);
964 
965  ApplyInversePermutation(Lmat, row_perm, col_perm);
966  ApplyInversePermutation(Umat, row_perm, col_perm);
967  }
968 
969  }
970 
971 
973 
984  template<class Prop, class Allocator>
985  void MatrixSuperLU<complex<double> >
986  ::GetLU(Matrix<complex<double>, Prop, RowSparse, Allocator>& Lmat,
987  Matrix<complex<double>, Prop, RowSparse, Allocator>& Umat,
988  bool permuted)
989  {
990  Lmat.Clear();
991  Umat.Clear();
992 
993  Matrix<complex<double>, Prop, ColSparse, Allocator> Lmat_col;
994  Matrix<complex<double>, Prop, ColSparse, Allocator> Umat_col;
995  GetLU(Lmat_col, Umat_col, permuted);
996 
997  Copy(Lmat_col, Lmat);
998  Lmat_col.Clear();
999  Copy(Umat_col, Umat);
1000  Umat_col.Clear();
1001  }
1002 
1003 
1005  size_t MatrixSuperLU<complex<double> >::GetMemorySize() const
1006  {
1007  size_t taille = sizeof(int)*(perm_r.GetM()+perm_c.GetM());
1008  if (this->n > 0)
1009  {
1010 #ifdef SELDON_WITH_SUPERLU_MT
1011  superlu::superlu_memusage_t mem_usage;
1012  int_t panel_size = superlu::sp_ienv(1);
1013  superlu::superlu_zQuerySpace(nprocs, const_cast<superlu::SuperMatrix*>(&L),
1014  const_cast<superlu::SuperMatrix*>(&U),
1015  panel_size, &mem_usage);
1016 #endif
1017 #ifdef SELDON_WITH_SUPERLU_DIST
1018  superlu::superlu_dist_mem_usage_t mem_usage;
1019  superlu::zQuerySpace(const_cast<superlu::SuperMatrix*>(&L),
1020  const_cast<superlu::SuperMatrix*>(&U), &mem_usage);
1021 #endif
1022 
1023 #ifdef SELDON_WITH_SUPERLU_DIST
1024  superlu::superlu_dist_mem_usage_t mem_usage;
1025  superlu::zQuerySpace(const_cast<superlu::SuperMatrix*>(&L),
1026  const_cast<superlu::SuperMatrix*>(&U), &mem_usage);
1027 #endif
1028 
1029 #ifdef SELDON_WITH_SUPERLU_SEQ
1030  superlu::mem_usage_t mem_usage;
1031  superlu::zQuerySpace(const_cast<superlu::SuperMatrix*>(&L),
1032  const_cast<superlu::SuperMatrix*>(&U), &mem_usage);
1033 #endif
1034 
1035  taille += mem_usage.total_needed;
1036  }
1037 
1038  return taille;
1039  }
1040 
1041 
1043  template<class T0, class Prop, class Storage, class Allocator>
1044  void MatrixSuperLU<complex<double> >::
1045  FactorizeMatrix(Matrix<T0, Prop, Storage, Allocator> & mat,
1046  bool keep_matrix)
1047  {
1048  // clearing previous factorization
1049  Clear();
1050 
1051  // conversion in CSC format
1052  General prop;
1053  Vector<superlu_int_t> Ptr; Vector<superlu_int_t> IndRow;
1054  Vector<complex<double> > Val;
1055 
1056  ConvertToCSC(mat, prop, Ptr, IndRow, Val, false);
1057  if (!keep_matrix)
1058  mat.Clear();
1059 
1060  FactorizeCSC(Ptr, IndRow, Val, false);
1061  }
1062 
1063 
1064  void MatrixSuperLU<complex<double> >
1065  ::FactorizeCSC(Vector<superlu_int_t>& Ptr, Vector<superlu_int_t>& IndRow,
1066  Vector<complex<double> >& Val, bool sym)
1067  {
1068  // initializing parameters
1069  int_t panel_size, relax;
1070  int_t lwork = 0;
1071  Init(Ptr.GetM()-1, panel_size, relax);
1072 
1073  superlu::SuperMatrix AA;
1074  int_t nnz = IndRow.GetM();
1075  superlu::
1076  zCreate_CompCol_Matrix(&AA, n, n, nnz,
1077  reinterpret_cast<superlu::doublecomplex*>(Val.GetData()),
1078  reinterpret_cast<int_t*>(IndRow.GetData()),
1079  reinterpret_cast<int_t*>(Ptr.GetData()),
1080  superlu::SLU_NC, superlu::SLU_Z, superlu::SLU_GE);
1081 
1082  // We get renumbering vectors perm_r and perm_c.
1083  options.ColPerm = permc_spec;
1084  if (permc_spec != superlu::MY_PERMC)
1085  {
1086  perm_r.Reallocate(n);
1087  perm_c.Reallocate(n);
1088  perm_r.Fill();
1089  perm_c.Fill();
1090 
1091  superlu::get_perm_c(permc_spec, &AA, perm_c.GetData());
1092  }
1093 
1094 #ifdef SELDON_WITH_SUPERLU_MT
1095  superlu::SuperMatrix AC;
1096  superlu::trans_t trans = superlu::NOTRANS;
1097  perm_r.Reallocate(n);
1098  perm_c.Reallocate(n);
1099 
1100  superlu::
1101  pzgstrf_init(nprocs, fact, trans, refact, panel_size, relax,
1102  diag_pivot_thresh, usepr, drop_tol,
1103  perm_c.GetData(), perm_r.GetData(),
1104  NULL, lwork, &AA, &AC, &options, &stat);
1105 
1106  int_t info;
1107  superlu::
1108  pzgstrf(&options, &AC, perm_r.GetData(), &L, &U, &stat, &info);
1109 
1110  info_facto = info;
1111  superlu::pxgstrf_finalize(&options, &AC);
1112 
1113  if (info_facto == 0 && display_info)
1114  {
1115  cout << "Memory used for factorization in MiB: "
1116  << this->GetMemorySize() / (1024. * 1024.) << endl;
1117  }
1118 
1119 #else
1120  superlu::SuperMatrix A;
1121  // permuting matrix
1122  Vector<int> etree(n);
1123  superlu::sp_preorder(&options, &AA, perm_c.GetData(), etree.GetData(), &A);
1124 
1125  // factorisation
1126  superlu::zgstrf(&options, &A, relax, panel_size, etree.GetData(),
1127  NULL, lwork, perm_c.GetData(), perm_r.GetData(), &L, &U,
1128  &Glu, &stat, &info_facto);
1129 
1130  if (info_facto == 0 && display_info)
1131  {
1132  superlu::mem_usage_t mem_usage;
1133  Lstore = (superlu::SCformat *) L.Store;
1134  Ustore = (superlu::NCformat *) U.Store;
1135  cout << "Number of nonzeros in factor L = " << Lstore->nnz<<endl;
1136  cout << "Number of nonzeros in factor U = " << Ustore->nnz<<endl;
1137  cout << "Number of nonzeros in L+U = "
1138  << Lstore->nnz + Ustore->nnz<<endl;
1139  superlu::zQuerySpace(&L, &U, &mem_usage);
1140  cout << "Memory used for factorization in MiB: "
1141  << mem_usage.total_needed / (1024. * 1024.) << endl;
1142  }
1143 
1144  superlu::Destroy_CompCol_Permuted(&A);
1145 
1146 #endif
1147 
1148  // clearing matrices
1149  superlu::Destroy_CompCol_Matrix(&AA);
1150 
1151  Ptr.Nullify(); IndRow.Nullify(); Val.Nullify();
1152  }
1153 
1154 
1156  template<class Allocator2>
1157  void MatrixSuperLU<complex<double> >::
1158  Solve(Vector<complex<double>, VectFull, Allocator2>& x)
1159  {
1160  Solve(SeldonNoTrans, x);
1161  }
1162 
1163 
1165  template<class Allocator2>
1166  void MatrixSuperLU<complex<double> >::
1167  Solve(const SeldonTranspose& TransA,
1168  Vector<complex<double>, VectFull, Allocator2>& x)
1169  {
1170  Solve(TransA, x.GetData(), 1);
1171  }
1172 
1173 
1175  void MatrixSuperLU<complex<double> >
1176  ::Solve(const SeldonTranspose& TransA, complex<double>* x_ptr, int nrhs_)
1177  {
1178  superlu::trans_t trans = superlu::NOTRANS;
1179  if (TransA.Trans())
1180  trans = superlu::TRANS;
1181 
1182  int_t nb_rhs = nrhs_, info;
1183  superlu::
1184  zCreate_Dense_Matrix(&B, n, nb_rhs,
1185  reinterpret_cast<superlu::doublecomplex*>(x_ptr),
1186  n, superlu::SLU_DN, superlu::SLU_Z, superlu::SLU_GE);
1187 
1188 #ifdef SELDON_WITH_SUPERLU_MT
1189  superlu::zgstrs(trans, &L, &U, perm_r.GetData(),
1190  perm_c.GetData(), &B, &stat, &info);
1191 #else
1192  superlu::zgstrs(trans, &L, &U, perm_c.GetData(),
1193  perm_r.GetData(), &B, &stat, &info);
1194 #endif
1195 
1196  superlu::Destroy_SuperMatrix_Store(&B);
1197  }
1198 
1199 
1201  template<class Allocator2>
1202  void MatrixSuperLU<complex<double> >::
1203  Solve(Matrix<complex<double>, General, ColMajor, Allocator2>& x)
1204  {
1205  Solve(SeldonNoTrans, x);
1206  }
1207 
1208 
1210  template<class Allocator2>
1211  void MatrixSuperLU<complex<double> >::
1212  Solve(const SeldonTranspose& TransA,
1213  Matrix<complex<double>, General, ColMajor, Allocator2>& x)
1214  {
1215  Solve(TransA, x.GetData(), x.GetN());
1216  }
1217 
1218 #else
1219 
1220 
1221  /*****************************************
1222  * Distributed version for complexdouble *
1223  *****************************************/
1224 
1225 
1227  size_t MatrixSuperLU<complex<double> >::GetMemorySize() const
1228  {
1229  size_t size = 0;
1230 #ifndef SELDON_WITH_SUPERLU_DOUBLE
1231  if (n > 0)
1232  {
1233  superlu::superlu_dist_mem_usage_t mem_usage;
1234  superlu::zQuerySpace_dist(n, const_cast<superlu::zLUstruct_t*>(&LUstruct),
1235  const_cast<superlu::gridinfo_t*>(&grid),
1236  const_cast<superlu::SuperLUStat_t*>(&stat), &mem_usage);
1237 
1238  size += mem_usage.total;
1239  }
1240 #endif
1241 
1242  return size;
1243  }
1244 
1245 
1247  template<class T0, class Prop, class Storage, class Allocator>
1248  void MatrixSuperLU<complex<double> >::
1249  FactorizeMatrix(Matrix<T0, Prop, Storage, Allocator> & mat,
1250  bool keep_matrix)
1251  {
1252  Vector<long> Ptr; Vector<int> IndRow;
1253  Vector<complex<double> > Val; General prop;
1254  ConvertToCSC(mat, prop, Ptr, IndRow, Val, false);
1255  if (!keep_matrix)
1256  mat.Clear();
1257 
1258  FactorizeCSC(Ptr, IndRow, Val, IsSymmetricMatrix(mat));
1259  }
1260 
1261 
1262  void MatrixSuperLU<complex<double> >
1263  ::FactorizeCSC(Vector<long>& Ptr, Vector<int>& IndRow,
1264  Vector<complex<double> >& Val, bool sym)
1265  {
1266  Vector<int> glob_num(1);
1267  glob_num(0) = 0;
1268  MPI_Comm comm = MPI_COMM_SELF;
1269  FactorizeDistributedMatrix(comm, Ptr, IndRow, Val,
1270  glob_num, sym, false);
1271  }
1272 
1273 
1275  template<class Allocator2>
1276  void MatrixSuperLU<complex<double> >::
1277  Solve(Vector<complex<double>, VectFull, Allocator2>& x)
1278  {
1279  Solve(SeldonNoTrans, x);
1280  }
1281 
1282 
1284  template<class Allocator2>
1285  void MatrixSuperLU<complex<double> >::
1286  Solve(const SeldonTranspose& TransA,
1287  Vector<complex<double>, VectFull, Allocator2>& x)
1288  {
1289  Vector<int> glob_num(1);
1290  glob_num(0) = 0;
1291  MPI_Comm comm = MPI_COMM_SELF;
1292  SolveDistributed(comm, TransA, x, glob_num);
1293  }
1294 
1295 
1296  void MatrixSuperLU<complex<double> >
1297  ::Solve(const SeldonTranspose& TransA,
1298  complex<double>* x_ptr, int nrhs_)
1299  {
1300  Vector<int> glob_num(1);
1301  glob_num(0) = 0;
1302  MPI_Comm comm = MPI_COMM_SELF;
1303  SolveDistributed(comm, TransA, x_ptr, nrhs_, glob_num);
1304  }
1305 
1306 
1308  template<class Allocator2>
1309  void MatrixSuperLU<complex<double> >::
1310  Solve(Matrix<complex<double>, General, ColMajor, Allocator2>& x)
1311  {
1312  Solve(SeldonNoTrans, x);
1313  }
1314 
1315 
1317  template<class Allocator2>
1318  void MatrixSuperLU<complex<double> >::
1319  Solve(const SeldonTranspose& TransA,
1320  Matrix<complex<double>, General, ColMajor, Allocator2>& x)
1321  {
1322  Vector<int> glob_num(1);
1323  glob_num(0) = 0;
1324  MPI_Comm comm = MPI_COMM_SELF;
1325  SolveDistributed(comm, TransA, x, glob_num);
1326  }
1327 
1328 
1329  void MatrixSuperLU<complex<double> >::
1330  FactorizeDistributedMatrix(MPI_Comm& comm_facto, Vector<long>& Ptr,
1331  Vector<int>& IndRow, Vector<complex<double> >& Val,
1332  const Vector<int>& glob_number,
1333  bool sym, bool keep_matrix)
1334  {
1335  Vector<int> PtrInt(Ptr.GetM());
1336  for (int i = 0; i < Ptr.GetM(); i++)
1337  PtrInt(i) = Ptr(i);
1338 
1339  FactorizeParallel(comm_facto, PtrInt, IndRow, Val,
1340  glob_number, sym, keep_matrix);
1341  }
1342 
1343 
1344  void MatrixSuperLU<complex<double> >::
1345  FactorizeDistributedMatrix(MPI_Comm& comm_facto, Vector<int64_t>& Ptr,
1346  Vector<int64_t>& IndRow, Vector<complex<double> >& Val,
1347  const Vector<int>& glob_number,
1348  bool sym, bool keep_matrix)
1349  {
1350  FactorizeParallel(comm_facto, Ptr, IndRow, Val,
1351  glob_number, sym, keep_matrix);
1352  }
1353 
1354 
1355  template<class Tint>
1356  void MatrixSuperLU<complex<double> >::
1357  FactorizeParallel(MPI_Comm& comm_facto,
1358  Vector<Tint>&, Vector<Tint>&,
1359  Vector<complex<double> >&,
1360  const Vector<int>& glob_num,
1361  bool sym, bool keep_matrix)
1362  {
1363  cout << "Not available for this type of integers " << endl;
1364  cout << "Size of Tint = " << sizeof(Tint) << endl;
1365  abort();
1366  }
1367 
1368 
1369  void MatrixSuperLU<complex<double> >::
1370  FactorizeParallel(MPI_Comm& comm_facto,
1371  Vector<superlu_int_t>& Ptr, Vector<superlu_int_t>& Row,
1372  Vector<complex<double> >& Val,
1373  const Vector<int>& glob_num,
1374  bool sym, bool keep_matrix)
1375  {
1376 #ifndef SELDON_WITH_SUPERLU_DOUBLE
1377  // previous factorization is cleared if present
1378  Clear();
1379 
1380  if (!sym)
1381  {
1382  cout << "Problem with option TRANS of SuperLU" << endl;
1383  abort();
1384  }
1385 
1386  // m_loc : local number of rows
1387  // N : global number of rows
1388  int m_loc_ = Ptr.GetM()-1;
1389  int N = m_loc_; nloc = m_loc_;
1390  MPI_Allreduce(&m_loc_, &N, 1, MPI_INTEGER, MPI_SUM, comm_facto);
1391  int_t m_loc = m_loc_;
1392 
1393  // structures are initialized with N
1394  int_t panel_size, relax;
1395  Init(N, panel_size, relax);
1396  superlu::zScalePermstructInit(N, N, &ScalePermstruct);
1397  superlu::zLUstructInit(n, &LUstruct);
1398 
1399  // 1-D grid
1400  MPI_Comm_size(comm_facto, &nprow);
1401  npcol = 1;
1402 
1403  // initialize the superlu process grid
1404  superlu::superlu_gridinit(comm_facto, nprow, npcol, &grid);
1405 
1406  // fills the superlu structure
1407  // global numbers are assumed to be consecutive, we provide the first row number
1408  int_t fst_row = glob_num(0);
1409  int_t nnz_loc = Row.GetM();
1410  //superlu::SuperMatrix A;
1411  superlu::
1412  zCreate_CompRowLoc_Matrix_dist(&A, n, n, nnz_loc, m_loc, fst_row,
1413  reinterpret_cast<superlu::doublecomplex*>
1414  (Val.GetData()),
1415  reinterpret_cast<int_t*>(Row.GetData()),
1416  reinterpret_cast<int_t*>(Ptr.GetData()),
1417  superlu::SLU_NR_loc, superlu::SLU_Z,
1418  superlu::SLU_GE);
1419 
1420  // completes factorization
1421  int_t nrhs = 0;
1422  options.Trans = superlu::NOTRANS;
1423  options.Fact = superlu::DOFACT;
1424  superlu::
1425  pzgssvx(&options, &A, &ScalePermstruct,
1426  NULL, m_loc, nrhs, &grid,
1427  &LUstruct, &SOLVEstruct, NULL, &stat, &info_facto);
1428 
1429  Ptr.Nullify(); Row.Nullify(); Val.Nullify();
1430 #else
1431  cout << "Recompile Seldon without SELDON_WITH_SUPERLU_DOUBLE" << endl;
1432  cout << "Double and complex<double> cannot be handled simultaneously" << endl;
1433  abort();
1434 #endif
1435  }
1436 
1437 
1438  template<class Allocator2>
1439  void MatrixSuperLU<complex<double> >::
1440  SolveDistributed(MPI_Comm& comm_facto,
1441  const SeldonTranspose& TransA,
1442  Vector<complex<double>, VectFull, Allocator2>& x,
1443  const Vector<int>& glob_num)
1444  {
1445  SolveDistributed(comm_facto, TransA, x.GetData(), 1, glob_num);
1446  }
1447 
1448 
1449  template<class Allocator2>
1450  void MatrixSuperLU<complex<double> >::
1451  SolveDistributed(MPI_Comm& comm_facto,
1452  const SeldonTranspose& TransA,
1453  Matrix<complex<double>, General, ColMajor, Allocator2>& x,
1454  const Vector<int>& glob_num)
1455  {
1456  SolveDistributed(comm_facto, TransA, x.GetData(), x.GetN(), glob_num);
1457  }
1458 
1459 
1460  void MatrixSuperLU<complex<double> >::
1461  SolveDistributed(MPI_Comm& comm_facto,
1462  const SeldonTranspose& TransA,
1463  complex<double>* x_ptr, int nrhs_,
1464  const IVect& glob_num)
1465  {
1466 #ifndef SELDON_WITH_SUPERLU_DOUBLE
1467  options.Fact = superlu::FACTORED;
1468  // inverting transpose because we have provided columns instead of rows
1469  if (TransA.NoTrans())
1470  options.Trans = superlu::TRANS;
1471  else
1472  options.Trans = superlu::NOTRANS;
1473 
1474  // TRANS does not work, we put NOTRANS
1475  options.Trans = superlu::NOTRANS;
1476  Vector<double> berr(nrhs_);
1477  for (int k = 0; k < nrhs_; k++)
1478  {
1479  // workaround, loop over right hand sides
1480  int_t nrhs = 1; int_t info;
1481  //int_t nrhs = nrhs_; int_t info;
1482  superlu::
1483  pzgssvx(&options, &A, &ScalePermstruct,
1484  reinterpret_cast<superlu::doublecomplex*>(x_ptr),
1485  n, nrhs, &grid, &LUstruct, &SOLVEstruct,
1486  berr.GetData(), &stat, &info);
1487 
1488  x_ptr += nloc;
1489  }
1490 #endif
1491  }
1492 #endif
1493 
1494 
1495  /***************************
1496  * GetLU/SolveLU functions *
1497  ***************************/
1498 
1499 
1501  template<class MatrixSparse, class T>
1502  void GetLU(MatrixSparse& A, MatrixSuperLU<T>& mat_lu, bool keep_matrix, T& x)
1503  {
1504  mat_lu.FactorizeMatrix(A, keep_matrix);
1505  }
1506 
1507 
1509  template<class MatrixSparse, class T>
1510  void GetLU(MatrixSparse& A, MatrixSuperLU<T>& mat_lu,
1511  bool keep_matrix, complex<T>& x)
1512  {
1513  throw WrongArgument("GetLU(Matrix<complex<T> >& A, MatrixSuperLU<T>& mat_lu, bool)",
1514  "The LU matrix must be complex");
1515  }
1516 
1517 
1519  template<class MatrixSparse, class T>
1520  void GetLU(MatrixSparse& A, MatrixSuperLU<complex<T> >& mat_lu, bool keep_matrix, T& x)
1521  {
1522  throw WrongArgument("GetLU(Matrix<T>& A, MatrixSuperLU<complex<T> >& mat_lu, bool)",
1523  "The sparse matrix must be complex");
1524  }
1525 
1526 
1528  template<class T0, class Prop, class Storage, class Allocator, class T>
1529  void GetLU(Matrix<T0, Prop, Storage, Allocator>& A, MatrixSuperLU<T>& mat_lu,
1530  bool keep_matrix)
1531  {
1532  // we check if the type of non-zero entries of matrix A
1533  // and of the SuperLU object (T) are different
1534  // we call one of the GetLUs written above
1535  // such a protection avoids to compile the factorisation of a complex
1536  // matrix with a real SuperLU object
1537  typename Matrix<T0, Prop, Storage, Allocator>::entry_type x;
1538  GetLU(A, mat_lu, keep_matrix, x);
1539  }
1540 
1541 
1543  template<class T, class Allocator>
1544  void SolveLU(MatrixSuperLU<T>& mat_lu, Vector<T, VectFull, Allocator>& x)
1545  {
1546  mat_lu.Solve(x);
1547  }
1548 
1549 
1552  template<class T, class Allocator>
1553  void SolveLU(const SeldonTranspose& TransA,
1554  MatrixSuperLU<T>& mat_lu, Vector<T, VectFull, Allocator>& x)
1555  {
1556  mat_lu.Solve(TransA, x);
1557  }
1558 
1559 
1561  template<class T, class Prop, class Allocator>
1562  void SolveLU(MatrixSuperLU<T>& mat_lu,
1563  Matrix<T, Prop, ColMajor, Allocator>& x)
1564  {
1565  mat_lu.Solve(SeldonNoTrans, x);
1566  }
1567 
1568 
1571  template<class T, class Prop, class Allocator>
1572  void SolveLU(const SeldonTranspose& TransA,
1573  MatrixSuperLU<T>& mat_lu, Matrix<T, Prop, ColMajor, Allocator>& x)
1574  {
1575  mat_lu.Solve(TransA, x);
1576  }
1577 
1578 
1580  template<class Allocator>
1581  void SolveLU(MatrixSuperLU<double>& mat_lu,
1582  Vector<complex<double>, VectFull, Allocator>& x)
1583  {
1584  Matrix<double, General, ColMajor> y(x.GetM(), 2);
1585 
1586  for (int i = 0; i < x.GetM(); i++)
1587  {
1588  y(i, 0) = real(x(i));
1589  y(i, 1) = imag(x(i));
1590  }
1591 
1592  SolveLU(mat_lu, y);
1593 
1594  for (int i = 0; i < x.GetM(); i++)
1595  x(i) = complex<double>(y(i, 0), y(i, 1));
1596  }
1597 
1598 
1600  template<class Allocator>
1601  void SolveLU(const SeldonTranspose& TransA,
1602  MatrixSuperLU<double>& mat_lu,
1603  Vector<complex<double>, VectFull, Allocator>& x)
1604  {
1605  Matrix<double, General, ColMajor> y(x.GetM(), 2);
1606 
1607  for (int i = 0; i < x.GetM(); i++)
1608  {
1609  y(i, 0) = real(x(i));
1610  y(i, 1) = imag(x(i));
1611  }
1612 
1613  SolveLU(TransA, mat_lu, y);
1614 
1615  for (int i = 0; i < x.GetM(); i++)
1616  x(i) = complex<double>(y(i, 0), y(i, 1));
1617 
1618  }
1619 
1620 
1622  template<class Allocator>
1623  void SolveLU(MatrixSuperLU<complex<double> >& mat_lu,
1624  Vector<double, VectFull, Allocator>& x)
1625  {
1626  throw WrongArgument("SolveLU(MatrixSuperLU<complex<double> >, Vector<double>)",
1627  "The result should be a complex vector");
1628  }
1629 
1630 
1632  template<class Allocator>
1633  void SolveLU(const SeldonTranspose& TransA,
1634  MatrixSuperLU<complex<double> >& mat_lu,
1635  Vector<double, VectFull, Allocator>& x)
1636  {
1637  throw WrongArgument("SolveLU(MatrixSuperLU<complex<double> >, Vector<double>)",
1638  "The result should be a complex vector");
1639  }
1640 
1641 }
1642 
1643 #define SELDON_FILE_SUPERLU_CXX
1644 #endif
Seldon::MatrixSuperLU< double >::FactorizeCSC
void FactorizeCSC(Vector< superlu_int_t > &Ptr, Vector< superlu_int_t > &IndRow, Vector< double > &Val, bool sym)
factorization of matrix given in CSC form
Definition: SuperLU.cxx:481
Seldon::SeldonTranspose
Definition: MatrixFlag.hxx:32
Seldon::ColSparse
Definition: Storage.hxx:79
Seldon::MatrixSuperLU_Base::~MatrixSuperLU_Base
~MatrixSuperLU_Base()
destructor
Definition: SuperLU.cxx:158
Seldon::MatrixSuperLU< double >::Solve
void Solve(Vector< double, VectFull, Allocator2 > &x)
Solves linear system A x = b.
Definition: SuperLU.cxx:566
Seldon::MatrixSuperLU_Base::Init
void Init(int size, int_t &panel_size, int_t &relax)
inits SuperLU computation
Definition: SuperLU.cxx:199
Seldon::MatrixSuperLU_Base::SetNumberOfThreadPerNode
void SetNumberOfThreadPerNode(int p)
Sets the number of threads per mpi process.
Definition: SuperLU.cxx:228
Seldon::Vector< int_t >
Seldon::IsSymmetricMatrix
bool IsSymmetricMatrix(const Matrix< T, Prop, Storage, Allocator > &A)
returns true if the matrix is symmetric
Definition: Functions_BaseInline.cxx:778
Seldon::MatrixSuperLU_Base::perm_r
Vector< int_t > perm_r
permutation array
Definition: SuperLU.hxx:92
Seldon::Matrix
Definition: SeldonHeader.hxx:226
Seldon::MatrixSuperLU_Base::HideMessages
void HideMessages()
no message from SuperLU
Definition: SuperLU.cxx:238
Seldon::MatrixSuperLU_Base::GetInfoFactorization
int GetInfoFactorization() const
returns status of factorisation
Definition: SuperLU.cxx:264
Seldon::MatrixSuperLU_Base::L
superlu::SuperMatrix L
objects of SuperLU
Definition: SuperLU.hxx:58
Seldon::MatrixSuperLU_Base::options
superlu::superlu_options_t options
options
Definition: SuperLU.hxx:76
Seldon::MatrixSuperLU_Base::MatrixSuperLU_Base
MatrixSuperLU_Base()
default constructor
Definition: SuperLU.cxx:119
Seldon::VectFull
Definition: StorageInline.cxx:74
Seldon::MatrixSuperLU_Base::GetColPermutation
const Vector< int_t > & GetColPermutation() const
Returns the permutation of columns.
Definition: SuperLU.cxx:292
Seldon::MatrixSuperLU_Base::stat
superlu::SuperLUStat_t stat
statistics
Definition: SuperLU.hxx:77
Seldon::MatrixSuperLU_Base::Clear
void Clear()
same effect as a call to the destructor
Definition: SuperLU.cxx:166
Seldon::General
Definition: Properties.hxx:26
Seldon::MatrixSuperLU
empty matrix
Definition: SuperLU.hxx:125
Seldon::Vector< T, VectFull, Allocator >
Full vector class.
Definition: Vector.hxx:88
Seldon::GetLU
void GetLU(Matrix< T0, Prop0, Storage0, Allocator0 > &A)
Returns the LU factorization of a matrix.
Definition: Functions_Matrix.cxx:2073
Seldon::MatrixSuperLU_Base::ShowMessages
void ShowMessages()
allows messages from SuperLU
Definition: SuperLU.cxx:246
Seldon::MatrixSuperLU_Base::GetRowPermutation
const Vector< int_t > & GetRowPermutation() const
Returns the permutation of rows.
Definition: SuperLU.cxx:278
Seldon::MatrixSuperLU< double >::GetLU
void GetLU(Matrix< double, Prop, ColSparse, Allocator > &Lmat, Matrix< double, Prop, ColSparse, Allocator > &Umat, bool permuted=true)
Returns the LU factorization.
Definition: SuperLU.cxx:347
Seldon::WrongArgument
Definition: Errors.hxx:76
Seldon::MatrixSuperLU_Base::SelectOrdering
void SelectOrdering(int type)
selects ordering to use in the interfaced solver
Definition: SuperLU.cxx:299
Seldon::RowSparse
Definition: Storage.hxx:91
Seldon::MatrixSuperLU_Base::n
int_t n
number of rows
Definition: SuperLU.hxx:95
Seldon
Seldon namespace.
Definition: Array.cxx:24
Seldon::MatrixSuperLU_Base
class interfacing SuperLU functions
Definition: SuperLU.hxx:53
Seldon::MatrixSuperLU< double >
class interfacing SuperLU functions in double precision
Definition: SuperLU.hxx:132
Seldon::Matrix< T, Prop, ColMajor, Allocator >
Column-major full-matrix class.
Definition: Matrix_Pointers.hxx:176
Seldon::ColMajor
Definition: Storage.hxx:33