source: issm/trunk/src/c/toolkits/objects/Matrix.h@ 24686

Last change on this file since 24686 was 24686, checked in by Mathieu Morlighem, 6 years ago

merged trunk-jpl and trunk for revision 24684
File size: 6.0 KB
Line 
1/*!\file: Matrix.h
2 * \brief wrapper to matrix objects. The goal is to control which API (PETSc,Scalpack, Plapack?)
3 * implements our underlying matrix format.
4 */
5
6#ifndef _MATRIX_H_
7#define _MATRIX_H_
8
9/*Headers:*/
10/*{{{*/
11#ifdef HAVE_CONFIG_H
12 #include <config.h>
13#else
14#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
15#endif
16#include <cstring>
17#include "../../shared/Enum/Enum.h"
18#include "../petsc/petscincludes.h"
19#include "../issm/issmtoolkit.h"
20/*}}}*/
21
22enum matrixtype { PetscMatType, IssmMatType };
23
24template <class doubletype> class Vector;
25
26template <class doubletype>
27class Matrix{
28
29 public:
30
31 int type;
32 #ifdef _HAVE_PETSC_
33 PetscMat *pmatrix;
34 #endif
35 IssmMat<doubletype> *imatrix;
36
37 /*Matrix constructors, destructors*/
38 Matrix(){/*{{{*/
39 InitCheckAndSetType();
40 }
41 /*}}}*/
42 Matrix(int M,int N){/*{{{*/
43
44 InitCheckAndSetType();
45
46 if(type==PetscMatType){
47 #ifdef _HAVE_PETSC_
48 this->pmatrix=new PetscMat(M,N);
49 #endif
50 }
51 else{
52 this->imatrix=new IssmMat<doubletype>(M,N);
53 }
54
55 }
56 /*}}}*/
57 Matrix(int m,int n,int M,int N,int* d_nnz,int* o_nnz){/*{{{*/
58
59 InitCheckAndSetType();
60
61 if(type==PetscMatType){
62 #ifdef _HAVE_PETSC_
63 this->pmatrix=new PetscMat(m,n,M,N,d_nnz,o_nnz);
64 #endif
65 }
66 else{
67 this->imatrix=new IssmMat<doubletype>(m,n,M,N,d_nnz,o_nnz);
68 }
69
70 }
71 /*}}}*/
72 Matrix(int M,int N,double sparsity){/*{{{*/
73
74 InitCheckAndSetType();
75
76 if(type==PetscMatType){
77 #ifdef _HAVE_PETSC_
78 this->pmatrix=new PetscMat(M,N,sparsity);
79 #endif
80 }
81 else{
82 this->imatrix=new IssmMat<doubletype>(M,N,sparsity);
83 }
84 }
85 /*}}}*/
86 Matrix(IssmPDouble* serial_mat,int M,int N,IssmPDouble sparsity){/*{{{*/
87
88 InitCheckAndSetType();
89
90 if(type==PetscMatType){
91 #ifdef _HAVE_PETSC_
92 this->pmatrix=new PetscMat(serial_mat,M,N,sparsity);
93 #endif
94 }
95 else{
96 this->imatrix=new IssmMat<doubletype>(serial_mat,M,N,sparsity);
97 }
98
99 }
100 /*}}}*/
101 Matrix(int M,int N,int connectivity,int numberofdofspernode){/*{{{*/
102
103 InitCheckAndSetType();
104
105 if(type==PetscMatType){
106 #ifdef _HAVE_PETSC_
107 this->pmatrix=new PetscMat(M,N,connectivity,numberofdofspernode);
108 #endif
109 }
110 else{
111 this->imatrix=new IssmMat<doubletype>(M,N,connectivity,numberofdofspernode);
112 }
113
114 }
115 /*}}}*/
116 ~Matrix(){/*{{{*/
117
118 if(type==PetscMatType){
119 #ifdef _HAVE_PETSC_
120 delete this->pmatrix;
121 #endif
122 }
123 else delete this->imatrix;
124
125 }
126 /*}}}*/
127 void InitCheckAndSetType(void){/*{{{*/
128
129 #ifdef _HAVE_PETSC_
130 pmatrix=NULL;
131 #endif
132 imatrix=NULL;
133
134 /*retrieve toolkittype: */
135 char* toolkittype=ToolkitOptions::GetToolkitType();
136
137 /*set matrix type: */
138 if (strcmp(toolkittype,"petsc")==0){
139 #ifdef _HAVE_PETSC_
140 type=PetscMatType;
141 #else
142 _error_("cannot create petsc matrix without PETSC compiled!");
143 #endif
144 }
145 else if(strcmp(toolkittype,"issm")==0){
146 /*let this choice stand:*/
147 type=IssmMatType;
148 }
149 else{
150 _error_("unknow toolkit type ");
151 }
152
153 /*Free ressources: */
154 xDelete<char>(toolkittype);
155 }
156 /*}}}*/
157
158 /*Matrix specific routines:*/
159 void Echo(void){/*{{{*/
160 _assert_(this);
161
162 if(type==PetscMatType){
163 #ifdef _HAVE_PETSC_
164 this->pmatrix->Echo();
165 #endif
166 }
167 else{
168 this->imatrix->Echo();
169 }
170
171 }
172 /*}}}*/
173 void AllocationInfo(void){/*{{{*/
174 _assert_(this);
175 if(type==PetscMatType){
176 #ifdef _HAVE_PETSC_
177 this->pmatrix->AllocationInfo();
178 #endif
179 }
180 else{
181 //this->imatrix->AllocationInfo();
182 _error_("not supported yet");
183 }
184 }/*}}}*/
185 void Assemble(void){/*{{{*/
186
187 if(type==PetscMatType){
188 #ifdef _HAVE_PETSC_
189 this->pmatrix->Assemble();
190 #endif
191 }
192 else{
193 this->imatrix->Assemble();
194 }
195 }
196 /*}}}*/
197 IssmDouble Norm(NormMode norm_type){/*{{{*/
198
199 IssmDouble norm=0;
200
201 if(type==PetscMatType){
202 #ifdef _HAVE_PETSC_
203 norm=this->pmatrix->Norm(norm_type);
204 #endif
205 }
206 else{
207 norm=this->imatrix->Norm(norm_type);
208 }
209
210 return norm;
211 }
212 /*}}}*/
213 void GetSize(int* pM,int* pN){/*{{{*/
214
215 if(type==PetscMatType){
216 #ifdef _HAVE_PETSC_
217 this->pmatrix->GetSize(pM,pN);
218 #endif
219 }
220 else{
221 this->imatrix->GetSize(pM,pN);
222 }
223
224 }
225 /*}}}*/
226 void GetLocalSize(int* pM,int* pN){/*{{{*/
227
228 if(type==PetscMatType){
229 #ifdef _HAVE_PETSC_
230 this->pmatrix->GetLocalSize(pM,pN);
231 #endif
232 }
233 else{
234 this->imatrix->GetLocalSize(pM,pN);
235 }
236
237 }
238 /*}}}*/
239 void MatMult(Vector<doubletype>* X,Vector<doubletype>* AX){/*{{{*/
240
241 if(type==PetscMatType){
242 #ifdef _HAVE_PETSC_
243 this->pmatrix->MatMult(X->pvector,AX->pvector);
244 #endif
245 }
246 else{
247 this->imatrix->MatMult(X->ivector,AX->ivector);
248 }
249
250 }
251 /*}}}*/
252 Matrix<doubletype>* Duplicate(void){/*{{{*/
253
254 Matrix<doubletype>* output=new Matrix<doubletype>();
255
256 if(type==PetscMatType){
257 #ifdef _HAVE_PETSC_
258 output->pmatrix=this->pmatrix->Duplicate();
259 #endif
260 }
261 else{
262 output->imatrix=this->imatrix->Duplicate();
263 }
264
265 return output;
266 }
267 /*}}}*/
268 doubletype* ToSerial(void){/*{{{*/
269
270 doubletype* output=NULL;
271
272 if(type==PetscMatType){
273 #ifdef _HAVE_PETSC_
274 output=this->pmatrix->ToSerial();
275 #endif
276 }
277 else{
278 output=this->imatrix->ToSerial();
279 }
280
281 return output;
282 }
283 /*}}}*/
284 void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){/*{{{*/
285
286 if(type==PetscMatType){
287 #ifdef _HAVE_PETSC_
288 this->pmatrix->SetValues(m,idxm,n,idxn,values,mode);
289 #endif
290 }
291 else{
292 this->imatrix->SetValues(m,idxm,n,idxn,values,mode);
293 }
294 }
295 /*}}}*/
296 void Convert(MatrixType newtype){/*{{{*/
297
298 if(type==PetscMatType){
299 #ifdef _HAVE_PETSC_
300 this->pmatrix->Convert(newtype);
301 #endif
302 }
303 else{
304 this->imatrix->Convert(newtype);
305 }
306
307 }
308 /*}}}*/
309 /*
310 * sets all values to 0 but keeps the structure of a sparse matrix
311 */
312 void SetZero(void) {/*{{{*/
313 if(type==PetscMatType){
314 #ifdef _HAVE_PETSC_
315 this->pmatrix->SetZero();
316 #endif
317 }
318 else{
319 this->imatrix->SetZero();
320 }
321 }
322 /*}}}*/
323};
324
325#endif //#ifndef _MATRIX_H_
Note: See TracBrowser for help on using the repository browser.