Wsmp.cxx
1 // Copyright (C) 2015-2015 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 #ifndef SELDON_FILE_WSMP_CXX
20 
21 #include "Wsmp.hxx"
22 
23 namespace Seldon
24 {
25 
26  template<>
27  void MatrixWsmp<double>::CallWssmp(int* n_, int* ia, int* ja, double* avals, double* diag,
28  int* perm, int* invp, double* b, int* ldb, int* nrhs,
29  double* aux, int* naux, int* mrp, int* iparam, double* dparam)
30  {
31 #ifdef SELDON_WITH_MPI
32  pwssmp_(n_, ia, ja, avals, diag, perm, invp,
33  b, ldb, nrhs, aux, naux, mrp, iparam, dparam);
34 #else
35  wssmp_(n_, ia, ja, avals, diag, perm, invp,
36  b, ldb, nrhs, aux, naux, mrp, iparam, dparam);
37 #endif
38  }
39 
40 
41  template<>
42  void MatrixWsmp<complex<double> >
43  ::CallWssmp(int* n_, int* ia, int* ja, complex<double>* avals, complex<double>* diag,
44  int* perm, int* invp, complex<double>* b, int* ldb, int* nrhs,
45  complex<double>* aux, int* naux, int* mrp, int* iparam, double* dparam)
46  {
47 #ifdef SELDON_WITH_MPI
48  pzssmp_(n_, ia, ja, avals, diag, perm, invp,
49  b, ldb, nrhs, aux, naux, mrp, iparam, dparam);
50 #else
51  zssmp_(n_, ia, ja, avals, diag, perm, invp,
52  b, ldb, nrhs, aux, naux, mrp, iparam, dparam);
53 #endif
54  }
55 
56 
57  template<>
58  void MatrixWsmp<double>::CallWgsmp(int* n_, int* ia, int* ja, double* avals,
59  double* b, int* ldb, int* nrhs,
60  double* rmisc, int* iparam, double* dparam)
61  {
62 #ifdef SELDON_WITH_MPI
63  pwgsmp_(n_, ia, ja, avals, b, ldb, nrhs, rmisc, iparam, dparam);
64 #else
65  wgsmp_(n_, ia, ja, avals, b, ldb, nrhs, rmisc, iparam, dparam);
66 #endif
67  }
68 
69 
70  template<>
71  void MatrixWsmp<complex<double> >
72  ::CallWgsmp(int* n_, int* ia, int* ja, complex<double>* avals,
73  complex<double>* b, int* ldb, int* nrhs,
74  complex<double>* rmisc, int* iparam, double* dparam)
75  {
76 #ifdef SELDON_WITH_MPI
77  pzgsmp_(n_, ia, ja, avals, b, ldb, nrhs, rmisc, iparam, dparam);
78 #else
79  zgsmp_(n_, ia, ja, avals, b, ldb, nrhs, rmisc, iparam, dparam);
80 #endif
81  }
82 
83 
84  template<class T>
85  MatrixWsmp<T>::MatrixWsmp() : iparm(64), dparm(64)
86  {
87  iparm.Zero();
88  dparm.Zero();
89  n = 0;
90  refine_solution = false;
91  use_pivoting = true;
92  cholesky = false;
93  symmetric = false;
94  distributed = false;
95  threshold_pivot = 0.001;
96 
97  int nb_threads = 1;
98  wsetmaxthrds_(&nb_threads);
99  }
100 
101 
102  template<class T>
103  MatrixWsmp<T>::~MatrixWsmp()
104  {
105  Clear();
106  }
107 
108 
109  template<class T>
110  bool MatrixWsmp<T>::UseInteger8() const
111  {
112  return false;
113  }
114 
115 
116  template<class T>
117  void MatrixWsmp<T>::Clear()
118  {
119  if (n > 0)
120  {
121  wsmp_clear_();
122  Ptr.Clear(); IndRow.Clear(); Val.Clear();
123  permut.Clear(); inverse_permut.Clear();
124  n = 0;
125  }
126  }
127 
128 
129  template<class T>
130  void MatrixWsmp<T>::ShowMessages()
131  {
132  }
133 
134 
135  template<class T>
136  void MatrixWsmp<T>::HideMessages()
137  {
138  }
139 
140 
141  template<class T>
142  void MatrixWsmp<T>::ShowFullHistory()
143  {
144  }
145 
146 
147  template<class T>
149  {
150  threshold_pivot = eps;
151  }
152 
153 
154  template<class T>
155  void MatrixWsmp<T>::SetCholeskyFacto(bool chol)
156  {
157  cholesky = chol;
158  }
159 
160 
161  template<class T>
162  int MatrixWsmp<T>::GetInfoFactorization() const
163  {
164  return iparm(63);
165  }
166 
167 
168  template<class T>
170  {
171  refine_solution = true;
172  }
173 
174 
175  template<class T>
177  {
178  refine_solution = false;
179  }
180 
181 
182  template<class T>
183  size_t MatrixWsmp<T>::GetMemorySize() const
184  {
185  if (n > 0)
186  return size_t(iparm(22))*1024*sizeof(T);
187 
188  return 0;
189  }
190 
191 
192  template<class T>
194  {
195  wsetmaxthrds_(&nb_threads);
196  }
197 
198 
200  template<class T> template<class Storage, class Allocator>
201  void MatrixWsmp<T>
203  bool keep_matrix)
204  {
205  Clear();
206 
207  Symmetric prop;
208  ConvertToCSR(mat, prop, Ptr, IndRow, Val);
209  if (!keep_matrix)
210  mat.Clear();
211 
212  FactorizeSymmetric();
213  }
214 
215 
216  template<class T>
218  Vector<T>& Val_, bool sym)
219  {
220  Ptr.SetData(Ptr_);
221  IndRow.SetData(IndRow_);
222  Val.SetData(Val_);
223  Ptr_.Nullify(); IndRow_.Nullify(); Val_.Nullify();
224 
225  if (sym)
226  FactorizeSymmetric();
227  else
228  FactorizeUnsymmetric();
229  }
230 
231 
232  template<class T>
233  void MatrixWsmp<T>::FactorizeSymmetric()
234  {
235  distributed = false;
236  symmetric = true;
237  n = Ptr.GetM()-1;
238 
239  wsmp_initialize_();
240 
241  iparm(0) = 0;
242  iparm(1) = 0;
243  iparm(2) = 0;
244 
245  // initialisation step
246  int nrhs = 0, naux = 0, mrp = 0;
247  CallWssmp(&n, NULL, NULL, NULL, NULL, NULL, NULL,
248  NULL, &n, &nrhs, NULL, &naux, &mrp, iparm.GetData(), dparm.GetData());
249 
250  // CSC/CSR format
251  iparm(3) = 0;
252  // C-style index
253  iparm(4) = 0;
254 
255  // wsmp can reuse the arrays IndRow and Val
256  iparm(13) = 3;
257 
258  // pivot threshold
259  dparm(10) = threshold_pivot;
260 
261  if (cholesky)
262  {
263  // L L^T factorization
264  iparm(30) = 0;
265  }
266  else
267  {
268  // L D L^T factorization
269  if (use_pivoting)
270  iparm(30) = 2;
271  else
272  iparm(30) = 1;
273 
274  if (IsComplexNumber(T()))
275  iparm(30) += 2;
276  }
277 
278  // ordering step, symbolic step
279  permut.Reallocate(n);
280  inverse_permut.Reallocate(n);
281 
282  iparm(1) = 1;
283  iparm(2) = 2;
284  CallWssmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(),
285  NULL, permut.GetData(), inverse_permut.GetData(), NULL, &n, &nrhs,
286  NULL, &naux, &mrp, iparm.GetData(), dparm.GetData());
287 
288  // then factorization
289  iparm(1) = 3;
290  iparm(2) = 3;
291  CallWssmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(),
292  NULL, permut.GetData(), inverse_permut.GetData(), NULL, &n, &nrhs,
293  NULL, &naux, &mrp, iparm.GetData(), dparm.GetData());
294 
295  }
296 
297 
299  template<class T> template<class Storage, class Allocator>
300  void MatrixWsmp<T>
302  bool keep_matrix)
303  {
304  Clear();
305 
306  General prop;
307  ConvertToCSR(mat, prop, Ptr, IndRow, Val);
308  if (!keep_matrix)
309  mat.Clear();
310 
311  FactorizeUnsymmetric();
312  }
313 
314 
315  template<class T>
317  {
318  distributed = false;
319  symmetric = false;
320  n = Ptr.GetM()-1;
321 
322  wsmp_initialize_();
323 
324  iparm(0) = 0;
325  iparm(1) = 0;
326  iparm(2) = 0;
327 
328  // initialisation step
329  int nrhs = 0; T darray;
330  CallWgsmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(),
331  NULL, &n, &nrhs, &darray, iparm.GetData(), dparm.GetData());
332 
333  // CSC/CSR format
334  iparm(3) = 0;
335  // C-style index
336  iparm(4) = 0;
337 
338  // pivot threshold
339  dparm(10) = threshold_pivot;
340 
341  // ordering step, symbolic step
342  iparm(1) = 1;
343  iparm(2) = 1;
344  CallWgsmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(),
345  NULL, &n, &nrhs, &darray, iparm.GetData(), dparm.GetData());
346 
347  // then factorization
348  iparm(1) = 2;
349  iparm(2) = 2;
350  CallWgsmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(),
351  NULL, &n, &nrhs, &darray, iparm.GetData(), dparm.GetData());
352 
353  }
354 
355 
356  template<class T>
357  void MatrixWsmp<T>::Solve(Vector<T>& b)
358  {
359  Solve(SeldonNoTrans, b);
360  }
361 
362 
363  template<class T>
364  void MatrixWsmp<T>::Solve(const SeldonTranspose& trans, Vector<T>& b)
365  {
366  Solve(trans, b.GetData(), 1);
367  }
368 
369 
370  template<class T>
371  void MatrixWsmp<T>::Solve(const SeldonTranspose& trans, T* x_ptr, int nrhs)
372  {
373  int naux = 0, mrp = 0;
374 
375  // solve
376  if (symmetric)
377  {
378  iparm(1) = 4;
379  iparm(2) = 4;
380  if (cholesky)
381  {
382  if (trans.Trans())
383  iparm(29) = 2;
384  else
385  iparm(29) = 1;
386  }
387  else
388  {
389  if (refine_solution)
390  iparm(2) = 5;
391  }
392 
393  CallWssmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(),
394  NULL, permut.GetData(), inverse_permut.GetData(), x_ptr, &n, &nrhs,
395  NULL, &naux, &mrp, iparm.GetData(), dparm.GetData());
396  }
397  else
398  {
399  iparm(1) = 3;
400  iparm(2) = 3;
401  if (refine_solution)
402  iparm(2) = 4;
403 
404  if (trans.Trans())
405  {
406  for (int i = 0; i < n*nrhs; i++)
407  x_ptr[i] = conjugate(x_ptr[i]);
408 
409  iparm(29) = 4;
410  }
411  else
412  iparm(29) = 0;
413 
414  CallWgsmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(),
415  x_ptr, &n, &nrhs,
416  NULL, iparm.GetData(), dparm.GetData());
417 
418  if (trans.Trans())
419  for (int i = 0; i < n*nrhs; i++)
420  x_ptr[i] = conjugate(x_ptr[i]);
421  }
422  }
423 
424 
425  template<class T>
426  void MatrixWsmp<T>::Solve(const SeldonTranspose& trans, Matrix<T, General, ColMajor>& b)
427  {
428  Solve(trans, b.GetData(), b.GetN());
429  }
430 
431 
432 #ifdef SELDON_WITH_MPI
433  template<class T>
434  void MatrixWsmp<T>::
435  FactorizeDistributedMatrix(MPI_Comm& comm_facto, Vector<long>& Ptr,
436  Vector<int>& IndRow, Vector<T>& Val,
437  const Vector<int>& glob_number,
438  bool sym, bool keep_matrix)
439  {
440  Vector<int> PtrInt(Ptr.GetM());
441  for (int i = 0; i < Ptr.GetM(); i++)
442  PtrInt(i) = Ptr(i);
443 
444  FactorizeParallel(comm_facto, PtrInt, IndRow, Val,
445  glob_number, sym, keep_matrix);
446  }
447 
448 
449  template<class T>
450  void MatrixWsmp<T>::
451  FactorizeDistributedMatrix(MPI_Comm& comm_facto, Vector<int64_t>& Ptr,
452  Vector<int64_t>& IndRow, Vector<T>& Val,
453  const Vector<int>& glob_number,
454  bool sym, bool keep_matrix)
455  {
456  FactorizeParallel(comm_facto, Ptr, IndRow, Val,
457  glob_number, sym, keep_matrix);
458  }
459 
460 
461  template<class T> template<class Tint>
462  void MatrixWsmp<T>::
463  FactorizeParallel(MPI_Comm& comm_facto,
464  Vector<Tint>&, Vector<Tint>&,
465  Vector<T>&,
466  const Vector<int>& glob_num,
467  bool sym, bool keep_matrix)
468  {
469  cout << "Not available for this type of integers " << endl;
470  cout << "Size of Tint = " << sizeof(Tint) << endl;
471  abort();
472  }
473 
474 
476  template<class T>
477  void MatrixWsmp<T>::
478  FactorizeParallel(MPI_Comm& comm_facto,
479  Vector<int>& Ptr_, Vector<int>& IndRow_,
480  Vector<T>& Val_, const Vector<int>& glob_number,
481  bool sym, bool keep_matrix)
482  {
483  // we clear previous factorization if present
484  Clear();
485 
486  Ptr.SetData(Ptr_);
487  IndRow.SetData(IndRow_);
488  Val.SetData(Val_);
489  Ptr_.Nullify(); IndRow_.Nullify(); Val_.Nullify();
490 
491  // size of local system
492  n = Ptr.GetM() - 1;
493  symmetric = sym;
494  distributed = true;
495 
496  wsmp_initialize_();
497 
498  // finding the size of the overall system
499  int nmax = 0, N = 0;
500  for (int i = 0; i < glob_number.GetM(); i++)
501  nmax = max(glob_number(i)+1, nmax);
502 
503  MPI_Allreduce(&nmax, &N, 1, MPI_INTEGER, MPI_MAX, comm_facto);
504 
505  int nrhs = 0, naux = 0, mrp = 0;
506  if (sym)
507  {
508  CallWssmp(&n, NULL, NULL, NULL, NULL, NULL, NULL,
509  NULL, &n, &nrhs, NULL, &naux, &mrp, iparm.GetData(), dparm.GetData());
510 
511  // CSC/CSR format
512  iparm(3) = 0;
513  // C-style index
514  iparm(4) = 0;
515 
516  // wsmp can reuse the arrays IndRow and Val
517  iparm(13) = 3;
518 
519  // pivot threshold
520  dparm(10) = threshold_pivot;
521 
522  if (cholesky)
523  {
524  // L L^T factorization
525  iparm(30) = 0;
526  }
527  else
528  {
529  // L D L^T factorization
530  if (use_pivoting)
531  iparm(30) = 2;
532  else
533  iparm(30) = 1;
534 
535  if (IsComplexNumber(T()))
536  iparm(30) += 2;
537  }
538 
539  permut.Reallocate(N);
540  inverse_permut.Reallocate(N);
541 
542  iparm(1) = 1;
543  iparm(2) = 3;
544  CallWssmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(), NULL,
545  permut.GetData(), inverse_permut.GetData(), NULL,
546  &n, &nrhs, NULL, &naux, &mrp, iparm.GetData(), dparm.GetData());
547  }
548  else
549  {
550 
551  iparm(0) = 0;
552  iparm(1) = 0;
553  iparm(2) = 0;
554 
555  // initialisation step
556  T darray;
557  CallWgsmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(),
558  NULL, &n, &nrhs, &darray, iparm.GetData(), dparm.GetData());
559 
560  // CSC/CSR format
561  iparm(3) = 1;
562  // C-style index
563  iparm(4) = 0;
564 
565  // pivot threshold
566  dparm(10) = threshold_pivot;
567 
568  // ordering step and factorization
569  iparm(1) = 1;
570  iparm(2) = 2;
571  CallWgsmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(),
572  NULL, &n, &nrhs, &darray, iparm.GetData(), dparm.GetData());
573 
574  }
575  }
576 
577 
579  template<class T>
580  void MatrixWsmp<T>::SolveDistributed(MPI_Comm& comm_facto,
581  const SeldonTranspose& TransA,
582  Vector<T>& x, const Vector<int>& glob_num)
583  {
584  SolveDistributed(comm_facto, TransA, x.GetData(), 1, glob_num);
585  }
586 
587 
588  template<class T>
589  void MatrixWsmp<T>::SolveDistributed(MPI_Comm& comm_facto,
590  const SeldonTranspose& TransA,
591  T* x_ptr, int nrhs, const Vector<int>& glob_num)
592  {
593  int naux = 0, mrp = 0;
594  if (symmetric)
595  {
596  iparm(1) = 4;
597  iparm(2) = 4;
598  if (cholesky)
599  {
600  if (TransA.Trans())
601  iparm(29) = 2;
602  else
603  iparm(29) = 1;
604  }
605  else
606  {
607  if (refine_solution)
608  iparm(2) = 5;
609  }
610 
611  CallWssmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(), NULL,
612  permut.GetData(), inverse_permut.GetData(), x_ptr,
613  &n, &nrhs, NULL, &naux, &mrp, iparm.GetData(), dparm.GetData());
614  }
615  else
616  {
617  iparm(1) = 3;
618  iparm(2) = 3;
619  if (refine_solution)
620  iparm(2) = 4;
621 
622  if (TransA.Trans())
623  {
624  for (int i = 0; i < n*nrhs; i++)
625  x_ptr[i] = conjugate(x_ptr[i]);
626 
627  iparm(29) = 4;
628  }
629  else
630  iparm(29) = 0;
631 
632  CallWgsmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(),
633  x_ptr, &n, &nrhs,
634  NULL, iparm.GetData(), dparm.GetData());
635 
636  if (TransA.Trans())
637  for (int i = 0; i < n*nrhs; i++)
638  x_ptr[i] = conjugate(x_ptr[i]);
639  }
640  }
641 
642 
644  template<class T>
645  void MatrixWsmp<T>::SolveDistributed(MPI_Comm& comm_facto,
646  const SeldonTranspose& TransA,
647  Matrix<T, General, ColMajor>& x,
648  const Vector<int>& glob_num)
649  {
650  SolveDistributed(comm_facto, TransA, x.GetData(), x.GetN(), glob_num);
651  }
652 #endif
653 
654 
656  template<class T, class Prop, class Storage, class Allocator>
658  bool keep_matrix)
659  {
660  mat_lu.FactorizeMatrix(A, keep_matrix);
661  }
662 
663 
665  template<class T, class Allocator>
667  {
668  mat_lu.Solve(x);
669  }
670 
671 
674  template<class T, class Allocator>
675  void SolveLU(const SeldonTranspose& TransA,
677  {
678  mat_lu.Solve(TransA, x);
679  }
680 
681 
683  template<class T, class Prop, class Allocator>
684  void SolveLU(MatrixWsmp<T>& mat_lu,
686  {
687  mat_lu.Solve(SeldonNoTrans, x);
688  }
689 
690 
693  template<class T, class Prop, class Allocator>
694  void SolveLU(const SeldonTranspose& TransA,
696  {
697  mat_lu.Solve(TransA, x);
698  }
699 
700 
702  template<class Allocator>
703  void SolveLU(MatrixWsmp<double>& mat_lu,
704  Vector<complex<double>, VectFull, Allocator>& x)
705  {
706  Matrix<double, General, ColMajor> y(x.GetM(), 2);
707 
708  for (int i = 0; i < x.GetM(); i++)
709  {
710  y(i, 0) = real(x(i));
711  y(i, 1) = imag(x(i));
712  }
713 
714  SolveLU(mat_lu, y);
715 
716  for (int i = 0; i < x.GetM(); i++)
717  x(i) = complex<double>(y(i, 0), y(i, 1));
718  }
719 
720 
722  template<class Allocator>
723  void SolveLU(const SeldonTranspose& TransA,
724  MatrixWsmp<double>& mat_lu,
725  Vector<complex<double>, VectFull, Allocator>& x)
726  {
727  Matrix<double, General, ColMajor> y(x.GetM(), 2);
728 
729  for (int i = 0; i < x.GetM(); i++)
730  {
731  y(i, 0) = real(x(i));
732  y(i, 1) = imag(x(i));
733  }
734 
735  SolveLU(TransA, mat_lu, y);
736 
737  for (int i = 0; i < x.GetM(); i++)
738  x(i) = complex<double>(y(i, 0), y(i, 1));
739 
740  }
741 
742 
744  template<class Allocator>
745  void SolveLU(MatrixWsmp<complex<double> >& mat_lu,
747  {
748  throw WrongArgument("SolveLU(MatrixSuperLU<complex<double> >, Vector<double>)",
749  "The result should be a complex vector");
750  }
751 
752 
754  template<class Allocator>
755  void SolveLU(const SeldonTranspose& TransA,
756  MatrixWsmp<complex<double> >& mat_lu,
758  {
759  throw WrongArgument("SolveLU(MatrixSuperLU<complex<double> >, Vector<double>)",
760  "The result should be a complex vector");
761  }
762 
763 }
764 
765 #define SELDON_FILE_WSMP_CXX
766 #endif
Seldon::SeldonTranspose
Definition: MatrixFlag.hxx:32
Seldon::MatrixWsmp::SetPivotThreshold
void SetPivotThreshold(double)
Sets the threshold for pivot.
Definition: Wsmp.cxx:148
Seldon::Vector< int >
Seldon::Matrix
Definition: SeldonHeader.hxx:226
Seldon::VectFull
Definition: StorageInline.cxx:74
Seldon::MatrixWsmp::DoNotRefineSolution
void DoNotRefineSolution()
Tells to the direct solver that no refinement is required.
Definition: Wsmp.cxx:176
Seldon::IsComplexNumber
bool IsComplexNumber(const T &number)
Returns true for a complex number.
Definition: CommonInline.cxx:293
Seldon::MatrixWsmp::RefineSolution
void RefineSolution()
Tells to the direct solver that refinement is required.
Definition: Wsmp.cxx:169
Seldon::General
Definition: Properties.hxx:26
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::MatrixUmfPack_Base::n
int n
number of rows in the matrix
Definition: UmfPack.hxx:45
Seldon::MatrixWsmp::SetNumberOfThreadPerNode
void SetNumberOfThreadPerNode(int nb_threads)
Sets the number of threads per mpi process.
Definition: Wsmp.cxx:193
Seldon::WrongArgument
Definition: Errors.hxx:76
Seldon::MatrixWsmp
Definition: Wsmp.hxx:58
Seldon::Symmetric
Definition: Properties.hxx:30
Seldon
Seldon namespace.
Definition: Array.cxx:24
Seldon::MatrixWsmp::FactorizeMatrix
void FactorizeMatrix(Matrix< T, Symmetric, Storage, Allocator > &mat, bool keep_matrix=false)
factorizes a symmetric matrix
Definition: Wsmp.cxx:202
Seldon::Matrix< T, Prop, ColMajor, Allocator >
Column-major full-matrix class.
Definition: Matrix_Pointers.hxx:176