19 #ifndef SELDON_FILE_WSMP_CXX
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)
31 #ifdef SELDON_WITH_MPI
32 pwssmp_(n_, ia, ja, avals, diag, perm, invp,
33 b, ldb, nrhs, aux, naux, mrp, iparam, dparam);
35 wssmp_(n_, ia, ja, avals, diag, perm, invp,
36 b, ldb, nrhs, aux, naux, mrp, iparam, dparam);
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)
47 #ifdef SELDON_WITH_MPI
48 pzssmp_(n_, ia, ja, avals, diag, perm, invp,
49 b, ldb, nrhs, aux, naux, mrp, iparam, dparam);
51 zssmp_(n_, ia, ja, avals, diag, perm, invp,
52 b, ldb, nrhs, aux, naux, mrp, iparam, dparam);
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)
62 #ifdef SELDON_WITH_MPI
63 pwgsmp_(n_, ia, ja, avals, b, ldb, nrhs, rmisc, iparam, dparam);
65 wgsmp_(n_, ia, ja, avals, b, ldb, nrhs, rmisc, iparam, dparam);
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)
76 #ifdef SELDON_WITH_MPI
77 pzgsmp_(n_, ia, ja, avals, b, ldb, nrhs, rmisc, iparam, dparam);
79 zgsmp_(n_, ia, ja, avals, b, ldb, nrhs, rmisc, iparam, dparam);
85 MatrixWsmp<T>::MatrixWsmp() : iparm(64), dparm(64)
90 refine_solution =
false;
95 threshold_pivot = 0.001;
98 wsetmaxthrds_(&nb_threads);
103 MatrixWsmp<T>::~MatrixWsmp()
110 bool MatrixWsmp<T>::UseInteger8()
const
117 void MatrixWsmp<T>::Clear()
122 Ptr.Clear(); IndRow.Clear(); Val.Clear();
123 permut.Clear(); inverse_permut.Clear();
130 void MatrixWsmp<T>::ShowMessages()
136 void MatrixWsmp<T>::HideMessages()
142 void MatrixWsmp<T>::ShowFullHistory()
150 threshold_pivot = eps;
162 int MatrixWsmp<T>::GetInfoFactorization()
const
171 refine_solution =
true;
178 refine_solution =
false;
186 return size_t(iparm(22))*1024*
sizeof(T);
195 wsetmaxthrds_(&nb_threads);
200 template<
class T>
template<
class Storage,
class Allocator>
208 ConvertToCSR(mat, prop, Ptr, IndRow, Val);
212 FactorizeSymmetric();
221 IndRow.SetData(IndRow_);
223 Ptr_.Nullify(); IndRow_.Nullify(); Val_.Nullify();
226 FactorizeSymmetric();
228 FactorizeUnsymmetric();
233 void MatrixWsmp<T>::FactorizeSymmetric()
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());
259 dparm(10) = threshold_pivot;
279 permut.Reallocate(n);
280 inverse_permut.Reallocate(n);
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());
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());
299 template<
class T>
template<
class Storage,
class Allocator>
307 ConvertToCSR(mat, prop, Ptr, IndRow, Val);
311 FactorizeUnsymmetric();
329 int nrhs = 0; T darray;
330 CallWgsmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(),
331 NULL, &n, &nrhs, &darray, iparm.GetData(), dparm.GetData());
339 dparm(10) = threshold_pivot;
344 CallWgsmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(),
345 NULL, &n, &nrhs, &darray, iparm.GetData(), dparm.GetData());
350 CallWgsmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(),
351 NULL, &n, &nrhs, &darray, iparm.GetData(), dparm.GetData());
357 void MatrixWsmp<T>::Solve(Vector<T>& b)
359 Solve(SeldonNoTrans, b);
364 void MatrixWsmp<T>::Solve(
const SeldonTranspose& trans, Vector<T>& b)
366 Solve(trans, b.GetData(), 1);
371 void MatrixWsmp<T>::Solve(
const SeldonTranspose& trans, T* x_ptr,
int nrhs)
373 int naux = 0, mrp = 0;
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());
406 for (
int i = 0; i < n*nrhs; i++)
407 x_ptr[i] = conjugate(x_ptr[i]);
414 CallWgsmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(),
416 NULL, iparm.GetData(), dparm.GetData());
419 for (
int i = 0; i < n*nrhs; i++)
420 x_ptr[i] = conjugate(x_ptr[i]);
426 void MatrixWsmp<T>::Solve(
const SeldonTranspose& trans, Matrix<T, General, ColMajor>& b)
428 Solve(trans, b.GetData(), b.GetN());
432 #ifdef SELDON_WITH_MPI
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)
440 Vector<int> PtrInt(Ptr.GetM());
441 for (
int i = 0; i < Ptr.GetM(); i++)
444 FactorizeParallel(comm_facto, PtrInt, IndRow, Val,
445 glob_number, sym, keep_matrix);
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)
456 FactorizeParallel(comm_facto, Ptr, IndRow, Val,
457 glob_number, sym, keep_matrix);
461 template<
class T>
template<
class T
int>
463 FactorizeParallel(MPI_Comm& comm_facto,
464 Vector<Tint>&, Vector<Tint>&,
466 const Vector<int>& glob_num,
467 bool sym,
bool keep_matrix)
469 cout <<
"Not available for this type of integers " << endl;
470 cout <<
"Size of Tint = " <<
sizeof(Tint) << endl;
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)
487 IndRow.SetData(IndRow_);
489 Ptr_.Nullify(); IndRow_.Nullify(); Val_.Nullify();
500 for (
int i = 0; i < glob_number.GetM(); i++)
501 nmax = max(glob_number(i)+1, nmax);
503 MPI_Allreduce(&nmax, &N, 1, MPI_INTEGER, MPI_MAX, comm_facto);
505 int nrhs = 0, naux = 0, mrp = 0;
508 CallWssmp(&n, NULL, NULL, NULL, NULL, NULL, NULL,
509 NULL, &n, &nrhs, NULL, &naux, &mrp, iparm.GetData(), dparm.GetData());
520 dparm(10) = threshold_pivot;
539 permut.Reallocate(N);
540 inverse_permut.Reallocate(N);
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());
557 CallWgsmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(),
558 NULL, &n, &nrhs, &darray, iparm.GetData(), dparm.GetData());
566 dparm(10) = threshold_pivot;
571 CallWgsmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(),
572 NULL, &n, &nrhs, &darray, iparm.GetData(), dparm.GetData());
580 void MatrixWsmp<T>::SolveDistributed(MPI_Comm& comm_facto,
581 const SeldonTranspose& TransA,
582 Vector<T>& x,
const Vector<int>& glob_num)
584 SolveDistributed(comm_facto, TransA, x.GetData(), 1, glob_num);
589 void MatrixWsmp<T>::SolveDistributed(MPI_Comm& comm_facto,
590 const SeldonTranspose& TransA,
591 T* x_ptr,
int nrhs,
const Vector<int>& glob_num)
593 int naux = 0, mrp = 0;
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());
624 for (
int i = 0; i < n*nrhs; i++)
625 x_ptr[i] = conjugate(x_ptr[i]);
632 CallWgsmp(&n, Ptr.GetData(), IndRow.GetData(), Val.GetData(),
634 NULL, iparm.GetData(), dparm.GetData());
637 for (
int i = 0; i < n*nrhs; i++)
638 x_ptr[i] = conjugate(x_ptr[i]);
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)
650 SolveDistributed(comm_facto, TransA, x.GetData(), x.GetN(), glob_num);
656 template<
class T,
class Prop,
class Storage,
class Allocator>
665 template<
class T,
class Allocator>
674 template<
class T,
class Allocator>
678 mat_lu.Solve(TransA, x);
683 template<
class T,
class Prop,
class Allocator>
687 mat_lu.Solve(SeldonNoTrans, x);
693 template<
class T,
class Prop,
class Allocator>
697 mat_lu.Solve(TransA, x);
702 template<
class Allocator>
708 for (
int i = 0; i < x.GetM(); i++)
710 y(i, 0) = real(x(i));
711 y(i, 1) = imag(x(i));
716 for (
int i = 0; i < x.GetM(); i++)
717 x(i) = complex<double>(y(i, 0), y(i, 1));
722 template<
class Allocator>
729 for (
int i = 0; i < x.GetM(); i++)
731 y(i, 0) = real(x(i));
732 y(i, 1) = imag(x(i));
735 SolveLU(TransA, mat_lu, y);
737 for (
int i = 0; i < x.GetM(); i++)
738 x(i) = complex<double>(y(i, 0), y(i, 1));
744 template<
class Allocator>
748 throw WrongArgument(
"SolveLU(MatrixSuperLU<complex<double> >, Vector<double>)",
749 "The result should be a complex vector");
754 template<
class Allocator>
759 throw WrongArgument(
"SolveLU(MatrixSuperLU<complex<double> >, Vector<double>)",
760 "The result should be a complex vector");
765 #define SELDON_FILE_WSMP_CXX