Wsmp.hxx
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_HXX
20 
21 extern "C"
22 {
23  void wsetmaxthrds_(int*);
24  void wsmp_initialize_();
25  void wsmp_clear_();
26 
27 #ifndef SELDON_WITH_MPI
28  void wssmp_(int*, int*, int*, double*, double*, int*, int*,
29  double*, int*, int*, double*, int*, int*, int*, double*);
30 
31  void zssmp_(int*, int*, int*, void*, void*, int*, int*,
32  void*, int*, int*, void*, int*, int*, int*, double*);
33 
34  void wgsmp_(int*, int*, int*, double*, double*, int*, int*,
35  double*, int*, double*);
36 
37  void zgsmp_(int*, int*, int*, void*, void*, int*, int*,
38  void*, int*, void*);
39 #else
40  void pwssmp_(int*, int*, int*, double*, double*, int*, int*,
41  double*, int*, int*, double*, int*, int*, int*, double*);
42 
43  void pzssmp_(int*, int*, int*, void*, void*, int*, int*,
44  void*, int*, int*, void*, int*, int*, int*, double*);
45 
46  void pwgsmp_(int*, int*, int*, double*, double*, int*, int*,
47  double*, int*, double*);
48 
49  void pzgsmp_(int*, int*, int*, void*, void*, int*, int*,
50  void*, int*, void*);
51 #endif
52 }
53 
54 namespace Seldon
55 {
56 
57  template<class T>
58  class MatrixWsmp : public VirtualSparseDirectSolver<T>
59  {
60  protected:
61  int n;
62  IVect iparm; Vector<double> dparm;
63  bool refine_solution, use_pivoting;
64  bool cholesky, symmetric, distributed;
65  IVect Ptr, IndRow; Vector<T> Val;
66  IVect permut, inverse_permut;
67  double threshold_pivot;
68 
69  void CallWssmp(int*, int*, int*, T*, T*, int*, int*,
70  T*, int*, int*, T*, int*, int*, int*, double*);
71 
72  void CallWgsmp(int*, int*, int*, T*, T*, int*, int*,
73  T*, int*, double*);
74 
75  public:
76  MatrixWsmp();
77  ~MatrixWsmp();
78 
79  bool UseInteger8() const;
80  void Clear();
81 
82  void ShowMessages();
83  void HideMessages();
84  void ShowFullHistory();
85 
86  void SetPivotThreshold(double);
87  void SetCholeskyFacto(bool chol);
88  int GetInfoFactorization() const;
89 
90  void RefineSolution();
91  void DoNotRefineSolution();
92  size_t GetMemorySize() const;
93 
94  void SetNumberOfThreadPerNode(int nb_threads);
95 
96  template<class Storage, class Allocator>
97  void FactorizeMatrix(Matrix<T, Symmetric, Storage, Allocator>& mat,
98  bool keep_matrix = false);
99 
100  template<class Storage, class Allocator>
101  void FactorizeMatrix(Matrix<T, General, Storage, Allocator>& mat,
102  bool keep_matrix = false);
103 
104  void FactorizeCSR(Vector<int>& Ptr, Vector<int>& IndRow,
105  Vector<T>& Val, bool sym);
106 
107  void FactorizeUnsymmetric();
108  void FactorizeSymmetric();
109 
110  void Solve(Vector<T>& b);
111  void Solve(const SeldonTranspose& trans, Vector<T>& b);
112  void Solve(const SeldonTranspose&, T* x, int nrhs);
113  void Solve(const SeldonTranspose& trans, Matrix<T, General, ColMajor>& b);
114 
115 #ifdef SELDON_WITH_MPI
116  void FactorizeDistributedMatrix(MPI_Comm& comm_facto, Vector<long>& Ptr,
117  Vector<int>& IndRow, Vector<T>& Val,
118  const Vector<int>& glob_number,
119  bool sym, bool keep_matrix = false);
120 
121  void FactorizeDistributedMatrix(MPI_Comm& comm_facto, Vector<int64_t>& Ptr,
122  Vector<int64_t>& IndRow, Vector<T>& Val,
123  const Vector<int>& glob_number,
124  bool sym, bool keep_matrix = false);
125 
126  template<class Tint>
127  void FactorizeParallel(MPI_Comm& comm_facto,
128  Vector<Tint>&, Vector<Tint>&,
129  Vector<T>&,
130  const Vector<int>& glob_num,
131  bool sym, bool keep_matrix = false);
132 
133  void FactorizeParallel(MPI_Comm& comm_facto,
134  Vector<int>&, Vector<int>&,
135  Vector<T>&, const Vector<int>& glob_number,
136  bool sym, bool keep_matrix = false);
137 
138  void SolveDistributed(MPI_Comm& comm_facto,
139  const SeldonTranspose& TransA,
140  Vector<T>& x, const Vector<int>& glob_num);
141 
142  void SolveDistributed(MPI_Comm& comm_facto,
143  const SeldonTranspose& TransA,
144  Matrix<T, General, ColMajor>& x, const Vector<int>& glob_num);
145 
146  void SolveDistributed(MPI_Comm& comm_facto,
147  const SeldonTranspose& TransA,
148  T* x_ptr, int nrhs, const Vector<int>& glob_num);
149 #endif
150 
151  };
152 
153 
154  template<class T, class Prop, class Storage, class Allocator>
155  void GetLU(Matrix<T, Prop, Storage, Allocator>& A, MatrixWsmp<T>& mat_lu,
156  bool keep_matrix = false);
157 
158  template<class T, class Allocator>
159  void SolveLU(MatrixWsmp<T>& mat_lu, Vector<T, VectFull, Allocator>& x);
160 
161  template<class T, class Allocator>
162  void SolveLU(const SeldonTranspose& TransA,
163  MatrixWsmp<T>& mat_lu, Vector<T, VectFull, Allocator>& x);
164 
165  template<class T, class Prop, class Allocator>
166  void SolveLU(MatrixWsmp<T>& mat_lu,
167  Matrix<T, Prop, ColMajor, Allocator>& x);
168 
169  template<class T, class Prop, class Allocator>
170  void SolveLU(const SeldonTranspose& TransA,
171  MatrixWsmp<T>& mat_lu, Matrix<T, Prop, ColMajor, Allocator>& x);
172 
173  template<class Allocator>
174  void SolveLU(MatrixWsmp<double>& mat_lu,
175  Vector<complex<double>, VectFull, Allocator>& x);
176 
177  template<class Allocator>
178  void SolveLU(const SeldonTranspose& TransA,
179  MatrixWsmp<double>& mat_lu,
180  Vector<complex<double>, VectFull, Allocator>& x);
181 
182  template<class Allocator>
183  void SolveLU(MatrixWsmp<complex<double> >& mat_lu,
184  Vector<double, VectFull, Allocator>& x);
185 
186  template<class Allocator>
187  void SolveLU(const SeldonTranspose& TransA,
188  MatrixWsmp<complex<double> >& mat_lu,
189  Vector<double, VectFull, Allocator>& x);
190 
191 }
192 
193 #define SELDON_FILE_WSMP_HXX
194 #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, VectFull >
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::MatrixWsmp::RefineSolution
void RefineSolution()
Tells to the direct solver that refinement is required.
Definition: Wsmp.cxx:169
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::VirtualSparseDirectSolver
Base class for an interface with a direct solver.
Definition: SparseSolver.hxx:30
Seldon::MatrixWsmp::SetNumberOfThreadPerNode
void SetNumberOfThreadPerNode(int nb_threads)
Sets the number of threads per mpi process.
Definition: Wsmp.cxx:193
Seldon::MatrixWsmp
Definition: Wsmp.hxx:58
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