source: issm/oecreview/Archive/14312-15392/ISSM-14791-14792.diff

Last change on this file was 15393, checked in by Mathieu Morlighem, 12 years ago

NEW: adding Archive/14312-15392 for oecreview
File size: 81.7 KB
RevLine 
[15393]1Index: ../trunk-jpl/src/c/modules/Solverx/SolverxPetsc.cpp
2===================================================================
3--- ../trunk-jpl/src/c/modules/Solverx/SolverxPetsc.cpp (revision 14791)
4+++ ../trunk-jpl/src/c/modules/Solverx/SolverxPetsc.cpp (revision 14792)
5@@ -1,159 +0,0 @@
6-/*!\file SolverxPetsc
7- * \brief Petsc implementation of solver
8- */
9-
10-#include "./Solverx.h"
11-#include "../../shared/shared.h"
12-#include "../../include/include.h"
13-#include "../../io/io.h"
14-
15-#ifdef HAVE_CONFIG_H
16- #include <config.h>
17-#else
18-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
19-#endif
20-
21-void SolverxPetsc(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters){
22-
23- PetscVec* uf=new PetscVec();
24-
25- Vec uf0_vector = NULL;
26- Vec df_vector = NULL;
27-
28- if(uf0) uf0_vector = uf0->vector;
29- if(df) df_vector = df->vector;
30-
31- SolverxPetsc(&uf->vector, Kff->matrix, pf->vector, uf0_vector, df_vector, parameters);
32-
33- *puf=uf;
34-
35-}
36-void SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){
37-
38- /*Output: */
39- Vec uf = NULL;
40-
41- /*Intermediary: */
42- int local_m,local_n,global_m,global_n;
43-
44- /*Solver */
45- KSP ksp = NULL;
46- PC pc = NULL;
47- int iteration_number;
48- int solver_type;
49- bool fromlocalsize = true;
50- #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
51- PetscTruth flag,flg;
52- #else
53- PetscBool flag,flg;
54- #endif
55-
56- /*Stokes: */
57- IS isv=NULL;
58- IS isp=NULL;
59-
60- #if _PETSC_MAJOR_ >= 3
61- char ksp_type[50];
62- #endif
63-
64- /*Display message*/
65- if(VerboseModule()) _pprintLine_(" Solving");
66- #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
67- if(VerboseSolver())PetscOptionsPrint(stdout);
68- #else
69- if(VerboseSolver())PetscOptionsView(PETSC_VIEWER_STDOUT_WORLD);
70- #endif
71-
72- /*First, check that f-set is not NULL, i.e. model is fully constrained:*/
73- _assert_(Kff);
74- MatGetSize(Kff,&global_m,&global_n); _assert_(global_m==global_n);
75- if(!global_n){
76- *puf=NewVec(0,IssmComm::GetComm()); return;
77- }
78-
79- /*Initial guess */
80- /*Now, check that we are not giving an initial guess to the solver, if we are running a direct solver: */
81- #if _PETSC_MAJOR_ >= 3
82- PetscOptionsGetString(PETSC_NULL,"-ksp_type",ksp_type,49,&flg);
83- if (strcmp(ksp_type,"preonly")==0)uf0=NULL;
84- #endif
85-
86- /*If initial guess for the solution exists, use it to create uf, otherwise,
87- * duplicate the right hand side so that the solution vector has the same structure*/
88- if(uf0){
89- VecDuplicate(uf0,&uf); VecCopy(uf0,uf);
90- }
91- else{
92- MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,IssmComm::GetComm(),fromlocalsize);
93- }
94-
95- /*Process petsc options to see if we are using special types of external solvers*/
96- PetscOptionsDetermineSolverType(&solver_type);
97-
98- /*Check the solver is available*/
99- if(solver_type==MUMPSPACKAGE_LU || solver_type==MUMPSPACKAGE_CHOL){
100- #if _PETSC_MAJOR_ >=3
101- #ifndef _HAVE_MUMPS_
102- _error_("requested MUMPS solver, which was not compiled into ISSM!\n");
103- #endif
104- #endif
105- }
106-
107- /*Prepare solver*/
108- KSPCreate(IssmComm::GetComm(),&ksp);
109- KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN);
110- KSPSetFromOptions(ksp);
111-
112- #if _PETSC_MAJOR_==3
113- /*Specific solver?: */
114- KSPGetPC(ksp,&pc);
115- if (solver_type==MUMPSPACKAGE_LU){
116- #if _PETSC_MINOR_==1
117- PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS);
118- #else
119- PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
120- #endif
121- }
122-
123- /*Stokes: */
124- if (solver_type==StokesSolverEnum){
125- /*Make indices out of doftypes: */
126- if(!df)_error_("need doftypes for Stokes solver!\n");
127- DofTypesToIndexSet(&isv,&isp,df,StokesSolverEnum);
128-
129- /*Set field splits: */
130- KSPGetPC(ksp,&pc);
131- #if _PETSC_MINOR_==1
132- PCFieldSplitSetIS(pc,isv);
133- PCFieldSplitSetIS(pc,isp);
134- #else
135- PCFieldSplitSetIS(pc,PETSC_NULL,isv);
136- PCFieldSplitSetIS(pc,PETSC_NULL,isp);
137- #endif
138-
139- }
140- #endif
141-
142- /*If there is an initial guess for the solution, use it
143- * except if we are using the MUMPS direct solver
144- * where any initial guess will crash Petsc*/
145- if (uf0){
146- if((solver_type!=MUMPSPACKAGE_LU) && (solver_type!=MUMPSPACKAGE_CHOL) && (solver_type!=SPOOLESPACKAGE_LU) && (solver_type!=SPOOLESPACKAGE_CHOL) && (solver_type!=SUPERLUDISTPACKAGE)){
147- KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
148- }
149- }
150-
151- /*Solve: */
152- if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
153- KSPSolve(ksp,pf,uf);
154-
155- /*Check convergence*/
156- KSPGetIterationNumber(ksp,&iteration_number);
157- if (iteration_number<0) _error_("Solver diverged at iteration number: " << -iteration_number);
158-
159- /*Free resources:*/
160- KSPFree(&ksp);
161-
162- /*Assign output pointers:*/
163- *puf=uf;
164-}
165Index: ../trunk-jpl/src/c/modules/Solverx/DofTypesToIndexSet.cpp
166===================================================================
167--- ../trunk-jpl/src/c/modules/Solverx/DofTypesToIndexSet.cpp (revision 14791)
168+++ ../trunk-jpl/src/c/modules/Solverx/DofTypesToIndexSet.cpp (revision 14792)
169@@ -1,82 +0,0 @@
170-/*!\file Solverx
171- * \brief solver
172- */
173-
174-#include "./Solverx.h"
175-#include "../../shared/shared.h"
176-#include "../../include/include.h"
177-
178-#ifdef HAVE_CONFIG_H
179- #include <config.h>
180-#else
181-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
182-#endif
183-
184-void DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum){
185-
186- /*output: */
187- IS isv=NULL;
188- IS isp=NULL;
189-
190- int start,end;
191- IssmDouble* df_local=NULL;
192- int df_local_size;
193- int i;
194-
195- int* pressure_indices=NULL;
196- int* velocity_indices=NULL;
197- int pressure_num=0;
198- int velocity_num=0;
199- int pressure_count=0;
200- int velocity_count=0;
201-
202- if(typeenum==StokesSolverEnum){
203-
204- /*Ok, recover doftypes vector values and indices: */
205- VecGetOwnershipRange(df,&start,&end);
206- VecGetLocalSize(df,&df_local_size);
207- VecGetArray(df,&df_local);
208-
209- pressure_num=0;
210- velocity_num=0;
211- for(i=0;i<df_local_size;i++){
212- if (df_local[i]==PressureEnum)pressure_num++;
213- else velocity_num++;
214- }
215-
216- /*Allocate indices: */
217- if(pressure_num)pressure_indices=xNew<int>(pressure_num);
218- if(velocity_num)velocity_indices=xNew<int>(velocity_num);
219-
220- pressure_count=0;
221- velocity_count=0;
222- for(i=0;i<df_local_size;i++){
223- if (df_local[i]==PressureEnum){
224- pressure_indices[pressure_count]=start+i;
225- pressure_count++;
226- }
227- if (df_local[i]==VelocityEnum){
228- velocity_indices[velocity_count]=start+i;
229- velocity_count++;
230- }
231- }
232- VecRestoreArray(df,&df_local);
233-
234- /*Create indices sets: */
235- #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
236- ISCreateGeneral(IssmComm::GetComm(),pressure_num,pressure_indices,&isp);
237- ISCreateGeneral(IssmComm::GetComm(),velocity_num,velocity_indices,&isv);
238- #else
239- ISCreateGeneral(IssmComm::GetComm(),pressure_num,pressure_indices,PETSC_COPY_VALUES,&isp);
240- ISCreateGeneral(IssmComm::GetComm(),velocity_num,velocity_indices,PETSC_COPY_VALUES,&isv);
241- #endif
242- }
243-
244- /*Free ressources:*/
245- xDelete<int>(pressure_indices);
246- xDelete<int>(velocity_indices);
247-
248- /*Assign output pointers:*/
249- *pisv=isv;
250- *pisp=isp;
251-}
252Index: ../trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp
253===================================================================
254--- ../trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp (revision 14791)
255+++ ../trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp (revision 14792)
256@@ -1,255 +0,0 @@
257-/*!\file SolverxSeq
258- * \brief implementation of sequential solver using the GSL librarie
259- */
260-
261-#ifdef HAVE_CONFIG_H
262- #include <config.h>
263-#else
264-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
265-#endif
266-#include <cstring>
267-
268-#include "./Solverx.h"
269-#include "../../shared/shared.h"
270-#include "../../classes/classes.h"
271-#include "../../include/include.h"
272-#include "../../io/io.h"
273-
274-#ifdef _HAVE_GSL_
275-#include <gsl/gsl_linalg.h>
276-#endif
277-
278-void SolverxSeq(IssmVec<IssmDouble>** pout,IssmMat<IssmDouble>* Kffin, IssmVec<IssmDouble>* pfin, Parameters* parameters){/*{{{*/
279-
280- /*First off, we assume that the type of IssmVec is IssmSeqVec and the type of IssmMat is IssmDenseMat. So downcast: */
281- IssmSeqVec<IssmDouble>* pf=(IssmSeqVec<IssmDouble>*)pfin->vector;
282- IssmDenseMat<IssmDouble>* Kff=(IssmDenseMat<IssmDouble>*)Kffin->matrix;
283-
284-#ifdef _HAVE_GSL_
285- /*Intermediary: */
286- int M,N,N2;
287- IssmSeqVec<IssmDouble> *uf = NULL;
288-
289- Kff->GetSize(&M,&N);
290- pf->GetSize(&N2);
291-
292- if(N!=N2)_error_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !");
293- if(M!=N)_error_("Stiffness matrix should be square!");
294-#ifdef _HAVE_ADOLC_
295- ensureContiguousLocations(N);
296-#endif
297- IssmDouble *x = xNew<IssmDouble>(N);
298-#ifdef _HAVE_ADOLC_
299- SolverxSeq(x,Kff->matrix,pf->vector,N,parameters);
300-#else
301- SolverxSeq(x,Kff->matrix,pf->vector,N);
302-#endif
303-
304- uf=new IssmSeqVec<IssmDouble>(x,N);
305- xDelete(x);
306-
307- /*Assign output pointers:*/
308- IssmVec<IssmDouble>* out=new IssmVec<IssmDouble>();
309- out->vector=(IssmAbsVec<IssmDouble>*)uf;
310- *pout=out;
311-
312-#else
313- _error_("GSL support not compiled in!");
314-#endif
315-
316-}/*}}}*/
317-void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
318-
319- /*Allocate output*/
320- double* X = xNew<double>(n);
321-
322- /*Solve*/
323- SolverxSeq(X,A,B,n);
324-
325- /*Assign output pointer*/
326- *pX=X;
327-}
328-/*}}}*/
329-
330-#ifdef _HAVE_ADOLC_
331-int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y){ /*{{{*/
332- SolverxSeq(y,x, x+m*m, m); // x is where the matrix starts, x+m*m is where the right-hand side starts
333- return 0;
334-} /*}}}*/
335-int EDF_fos_forward_for_solverx(int n, IssmPDouble *inVal, IssmPDouble *inDeriv, int m, IssmPDouble *outVal, IssmPDouble *outDeriv) { /*{{{*/
336-#ifdef _HAVE_GSL_
337- // for (int i=0; i<m*m; ++i) std::cout << "EDF_fos_forward_for_solverx A["<< i << "]=" << inVal[i] << std::endl;
338- // for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx b["<< i << "]=" << inVal[i+m*m] << std::endl;
339- // the matrix will be modified by LU decomposition. Use gsl_A copy
340- double* Acopy = xNew<double>(m*m);
341- xMemCpy(Acopy,inVal,m*m);
342- /*Initialize gsl matrices and vectors: */
343- gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
344- gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
345- gsl_permutation *perm_p = gsl_permutation_alloc (m);
346- int signPerm;
347- // factorize
348- gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
349- gsl_vector *gsl_x_p = gsl_vector_alloc (m);
350- // solve for the value
351- gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
352- /*Copy result*/
353- xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
354- gsl_vector_free(gsl_x_p);
355- // for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx x["<< i << "]=" << outVal[i] << std::endl;
356- // solve for the derivatives acc. to A * dx = r with r=db - dA * x
357- // compute the RHS
358- double* r=xNew<double>(m);
359- for (int i=0; i<m; i++) {
360- r[i]=inDeriv[m*m+i]; // this is db[i]
361- for (int j=0;j<m; j++) {
362- r[i]-=inDeriv[i*m+j]*outVal[j]; // this is dA[i][j]*x[j]
363- }
364- }
365- gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
366- gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
367- gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
368- xMemCpy(outDeriv,gsl_vector_ptr(gsl_dx_p,0),m);
369- gsl_vector_free(gsl_dx_p);
370- xDelete(r);
371- gsl_permutation_free(perm_p);
372- xDelete(Acopy);
373-#endif
374- return 0;
375-} /*}}}*/
376-int EDF_fov_forward_for_solverx(int n, IssmPDouble *inVal, int directionCount, IssmPDouble **inDeriv, int m, IssmPDouble *outVal, IssmPDouble **outDeriv) { /*{{{*/
377-#ifdef _HAVE_GSL_
378- // the matrix will be modified by LU decomposition. Use gsl_A copy
379- double* Acopy = xNew<double>(m*m);
380- xMemCpy(Acopy,inVal,m*m);
381- /*Initialize gsl matrices and vectors: */
382- gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
383- gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
384- gsl_permutation *perm_p = gsl_permutation_alloc (m);
385- int signPerm;
386- // factorize
387- gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
388- gsl_vector *gsl_x_p = gsl_vector_alloc (m);
389- // solve for the value
390- gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
391- /*Copy result*/
392- xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
393- gsl_vector_free(gsl_x_p);
394- // solve for the derivatives acc. to A * dx = r with r=db - dA * x
395- double* r=xNew<double>(m);
396- gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
397- for (int dir=0;dir<directionCount;++dir) {
398- // compute the RHS
399- for (int i=0; i<m; i++) {
400- r[i]=inDeriv[m*m+i][dir]; // this is db[i]
401- for (int j=0;j<m; j++) {
402- r[i]-=inDeriv[i*m+j][dir]*outVal[j]; // this is dA[i][j]*x[j]
403- }
404- }
405- gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
406- gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
407- // reuse r
408- xMemCpy(r,gsl_vector_ptr(gsl_dx_p,0),m);
409- for (int i=0; i<m; i++) {
410- outDeriv[i][dir]=r[i];
411- }
412- }
413- gsl_vector_free(gsl_dx_p);
414- xDelete(r);
415- gsl_permutation_free(perm_p);
416- xDelete(Acopy);
417-#endif
418- return 0;
419-}
420-/*}}}*/
421-int EDF_fos_reverse_for_solverx(int m, double *dp_U, int n, double *dp_Z, double* dp_x, double* dp_y) { /*{{{*/
422- // copy to transpose the matrix
423- double* transposed=xNew<double>(m*m);
424- for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j];
425- gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m);
426- // the adjoint of the solution is our right-hand side
427- gsl_vector_view x_bar=gsl_vector_view_array(dp_U,m);
428- // the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute:
429- gsl_vector_view b_bar=gsl_vector_view_array(dp_Z+m*m,m);
430- gsl_permutation *perm_p = gsl_permutation_alloc (m);
431- int permSign;
432- gsl_linalg_LU_decomp (&aTransposed.matrix, perm_p, &permSign);
433- gsl_linalg_LU_solve (&aTransposed.matrix, perm_p, &x_bar.vector, &b_bar.vector);
434- // now do the adjoint of the matrix
435- for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) dp_Z[i*m+j]-=dp_Z[m*m+i]*dp_y[j];
436- gsl_permutation_free(perm_p);
437- xDelete(transposed);
438- return 0;
439-}
440-/*}}}*/
441-int EDF_fov_reverse_for_solverx(int m, int p, double **dpp_U, int n, double **dpp_Z, double* dp_x, double* dp_y) { /*{{{*/
442- // copy to transpose the matrix
443- double* transposed=xNew<double>(m*m);
444- for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j];
445- gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m);
446- gsl_permutation *perm_p = gsl_permutation_alloc (m);
447- int permSign;
448- gsl_linalg_LU_decomp (&aTransposed.matrix, perm_p, &permSign);
449- for (int weightsRowIndex=0;weightsRowIndex<p;++weightsRowIndex) {
450- // the adjoint of the solution is our right-hand side
451- gsl_vector_view x_bar=gsl_vector_view_array(dpp_U[weightsRowIndex],m);
452- // the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute:
453- gsl_vector_view b_bar=gsl_vector_view_array(dpp_Z[weightsRowIndex]+m*m,m);
454- gsl_linalg_LU_solve (&aTransposed.matrix, perm_p, &x_bar.vector, &b_bar.vector);
455- // now do the adjoint of the matrix
456- for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) dpp_Z[weightsRowIndex][i*m+j]-=dpp_Z[weightsRowIndex][m*m+i]*dp_y[j];
457- }
458- gsl_permutation_free(perm_p);
459- xDelete(transposed);
460- return 0;
461-}
462-/*}}}*/
463-void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
464- // pack inputs to conform to the EDF-prescribed interface
465- // ensure a contiguous block of locations:
466- ensureContiguousLocations(n*(n+1));
467- IssmDouble* adoubleEDFin=xNew<IssmDouble>(n*(n+1)); // packed inputs, i.e. matrix and right hand side
468- for(int i=0; i<n*n;i++)adoubleEDFin[i] =A[i]; // pack matrix
469- for(int i=0; i<n; i++)adoubleEDFin[i+n*n]=B[i]; // pack the right hand side
470- IssmPDouble* pdoubleEDFin=xNew<IssmPDouble>(n*(n+1)); // provide space to transfer inputs during call_ext_fct
471- IssmPDouble* pdoubleEDFout=xNew<IssmPDouble>(n); // provide space to transfer outputs during call_ext_fct
472- // call the wrapped solver through the registry entry we retrieve from parameters
473- call_ext_fct(dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p,
474- n*(n+1), pdoubleEDFin, adoubleEDFin,
475- n, pdoubleEDFout,X);
476- // for(int i=0; i<n; i++) {ADOLC_DUMP_MACRO(X[i]);}
477- xDelete(adoubleEDFin);
478- xDelete(pdoubleEDFin);
479- xDelete(pdoubleEDFout);
480-}
481-/*}}}*/
482-#endif
483-void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
484-#ifdef _HAVE_GSL_
485- /*GSL Matrices and vectors: */
486- int s;
487- gsl_matrix_view a;
488- gsl_vector_view b,x;
489- gsl_permutation *p = NULL;
490-// for (int i=0; i<n*n; ++i) std::cout << "SolverxSeq A["<< i << "]=" << A[i] << std::endl;
491-// for (int i=0; i<n; ++i) std::cout << "SolverxSeq b["<< i << "]=" << B[i] << std::endl;
492- /*A will be modified by LU decomposition. Use copy*/
493- double* Acopy = xNew<double>(n*n);
494- xMemCpy(Acopy,A,n*n);
495-
496- /*Initialize gsl matrices and vectors: */
497- a = gsl_matrix_view_array (Acopy,n,n);
498- b = gsl_vector_view_array (B,n);
499- x = gsl_vector_view_array (X,n);
500-
501- /*Run LU and solve: */
502- p = gsl_permutation_alloc (n);
503- gsl_linalg_LU_decomp (&a.matrix, p, &s);
504- gsl_linalg_LU_solve (&a.matrix, p, &b.vector, &x.vector);
505-
506- /*Clean up and assign output pointer*/
507- xDelete(Acopy);
508- gsl_permutation_free(p);
509-#endif
510-}
511-/*}}}*/
512Index: ../trunk-jpl/src/c/modules/Solverx/Solverx.cpp
513===================================================================
514--- ../trunk-jpl/src/c/modules/Solverx/Solverx.cpp (revision 14791)
515+++ ../trunk-jpl/src/c/modules/Solverx/Solverx.cpp (revision 14792)
516@@ -2,48 +2,30 @@
517 * \brief solver
518 */
519
520-#include "./Solverx.h"
521-#include "../../shared/shared.h"
522-#include "../../include/include.h"
523-#include "../../io/io.h"
524-
525 #ifdef HAVE_CONFIG_H
526 #include <config.h>
527 #else
528 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
529 #endif
530
531+#include "./Solverx.h"
532+#include "../../shared/shared.h"
533+
534 void Solverx(Vector<IssmDouble>** puf, Matrix<IssmDouble>* Kff, Vector<IssmDouble>* pf, Vector<IssmDouble>* uf0,Vector<IssmDouble>* df, Parameters* parameters){
535
536+ /*intermediary: */
537+ Solver<IssmDouble> *solver=NULL;
538+
539 /*output: */
540 Vector<IssmDouble> *uf=NULL;
541
542- /*In debugging mode, check that stiffness matrix and load vectors are not NULL (they can be empty)*/
543- _assert_(Kff);
544- _assert_(pf);
545-
546 if(VerboseModule()) _pprintLine_(" Solving matrix system");
547
548- /*Initialize vector: */
549- uf=new Vector<IssmDouble>();
550+ /*Initialize solver: */
551+ solver=new Solver<IssmDouble>(Kff,pf,uf0,df,parameters);
552
553- /*According to matrix type, use specific solvers: */
554- switch(Kff->type){
555- #ifdef _HAVE_PETSC_
556- case PetscMatType:{
557- PetscVec* uf0_vector = NULL;
558- PetscVec* df_vector = NULL;
559- if(uf0) uf0_vector = uf0->pvector;
560- if(df) df_vector = df->pvector;
561- SolverxPetsc(&uf->pvector,Kff->pmatrix,pf->pvector,uf0_vector,df_vector,parameters);
562- break;}
563- #endif
564- case IssmMatType:{
565- SolverxSeq(&uf->ivector,Kff->imatrix,pf->ivector,parameters);
566- break;}
567- default:
568- _error_("Matrix type: " << Kff->type << " not supported yet!");
569- }
570+ /*Solve:*/
571+ uf=solver->Solve();
572
573 /*Assign output pointers:*/
574 *puf=uf;
575Index: ../trunk-jpl/src/c/modules/Solverx/Solverx.h
576===================================================================
577--- ../trunk-jpl/src/c/modules/Solverx/Solverx.h (revision 14791)
578+++ ../trunk-jpl/src/c/modules/Solverx/Solverx.h (revision 14792)
579@@ -11,30 +11,9 @@
580 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
581 #endif
582
583-#include "../../classes/objects/objects.h"
584-#include "../../toolkits/toolkits.h"
585+#include "../../classes/toolkits/toolkitobjects.h"
586
587 /* local prototypes: */
588 void Solverx(Vector<IssmDouble>** puf, Matrix<IssmDouble>* Kff, Vector<IssmDouble>* pf, Vector<IssmDouble>* uf0,Vector<IssmDouble>* df, Parameters* parameters);
589
590-#ifdef _HAVE_PETSC_
591-void SolverxPetsc(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters);
592-void SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters);
593-void DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum);
594-#endif
595-
596-void SolverxSeq(IssmVec<IssmDouble>** puf,IssmMat<IssmDouble>* Kff, IssmVec<IssmDouble>* pf,Parameters* parameters);
597-void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n);
598-void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n);
599-
600-#if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_)
601-void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters);
602-// call back functions:
603-ADOLC_ext_fct EDF_for_solverx;
604-ADOLC_ext_fct_fos_forward EDF_fos_forward_for_solverx;
605-ADOLC_ext_fct_fos_reverse EDF_fos_reverse_for_solverx;
606-ADOLC_ext_fct_fov_forward EDF_fov_forward_for_solverx;
607-ADOLC_ext_fct_fov_reverse EDF_fov_reverse_for_solverx;
608-#endif
609-
610 #endif /* _SOLVERX_H */
611Index: ../trunk-jpl/src/c/modules/Solverx/CMakeLists.txt
612===================================================================
613--- ../trunk-jpl/src/c/modules/Solverx/CMakeLists.txt (revision 14791)
614+++ ../trunk-jpl/src/c/modules/Solverx/CMakeLists.txt (revision 14792)
615@@ -4,6 +4,5 @@
616 include_directories(AFTER $ENV{ISSM_DIR}/src/c/modules/Solverx)
617 # }}}
618 # CORE_SOURCES {{{
619-set(CORE_SOURCES $ENV{ISSM_DIR}/src/c/modules/Solverx/Solverx.cpp
620- $ENV{ISSM_DIR}/src/c/modules/Solverx/SolverxSeq.cpp PARENT_SCOPE)
621+set(CORE_SOURCES $ENV{ISSM_DIR}/src/c/modules/Solverx/Solverx.cpp PARENT_SCOPE)
622 # }}}
623Index: ../trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h
624===================================================================
625--- ../trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h (revision 14791)
626+++ ../trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h (revision 14792)
627@@ -21,6 +21,7 @@
628 #include "../../shared/MemOps/xMemCpy.h"
629 #include "../../shared/Alloc/alloc.h"
630 #include "../../include/macros.h"
631+#include "../../io/io.h"
632 #ifdef _HAVE_MPI_
633 #include "../mpi/mpiincludes.h"
634 #endif
635@@ -302,22 +303,77 @@
636 /*}}}*/
637 /*FUNCTION AXPY{{{*/
638 void AXPY(IssmAbsVec<doubletype>* Xin, doubletype a){
639- printf("AXPY not implemented yet!");
640- exit(1);
641+
642+ int i;
643+
644+ /*Assume X is of the correct type, and downcast: */
645+ IssmMpiVec* X=NULL;
646+
647+ X=(IssmMpiVec<doubletype>*)Xin;
648+
649+ /*y=a*x+y where this->vector is y*/
650+ for(i=0;i<this->m;i++)this->vector[i]=a*X->vector[i]+this->vector[i];
651+
652+
653 }
654 /*}}}*/
655 /*FUNCTION AYPX{{{*/
656 void AYPX(IssmAbsVec<doubletype>* Xin, doubletype a){
657- printf("AYPX not implemented yet!");
658- exit(1);
659+ int i;
660+
661+ /*Assume X is of the correct type, and downcast: */
662+ IssmMpiVec* X=NULL;
663+
664+ X=(IssmMpiVec<doubletype>*)Xin;
665+
666+ /*y=x+a*y where this->vector is y*/
667+ for(i=0;i<this->m;i++)this->vector[i]=X->vector[i]+a*this->vector[i];
668+
669 }
670 /*}}}*/
671 /*FUNCTION ToMPISerial{{{*/
672 doubletype* ToMPISerial(void){
673+
674+ /*communicator info: */
675+ MPI_Comm comm;
676+ int num_procs;
677+
678+ /*MPI_Allgatherv info: */
679+ int lower_row,upper_row;
680+ int* recvcounts=NULL;
681+ int* displs=NULL;
682
683- printf("IssmMpiVec ToMPISerial not implemented yet!");
684- exit(1);
685+ /*output: */
686+ doubletype* buffer=NULL;
687
688+ /*initialize comm info: */
689+ comm=IssmComm::GetComm();
690+ num_procs=IssmComm::GetSize();
691+
692+ /*Allocate: */
693+ buffer=xNew<doubletype>(M);
694+ recvcounts=xNew<int>(num_procs);
695+ displs=xNew<int>(num_procs);
696+
697+ /*recvcounts:*/
698+ MPI_Allgather(&this->m,1,MPI_INT,recvcounts,1,MPI_INT,comm);
699+
700+ /*get lower_row: */
701+ GetOwnershipBoundariesFromRange(&lower_row,&upper_row,this->m,comm);
702+
703+ /*displs: */
704+ MPI_Allgather(&lower_row,1,MPI_INT,displs,1,MPI_INT,comm);
705+
706+ /*All gather:*/
707+ MPI_Allgatherv(this->vector, this->m, MPI_DOUBLE, buffer, recvcounts, displs, MPI_DOUBLE,comm);
708+
709+ /*free ressources: */
710+ xDelete<int>(recvcounts);
711+ xDelete<int>(displs);
712+
713+ /*return: */
714+ return buffer;
715+
716 }
717 /*}}}*/
718 /*FUNCTION Copy{{{*/
719@@ -337,9 +393,28 @@
720 /*}}}*/
721 /*FUNCTION Norm{{{*/
722 doubletype Norm(NormMode mode){
723+
724+ doubletype local_norm;
725+ doubletype norm;
726+ int i;
727
728- printf("IssmMpiVec Norm not implemented yet!");
729- exit(1);
730+ switch(mode){
731+ case NORM_INF:
732+ //local_norm=0; for(i=0;i<this->m;i++)local_norm=max(local_norm,fabs(this->vector[i]));
733+ local_norm=0; for(i=0;i<this->m;i++)local_norm=max(local_norm,this->vector[i]);
734+ MPI_Reduce(&local_norm, &norm, 1, MPI_DOUBLE, MPI_MAX, 0, IssmComm::GetComm());
735+ return norm;
736+ break;
737+ case NORM_TWO:
738+ local_norm=0;
739+ for(i=0;i<this->m;i++)local_norm+=pow(this->vector[i],2);
740+ MPI_Reduce(&local_norm, &norm, 1, MPI_DOUBLE, MPI_SUM, 0, IssmComm::GetComm());
741+ return sqrt(norm);
742+ break;
743+ default:
744+ _error_("unknown norm !");
745+ break;
746+ }
747 }
748 /*}}}*/
749 /*FUNCTION Scale{{{*/
750@@ -391,4 +466,4 @@
751 }
752 /*}}}*/
753 };
754-#endif //#ifndef _ISSM_MPI_VEC_H_
755+#endif //#ifndef _ISSM_MPI_VEC_H_
756Index: ../trunk-jpl/src/c/toolkits/issm/IssmSolver.cpp
757===================================================================
758--- ../trunk-jpl/src/c/toolkits/issm/IssmSolver.cpp (revision 0)
759+++ ../trunk-jpl/src/c/toolkits/issm/IssmSolver.cpp (revision 14792)
760@@ -0,0 +1,255 @@
761+/*!\file SolverxSeq
762+ * \brief implementation of sequential solver using the GSL librarie
763+ */
764+
765+#ifdef HAVE_CONFIG_H
766+ #include <config.h>
767+#else
768+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
769+#endif
770+#include <cstring>
771+
772+#include "./IssmSolver.h"
773+#include "../../shared/shared.h"
774+#include "../../classes/classes.h"
775+#include "../../include/include.h"
776+#include "../../io/io.h"
777+
778+#ifdef _HAVE_GSL_
779+#include <gsl/gsl_linalg.h>
780+#endif
781+
782+void IssmSolve(IssmVec<IssmDouble>** pout,IssmMat<IssmDouble>* Kffin, IssmVec<IssmDouble>* pfin, Parameters* parameters){/*{{{*/
783+
784+ /*First off, we assume that the type of IssmVec is IssmSeqVec and the type of IssmMat is IssmDenseMat. So downcast: */
785+ IssmSeqVec<IssmDouble>* pf=(IssmSeqVec<IssmDouble>*)pfin->vector;
786+ IssmDenseMat<IssmDouble>* Kff=(IssmDenseMat<IssmDouble>*)Kffin->matrix;
787+
788+#ifdef _HAVE_GSL_
789+ /*Intermediary: */
790+ int M,N,N2;
791+ IssmSeqVec<IssmDouble> *uf = NULL;
792+
793+ Kff->GetSize(&M,&N);
794+ pf->GetSize(&N2);
795+
796+ if(N!=N2)_error_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !");
797+ if(M!=N)_error_("Stiffness matrix should be square!");
798+#ifdef _HAVE_ADOLC_
799+ ensureContiguousLocations(N);
800+#endif
801+ IssmDouble *x = xNew<IssmDouble>(N);
802+#ifdef _HAVE_ADOLC_
803+ SolverxSeq(x,Kff->matrix,pf->vector,N,parameters);
804+#else
805+ SolverxSeq(x,Kff->matrix,pf->vector,N);
806+#endif
807+
808+ uf=new IssmSeqVec<IssmDouble>(x,N);
809+ xDelete(x);
810+
811+ /*Assign output pointers:*/
812+ IssmVec<IssmDouble>* out=new IssmVec<IssmDouble>();
813+ out->vector=(IssmAbsVec<IssmDouble>*)uf;
814+ *pout=out;
815+
816+#else
817+ _error_("GSL support not compiled in!");
818+#endif
819+
820+}/*}}}*/
821+void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
822+
823+ /*Allocate output*/
824+ double* X = xNew<double>(n);
825+
826+ /*Solve*/
827+ SolverxSeq(X,A,B,n);
828+
829+ /*Assign output pointer*/
830+ *pX=X;
831+}
832+/*}}}*/
833+void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
834+#ifdef _HAVE_GSL_
835+ /*GSL Matrices and vectors: */
836+ int s;
837+ gsl_matrix_view a;
838+ gsl_vector_view b,x;
839+ gsl_permutation *p = NULL;
840+// for (int i=0; i<n*n; ++i) std::cout << "SolverxSeq A["<< i << "]=" << A[i] << std::endl;
841+// for (int i=0; i<n; ++i) std::cout << "SolverxSeq b["<< i << "]=" << B[i] << std::endl;
842+ /*A will be modified by LU decomposition. Use copy*/
843+ double* Acopy = xNew<double>(n*n);
844+ xMemCpy(Acopy,A,n*n);
845+
846+ /*Initialize gsl matrices and vectors: */
847+ a = gsl_matrix_view_array (Acopy,n,n);
848+ b = gsl_vector_view_array (B,n);
849+ x = gsl_vector_view_array (X,n);
850+
851+ /*Run LU and solve: */
852+ p = gsl_permutation_alloc (n);
853+ gsl_linalg_LU_decomp (&a.matrix, p, &s);
854+ gsl_linalg_LU_solve (&a.matrix, p, &b.vector, &x.vector);
855+
856+ /*Clean up and assign output pointer*/
857+ xDelete(Acopy);
858+ gsl_permutation_free(p);
859+#endif
860+}
861+/*}}}*/
862+
863+#ifdef _HAVE_ADOLC_
864+int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y){ /*{{{*/
865+ SolverxSeq(y,x, x+m*m, m); // x is where the matrix starts, x+m*m is where the right-hand side starts
866+ return 0;
867+} /*}}}*/
868+int EDF_fos_forward_for_solverx(int n, IssmPDouble *inVal, IssmPDouble *inDeriv, int m, IssmPDouble *outVal, IssmPDouble *outDeriv) { /*{{{*/
869+#ifdef _HAVE_GSL_
870+ // for (int i=0; i<m*m; ++i) std::cout << "EDF_fos_forward_for_solverx A["<< i << "]=" << inVal[i] << std::endl;
871+ // for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx b["<< i << "]=" << inVal[i+m*m] << std::endl;
872+ // the matrix will be modified by LU decomposition. Use gsl_A copy
873+ double* Acopy = xNew<double>(m*m);
874+ xMemCpy(Acopy,inVal,m*m);
875+ /*Initialize gsl matrices and vectors: */
876+ gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
877+ gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
878+ gsl_permutation *perm_p = gsl_permutation_alloc (m);
879+ int signPerm;
880+ // factorize
881+ gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
882+ gsl_vector *gsl_x_p = gsl_vector_alloc (m);
883+ // solve for the value
884+ gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
885+ /*Copy result*/
886+ xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
887+ gsl_vector_free(gsl_x_p);
888+ // for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx x["<< i << "]=" << outVal[i] << std::endl;
889+ // solve for the derivatives acc. to A * dx = r with r=db - dA * x
890+ // compute the RHS
891+ double* r=xNew<double>(m);
892+ for (int i=0; i<m; i++) {
893+ r[i]=inDeriv[m*m+i]; // this is db[i]
894+ for (int j=0;j<m; j++) {
895+ r[i]-=inDeriv[i*m+j]*outVal[j]; // this is dA[i][j]*x[j]
896+ }
897+ }
898+ gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
899+ gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
900+ gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
901+ xMemCpy(outDeriv,gsl_vector_ptr(gsl_dx_p,0),m);
902+ gsl_vector_free(gsl_dx_p);
903+ xDelete(r);
904+ gsl_permutation_free(perm_p);
905+ xDelete(Acopy);
906+#endif
907+ return 0;
908+} /*}}}*/
909+int EDF_fov_forward_for_solverx(int n, IssmPDouble *inVal, int directionCount, IssmPDouble **inDeriv, int m, IssmPDouble *outVal, IssmPDouble **outDeriv) { /*{{{*/
910+#ifdef _HAVE_GSL_
911+ // the matrix will be modified by LU decomposition. Use gsl_A copy
912+ double* Acopy = xNew<double>(m*m);
913+ xMemCpy(Acopy,inVal,m*m);
914+ /*Initialize gsl matrices and vectors: */
915+ gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
916+ gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
917+ gsl_permutation *perm_p = gsl_permutation_alloc (m);
918+ int signPerm;
919+ // factorize
920+ gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
921+ gsl_vector *gsl_x_p = gsl_vector_alloc (m);
922+ // solve for the value
923+ gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
924+ /*Copy result*/
925+ xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
926+ gsl_vector_free(gsl_x_p);
927+ // solve for the derivatives acc. to A * dx = r with r=db - dA * x
928+ double* r=xNew<double>(m);
929+ gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
930+ for (int dir=0;dir<directionCount;++dir) {
931+ // compute the RHS
932+ for (int i=0; i<m; i++) {
933+ r[i]=inDeriv[m*m+i][dir]; // this is db[i]
934+ for (int j=0;j<m; j++) {
935+ r[i]-=inDeriv[i*m+j][dir]*outVal[j]; // this is dA[i][j]*x[j]
936+ }
937+ }
938+ gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
939+ gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
940+ // reuse r
941+ xMemCpy(r,gsl_vector_ptr(gsl_dx_p,0),m);
942+ for (int i=0; i<m; i++) {
943+ outDeriv[i][dir]=r[i];
944+ }
945+ }
946+ gsl_vector_free(gsl_dx_p);
947+ xDelete(r);
948+ gsl_permutation_free(perm_p);
949+ xDelete(Acopy);
950+#endif
951+ return 0;
952+}
953+/*}}}*/
954+int EDF_fos_reverse_for_solverx(int m, double *dp_U, int n, double *dp_Z, double* dp_x, double* dp_y) { /*{{{*/
955+ // copy to transpose the matrix
956+ double* transposed=xNew<double>(m*m);
957+ for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j];
958+ gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m);
959+ // the adjoint of the solution is our right-hand side
960+ gsl_vector_view x_bar=gsl_vector_view_array(dp_U,m);
961+ // the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute:
962+ gsl_vector_view b_bar=gsl_vector_view_array(dp_Z+m*m,m);
963+ gsl_permutation *perm_p = gsl_permutation_alloc (m);
964+ int permSign;
965+ gsl_linalg_LU_decomp (&aTransposed.matrix, perm_p, &permSign);
966+ gsl_linalg_LU_solve (&aTransposed.matrix, perm_p, &x_bar.vector, &b_bar.vector);
967+ // now do the adjoint of the matrix
968+ for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) dp_Z[i*m+j]-=dp_Z[m*m+i]*dp_y[j];
969+ gsl_permutation_free(perm_p);
970+ xDelete(transposed);
971+ return 0;
972+}
973+/*}}}*/
974+int EDF_fov_reverse_for_solverx(int m, int p, double **dpp_U, int n, double **dpp_Z, double* dp_x, double* dp_y) { /*{{{*/
975+ // copy to transpose the matrix
976+ double* transposed=xNew<double>(m*m);
977+ for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j];
978+ gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m);
979+ gsl_permutation *perm_p = gsl_permutation_alloc (m);
980+ int permSign;
981+ gsl_linalg_LU_decomp (&aTransposed.matrix, perm_p, &permSign);
982+ for (int weightsRowIndex=0;weightsRowIndex<p;++weightsRowIndex) {
983+ // the adjoint of the solution is our right-hand side
984+ gsl_vector_view x_bar=gsl_vector_view_array(dpp_U[weightsRowIndex],m);
985+ // the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute:
986+ gsl_vector_view b_bar=gsl_vector_view_array(dpp_Z[weightsRowIndex]+m*m,m);
987+ gsl_linalg_LU_solve (&aTransposed.matrix, perm_p, &x_bar.vector, &b_bar.vector);
988+ // now do the adjoint of the matrix
989+ for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) dpp_Z[weightsRowIndex][i*m+j]-=dpp_Z[weightsRowIndex][m*m+i]*dp_y[j];
990+ }
991+ gsl_permutation_free(perm_p);
992+ xDelete(transposed);
993+ return 0;
994+}
995+/*}}}*/
996+void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
997+ // pack inputs to conform to the EDF-prescribed interface
998+ // ensure a contiguous block of locations:
999+ ensureContiguousLocations(n*(n+1));
1000+ IssmDouble* adoubleEDFin=xNew<IssmDouble>(n*(n+1)); // packed inputs, i.e. matrix and right hand side
1001+ for(int i=0; i<n*n;i++)adoubleEDFin[i] =A[i]; // pack matrix
1002+ for(int i=0; i<n; i++)adoubleEDFin[i+n*n]=B[i]; // pack the right hand side
1003+ IssmPDouble* pdoubleEDFin=xNew<IssmPDouble>(n*(n+1)); // provide space to transfer inputs during call_ext_fct
1004+ IssmPDouble* pdoubleEDFout=xNew<IssmPDouble>(n); // provide space to transfer outputs during call_ext_fct
1005+ // call the wrapped solver through the registry entry we retrieve from parameters
1006+ call_ext_fct(dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p,
1007+ n*(n+1), pdoubleEDFin, adoubleEDFin,
1008+ n, pdoubleEDFout,X);
1009+ // for(int i=0; i<n; i++) {ADOLC_DUMP_MACRO(X[i]);}
1010+ xDelete(adoubleEDFin);
1011+ xDelete(pdoubleEDFin);
1012+ xDelete(pdoubleEDFout);
1013+}
1014+/*}}}*/
1015+#endif
1016Index: ../trunk-jpl/src/c/toolkits/issm/issmtoolkit.h
1017===================================================================
1018--- ../trunk-jpl/src/c/toolkits/issm/issmtoolkit.h (revision 14791)
1019+++ ../trunk-jpl/src/c/toolkits/issm/issmtoolkit.h (revision 14792)
1020@@ -17,6 +17,7 @@
1021 #include "./IssmMat.h"
1022 #include "./IssmSeqVec.h"
1023 #include "./IssmVec.h"
1024+#include "./IssmSolver.h"
1025
1026 #ifdef _HAVE_MPI_
1027 #include "./IssmMpiDenseMat.h"
1028Index: ../trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h
1029===================================================================
1030--- ../trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h (revision 14791)
1031+++ ../trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h (revision 14792)
1032@@ -271,27 +271,80 @@
1033 /*}}}*/
1034 /*FUNCTION Norm{{{*/
1035 doubletype Norm(NormMode mode){
1036- _error_("not supported yet!");
1037+
1038+
1039+ doubletype norm,local_norm;
1040+ doubletype absolute;
1041+ int i,j;
1042+
1043+ switch(mode){
1044+ case NORM_INF:
1045+ local_norm=0;
1046+ for(i=0;i<this->M;i++){
1047+ absolute=0;
1048+ for(j=0;j<this->N;j++){
1049+ absolute+=fabs(this->matrix[N*i+j]);
1050+ }
1051+ local_norm=max(local_norm,absolute);
1052+ }
1053+ MPI_Reduce(&local_norm, &norm, 1, MPI_DOUBLE, MPI_MAX, 0, IssmComm::GetComm());
1054+ MPI_Bcast(&norm,1,MPI_DOUBLE,0,IssmComm::GetComm());
1055+ return norm;
1056+ break;
1057+ default:
1058+ _error_("unknown norm !");
1059+ break;
1060+ }
1061 }
1062 /*}}}*/
1063 /*FUNCTION GetSize{{{*/
1064 void GetSize(int* pM,int* pN){
1065- _error_("not supported yet!");
1066+ *pM=M;
1067+ *pN=N;
1068 }
1069 /*}}}*/
1070 /*FUNCTION GetLocalSize{{{*/
1071 void GetLocalSize(int* pM,int* pN){
1072- _error_("not supported yet!");
1073+ *pM=m;
1074+ *pN=N;
1075 }
1076 /*}}}*/
1077 /*FUNCTION MatMult{{{*/
1078 void MatMult(IssmAbsVec<doubletype>* Xin,IssmAbsVec<doubletype>* AXin){
1079- _error_("not supported yet!");
1080+
1081+
1082+ int i,j;
1083+ doubletype *X_serial = NULL;
1084+
1085+
1086+ /*A check on the types: */
1087+ if(IssmVecTypeFromToolkitOptions()!=MpiEnum)_error_("MatMult operation only possible with 'mpi' vectors");
1088+
1089+ /*Now that we are sure, cast vectors: */
1090+ IssmMpiVec<doubletype>* X=(IssmMpiVec<doubletype>*)Xin;
1091+ IssmMpiVec<doubletype>* AX=(IssmMpiVec<doubletype>*)AXin;
1092+
1093+ /*Serialize input Xin: */
1094+ X_serial=X->ToMPISerial();
1095+
1096+ /*Every cpu has a serial version of the input vector. Use it to do the Matrix-Vector multiply
1097+ *locally and plug it into AXin: */
1098+ for(i=0;i<this->m;i++){
1099+ for(j=0;j<this->N;j++){
1100+ AX->vector[i]+=this->matrix[i*N+j]*X_serial[j];
1101+ }
1102+ }
1103+
1104+ /*Free ressources: */
1105+ xDelete<doubletype>(X_serial);
1106 }
1107 /*}}}*/
1108 /*FUNCTION Duplicate{{{*/
1109 IssmMpiDenseMat<doubletype>* Duplicate(void){
1110- _error_("not supported yet!");
1111+
1112+ IssmMpiDenseMat<doubletype>* dup=new IssmMpiDenseMat<doubletype>(this->matrix,this->M,this->N,0);
1113+ return dup;
1114+
1115 }
1116 /*}}}*/
1117 /*FUNCTION ToSerial{{{*/
1118@@ -313,7 +366,7 @@
1119 /*}}}*/
1120 /*FUNCTION Convert{{{*/
1121 void Convert(MatrixType type){
1122- /*do nothing: */
1123+ _error_("not supported yet!");
1124 }
1125 /*}}}*/
1126 };
1127Index: ../trunk-jpl/src/c/toolkits/issm/IssmSolver.h
1128===================================================================
1129--- ../trunk-jpl/src/c/toolkits/issm/IssmSolver.h (revision 0)
1130+++ ../trunk-jpl/src/c/toolkits/issm/IssmSolver.h (revision 14792)
1131@@ -0,0 +1,38 @@
1132+/*!\file: IssmSolver.h
1133+ * \brief main hook up from Solver toolkit object (in src/c/classes/toolkits) to the ISSM toolkit
1134+ */
1135+
1136+#ifndef _ISSM_SOLVER_H_
1137+#define _ISSM_SOLVER_H_
1138+
1139+/*Headers:*/
1140+/*{{{*/
1141+#ifdef HAVE_CONFIG_H
1142+ #include <config.h>
1143+#else
1144+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
1145+#endif
1146+
1147+#include "../../include/types.h"
1148+
1149+/*}}}*/
1150+
1151+template <class doubletype> class IssmVec;
1152+template <class doubletype> class IssmMat;
1153+class Parameters;
1154+
1155+void IssmSolve(IssmVec<IssmDouble>** puf,IssmMat<IssmDouble>* Kff, IssmVec<IssmDouble>* pf,Parameters* parameters);
1156+void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n);
1157+void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n);
1158+
1159+#if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_)
1160+void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters);
1161+// call back functions:
1162+ADOLC_ext_fct EDF_for_solverx;
1163+ADOLC_ext_fct_fos_forward EDF_fos_forward_for_solverx;
1164+ADOLC_ext_fct_fos_reverse EDF_fos_reverse_for_solverx;
1165+ADOLC_ext_fct_fov_forward EDF_fov_forward_for_solverx;
1166+ADOLC_ext_fct_fov_reverse EDF_fov_reverse_for_solverx;
1167+#endif
1168+
1169+#endif
1170Index: ../trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.h
1171===================================================================
1172--- ../trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.h (revision 0)
1173+++ ../trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.h (revision 14792)
1174@@ -0,0 +1,24 @@
1175+/*!\file: PetscMat.h
1176+ */
1177+
1178+#ifndef _PETSC_SOLVER_H_
1179+#define _PETSC_SOLVER_H_
1180+
1181+/*Headers:*/
1182+/*{{{*/
1183+#ifdef HAVE_CONFIG_H
1184+ #include <config.h>
1185+#else
1186+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
1187+#endif
1188+
1189+#include "../petscincludes.h"
1190+class Parameters;
1191+
1192+/*}}}*/
1193+
1194+void PetscSolve(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters);
1195+void SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters);
1196+void DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum);
1197+
1198+#endif //#ifndef _PETSCSOLVER_H_
1199Index: ../trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.cpp
1200===================================================================
1201--- ../trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.cpp (revision 0)
1202+++ ../trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.cpp (revision 14792)
1203@@ -0,0 +1,232 @@
1204+/*!\file SolverxPetsc
1205+ * \brief Petsc implementation of solver
1206+ */
1207+
1208+
1209+#ifdef HAVE_CONFIG_H
1210+ #include <config.h>
1211+#else
1212+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
1213+#endif
1214+
1215+#include "./PetscSolver.h"
1216+#include "../../../shared/Numerics/Verbosity.h"
1217+#include "../../../shared/Alloc/alloc.h"
1218+#include "../../../shared/Exceptions/exceptions.h"
1219+#include "../../../classes/IssmComm.h"
1220+#include "../../../EnumDefinitions/EnumDefinitions.h"
1221+
1222+void PetscSolve(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters){ /*{{{*/
1223+
1224+ PetscVec* uf=new PetscVec();
1225+
1226+ Vec uf0_vector = NULL;
1227+ Vec df_vector = NULL;
1228+
1229+ if(uf0) uf0_vector = uf0->vector;
1230+ if(df) df_vector = df->vector;
1231+
1232+ SolverxPetsc(&uf->vector, Kff->matrix, pf->vector, uf0_vector, df_vector, parameters);
1233+
1234+ *puf=uf;
1235+
1236+}
1237+/*}}}*/
1238+void SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){ /*{{{*/
1239+
1240+ /*Output: */
1241+ Vec uf = NULL;
1242+
1243+ /*Intermediary: */
1244+ int local_m,local_n,global_m,global_n;
1245+
1246+ /*Solver */
1247+ KSP ksp = NULL;
1248+ PC pc = NULL;
1249+ int iteration_number;
1250+ int solver_type;
1251+ bool fromlocalsize = true;
1252+ #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
1253+ PetscTruth flag,flg;
1254+ #else
1255+ PetscBool flag,flg;
1256+ #endif
1257+
1258+ /*Stokes: */
1259+ IS isv=NULL;
1260+ IS isp=NULL;
1261+
1262+ #if _PETSC_MAJOR_ >= 3
1263+ char ksp_type[50];
1264+ #endif
1265+
1266+ /*Display message*/
1267+ #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
1268+ if(VerboseSolver())PetscOptionsPrint(stdout);
1269+ #else
1270+ if(VerboseSolver())PetscOptionsView(PETSC_VIEWER_STDOUT_WORLD);
1271+ #endif
1272+
1273+ /*First, check that f-set is not NULL, i.e. model is fully constrained:*/
1274+ _assert_(Kff);
1275+ MatGetSize(Kff,&global_m,&global_n); _assert_(global_m==global_n);
1276+ if(!global_n){
1277+ *puf=NewVec(0,IssmComm::GetComm()); return;
1278+ }
1279+
1280+ /*Initial guess */
1281+ /*Now, check that we are not giving an initial guess to the solver, if we are running a direct solver: */
1282+ #if _PETSC_MAJOR_ >= 3
1283+ PetscOptionsGetString(PETSC_NULL,"-ksp_type",ksp_type,49,&flg);
1284+ if (strcmp(ksp_type,"preonly")==0)uf0=NULL;
1285+ #endif
1286+
1287+ /*If initial guess for the solution exists, use it to create uf, otherwise,
1288+ * duplicate the right hand side so that the solution vector has the same structure*/
1289+ if(uf0){
1290+ VecDuplicate(uf0,&uf); VecCopy(uf0,uf);
1291+ }
1292+ else{
1293+ MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,IssmComm::GetComm(),fromlocalsize);
1294+ }
1295+
1296+ /*Process petsc options to see if we are using special types of external solvers*/
1297+ PetscOptionsDetermineSolverType(&solver_type);
1298+
1299+ /*Check the solver is available*/
1300+ if(solver_type==MUMPSPACKAGE_LU || solver_type==MUMPSPACKAGE_CHOL){
1301+ #if _PETSC_MAJOR_ >=3
1302+ #ifndef _HAVE_MUMPS_
1303+ _error_("requested MUMPS solver, which was not compiled into ISSM!\n");
1304+ #endif
1305+ #endif
1306+ }
1307+
1308+ /*Prepare solver*/
1309+ KSPCreate(IssmComm::GetComm(),&ksp);
1310+ KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN);
1311+ KSPSetFromOptions(ksp);
1312+
1313+ #if _PETSC_MAJOR_==3
1314+ /*Specific solver?: */
1315+ KSPGetPC(ksp,&pc);
1316+ if (solver_type==MUMPSPACKAGE_LU){
1317+ #if _PETSC_MINOR_==1
1318+ PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS);
1319+ #else
1320+ PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
1321+ #endif
1322+ }
1323+
1324+ /*Stokes: */
1325+ if (solver_type==StokesSolverEnum){
1326+ /*Make indices out of doftypes: */
1327+ if(!df)_error_("need doftypes for Stokes solver!\n");
1328+ DofTypesToIndexSet(&isv,&isp,df,StokesSolverEnum);
1329+
1330+ /*Set field splits: */
1331+ KSPGetPC(ksp,&pc);
1332+ #if _PETSC_MINOR_==1
1333+ PCFieldSplitSetIS(pc,isv);
1334+ PCFieldSplitSetIS(pc,isp);
1335+ #else
1336+ PCFieldSplitSetIS(pc,PETSC_NULL,isv);
1337+ PCFieldSplitSetIS(pc,PETSC_NULL,isp);
1338+ #endif
1339+
1340+ }
1341+ #endif
1342+
1343+ /*If there is an initial guess for the solution, use it
1344+ * except if we are using the MUMPS direct solver
1345+ * where any initial guess will crash Petsc*/
1346+ if (uf0){
1347+ if((solver_type!=MUMPSPACKAGE_LU) && (solver_type!=MUMPSPACKAGE_CHOL) && (solver_type!=SPOOLESPACKAGE_LU) && (solver_type!=SPOOLESPACKAGE_CHOL) && (solver_type!=SUPERLUDISTPACKAGE)){
1348+ KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
1349+ }
1350+ }
1351+
1352+ /*Solve: */
1353+ if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
1354+ KSPSolve(ksp,pf,uf);
1355+
1356+ /*Check convergence*/
1357+ KSPGetIterationNumber(ksp,&iteration_number);
1358+ if (iteration_number<0) _error_("Solver diverged at iteration number: " << -iteration_number);
1359+
1360+ /*Free resources:*/
1361+ KSPFree(&ksp);
1362+
1363+ /*Assign output pointers:*/
1364+ *puf=uf;
1365+}
1366+/*}}}*/
1367+void DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum){ /*{{{*/
1368+
1369+ /*output: */
1370+ IS isv=NULL;
1371+ IS isp=NULL;
1372+
1373+ int start,end;
1374+ IssmDouble* df_local=NULL;
1375+ int df_local_size;
1376+ int i;
1377+
1378+ int* pressure_indices=NULL;
1379+ int* velocity_indices=NULL;
1380+ int pressure_num=0;
1381+ int velocity_num=0;
1382+ int pressure_count=0;
1383+ int velocity_count=0;
1384+
1385+ if(typeenum==StokesSolverEnum){
1386+
1387+ /*Ok, recover doftypes vector values and indices: */
1388+ VecGetOwnershipRange(df,&start,&end);
1389+ VecGetLocalSize(df,&df_local_size);
1390+ VecGetArray(df,&df_local);
1391+
1392+ pressure_num=0;
1393+ velocity_num=0;
1394+ for(i=0;i<df_local_size;i++){
1395+ if (df_local[i]==PressureEnum)pressure_num++;
1396+ else velocity_num++;
1397+ }
1398+
1399+ /*Allocate indices: */
1400+ if(pressure_num)pressure_indices=xNew<int>(pressure_num);
1401+ if(velocity_num)velocity_indices=xNew<int>(velocity_num);
1402+
1403+ pressure_count=0;
1404+ velocity_count=0;
1405+ for(i=0;i<df_local_size;i++){
1406+ if (df_local[i]==PressureEnum){
1407+ pressure_indices[pressure_count]=start+i;
1408+ pressure_count++;
1409+ }
1410+ if (df_local[i]==VelocityEnum){
1411+ velocity_indices[velocity_count]=start+i;
1412+ velocity_count++;
1413+ }
1414+ }
1415+ VecRestoreArray(df,&df_local);
1416+
1417+ /*Create indices sets: */
1418+ #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
1419+ ISCreateGeneral(IssmComm::GetComm(),pressure_num,pressure_indices,&isp);
1420+ ISCreateGeneral(IssmComm::GetComm(),velocity_num,velocity_indices,&isv);
1421+ #else
1422+ ISCreateGeneral(IssmComm::GetComm(),pressure_num,pressure_indices,PETSC_COPY_VALUES,&isp);
1423+ ISCreateGeneral(IssmComm::GetComm(),velocity_num,velocity_indices,PETSC_COPY_VALUES,&isv);
1424+ #endif
1425+ }
1426+
1427+ /*Free ressources:*/
1428+ xDelete<int>(pressure_indices);
1429+ xDelete<int>(velocity_indices);
1430+
1431+ /*Assign output pointers:*/
1432+ *pisv=isv;
1433+ *pisp=isp;
1434+}
1435+/*}}}*/
1436Index: ../trunk-jpl/src/c/toolkits/petsc/objects/petscobjects.h
1437===================================================================
1438--- ../trunk-jpl/src/c/toolkits/petsc/objects/petscobjects.h (revision 14791)
1439+++ ../trunk-jpl/src/c/toolkits/petsc/objects/petscobjects.h (revision 14792)
1440@@ -7,5 +7,6 @@
1441
1442 #include "./PetscMat.h"
1443 #include "./PetscVec.h"
1444+#include "./PetscSolver.h"
1445
1446 #endif
1447Index: ../trunk-jpl/src/c/Makefile.am
1448===================================================================
1449--- ../trunk-jpl/src/c/Makefile.am (revision 14791)
1450+++ ../trunk-jpl/src/c/Makefile.am (revision 14792)
1451@@ -118,8 +118,10 @@
1452 ./classes/matrix/ElementMatrix.cpp\
1453 ./classes/matrix/ElementVector.h\
1454 ./classes/matrix/ElementVector.cpp\
1455- ./classes/matrix/Matrix.h\
1456- ./classes/matrix/Vector.h\
1457+ ./classes/toolkits/toolkitobjects.h\
1458+ ./classes/toolkits/Matrix.h\
1459+ ./classes/toolkits/Vector.h\
1460+ ./classes/toolkits/Solver.h\
1461 ./classes/objects/Params/Param.h\
1462 ./classes/objects/Params/GenericParam.h\
1463 ./classes/objects/Params/BoolParam.cpp\
1464@@ -228,6 +230,8 @@
1465 ./toolkits/issm/IssmMat.h\
1466 ./toolkits/issm/IssmSeqVec.h\
1467 ./toolkits/issm/IssmVec.h\
1468+ ./toolkits/issm/IssmSolver.h\
1469+ ./toolkits/issm/IssmSolver.cpp\
1470 ./toolkits/triangle/triangleincludes.h\
1471 ./toolkitsenums.h\
1472 ./toolkits.h\
1473@@ -319,7 +323,6 @@
1474 ./modules/ResetCoordinateSystemx/ResetCoordinateSystemx.cpp\
1475 ./modules/Solverx/Solverx.cpp\
1476 ./modules/Solverx/Solverx.h\
1477- ./modules/Solverx/SolverxSeq.cpp\
1478 ./modules/VecMergex/VecMergex.cpp\
1479 ./modules/VecMergex/VecMergex.h\
1480 ./modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp\
1481@@ -759,9 +762,9 @@
1482 ./toolkits/petsc/objects/PetscMat.cpp\
1483 ./toolkits/petsc/objects/PetscVec.h\
1484 ./toolkits/petsc/objects/PetscVec.cpp\
1485- ./toolkits/petsc/petscincludes.h\
1486- ./modules/Solverx/SolverxPetsc.cpp\
1487- ./modules/Solverx/DofTypesToIndexSet.cpp
1488+ ./toolkits/petsc/objects/PetscSolver.cpp\
1489+ ./toolkits/petsc/objects/PetscSolver.h\
1490+ ./toolkits/petsc/petscincludes.h
1491
1492 #}}}
1493 #Mpi sources {{{
1494Index: ../trunk-jpl/src/c/classes/classes.h
1495===================================================================
1496--- ../trunk-jpl/src/c/classes/classes.h (revision 14791)
1497+++ ../trunk-jpl/src/c/classes/classes.h (revision 14792)
1498@@ -11,6 +11,9 @@
1499 /*matrix: */
1500 #include "./matrix/matrixobjects.h"
1501
1502+/*toolkit objects: */
1503+#include "./toolkits/toolkitobjects.h"
1504+
1505 /*bamg: */
1506 #include "./bamg/bamgobjects.h"
1507
1508Index: ../trunk-jpl/src/c/classes/matrix/Vector.h
1509===================================================================
1510--- ../trunk-jpl/src/c/classes/matrix/Vector.h (revision 14791)
1511+++ ../trunk-jpl/src/c/classes/matrix/Vector.h (revision 14792)
1512@@ -1,357 +0,0 @@
1513-/*!\file: Vector.h
1514- * \brief wrapper to vector objects. The goal is to control which API (PETSc,Scalpack, Plapack?)
1515- * implements our underlying vector format.
1516- */
1517-
1518-#ifndef _VECTOR_H_
1519-#define _VECTOR_H_
1520-
1521-/*Headers:*/
1522-/*{{{*/
1523-#ifdef HAVE_CONFIG_H
1524- #include <config.h>
1525-#else
1526-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
1527-#endif
1528-#include <cstring>
1529-#include "../../toolkits/toolkits.h"
1530-#include "../../EnumDefinitions/EnumDefinitions.h"
1531-/*}}}*/
1532-
1533-enum vectortype { PetscVecType, IssmVecType };
1534-
1535-template <class doubletype>
1536-class Vector{
1537-
1538- public:
1539-
1540- int type;
1541- #ifdef _HAVE_PETSC_
1542- PetscVec* pvector;
1543- #endif
1544- IssmVec<doubletype>* ivector;
1545-
1546- /*Vector constructors, destructors */
1547- Vector(){ /*{{{*/
1548-
1549- InitCheckAndSetType();
1550- }
1551- /*}}}*/
1552- Vector(int M,bool fromlocalsize=false){ /*{{{*/
1553-
1554- InitCheckAndSetType();
1555-
1556- if(type==PetscVecType){
1557- #ifdef _HAVE_PETSC_
1558- this->pvector=new PetscVec(M,fromlocalsize);
1559- #endif
1560- }
1561- else this->ivector=new IssmVec<doubletype>(M,fromlocalsize);
1562-
1563- }
1564- /*}}}*/
1565- Vector(int m,int M){ /*{{{*/
1566-
1567- InitCheckAndSetType();
1568-
1569- if(type==PetscVecType){
1570- #ifdef _HAVE_PETSC_
1571- this->pvector=new PetscVec(m,M);
1572- #endif
1573- }
1574- else this->ivector=new IssmVec<doubletype>(m,M);
1575- }
1576- /*}}}*/
1577- Vector(doubletype* serial_vec,int M){ /*{{{*/
1578-
1579- InitCheckAndSetType();
1580-
1581- if(type==PetscVecType){
1582- #ifdef _HAVE_PETSC_
1583- this->pvector=new PetscVec(serial_vec,M);
1584- #endif
1585- }
1586- else this->ivector=new IssmVec<doubletype>(serial_vec,M);
1587- }
1588- /*}}}*/
1589- ~Vector(){ /*{{{*/
1590-
1591- if(type==PetscVecType){
1592- #ifdef _HAVE_PETSC_
1593- delete this->pvector;
1594- #endif
1595- }
1596- else delete this->ivector;
1597- }
1598- /*}}}*/
1599- #ifdef _HAVE_PETSC_
1600- Vector(Vec petsc_vector){ /*{{{*/
1601-
1602- this->type=PetscVecType;
1603- this->ivector=NULL;
1604- this->pvector=new PetscVec(petsc_vector);
1605-
1606- }
1607- /*}}}*/
1608- #endif
1609- void InitCheckAndSetType(void){ /*{{{*/
1610-
1611- char* toolkittype=NULL;
1612-
1613- #ifdef _HAVE_PETSC_
1614- pvector=NULL;
1615- #endif
1616- ivector=NULL;
1617-
1618- /*retrieve toolkittype: */
1619- toolkittype=ToolkitOptions::GetToolkitType();
1620-
1621- /*set vector type: */
1622- if (strcmp(toolkittype,"petsc")==0){
1623- #ifdef _HAVE_PETSC_
1624- type=PetscVecType;
1625- #else
1626- _error_("cannot create petsc vector without PETSC compiled!");
1627- #endif
1628- }
1629- else if(strcmp(toolkittype,"issm")==0){
1630- /*let this choice stand:*/
1631- type=IssmVecType;
1632- }
1633- else _error_("unknow toolkit type ");
1634-
1635- /*Free ressources: */
1636- xDelete<char>(toolkittype);
1637- }
1638- /*}}}*/
1639-
1640- /*Vector specific routines*/
1641- /*FUNCTION Echo{{{*/
1642- void Echo(void){_assert_(this);
1643-
1644- if(type==PetscVecType){
1645- #ifdef _HAVE_PETSC_
1646- this->pvector->Echo();
1647- #endif
1648- }
1649- else this->ivector->Echo();
1650-
1651- }
1652- /*}}}*/
1653- /*FUNCTION Assemble{{{*/
1654- void Assemble(void){_assert_(this);
1655-
1656- if(type==PetscVecType){
1657- #ifdef _HAVE_PETSC_
1658- this->pvector->Assemble();
1659- #endif
1660- }
1661- else this->ivector->Assemble();
1662-
1663- }
1664- /*}}}*/
1665- /*FUNCTION SetValues{{{*/
1666- void SetValues(int ssize, int* list, doubletype* values, InsMode mode){ _assert_(this);
1667- if(type==PetscVecType){
1668- #ifdef _HAVE_PETSC_
1669- this->pvector->SetValues(ssize,list,values,mode);
1670- #endif
1671- }
1672- else this->ivector->SetValues(ssize,list,values,mode);
1673-
1674- }
1675- /*}}}*/
1676- /*FUNCTION SetValue{{{*/
1677- void SetValue(int dof, doubletype value, InsMode mode){_assert_(this);
1678-
1679- if(type==PetscVecType){
1680- #ifdef _HAVE_PETSC_
1681- this->pvector->SetValue(dof,value,mode);
1682- #endif
1683- }
1684- else this->ivector->SetValue(dof,value,mode);
1685-
1686- }
1687- /*}}}*/
1688- /*FUNCTION GetValue{{{*/
1689- void GetValue(doubletype* pvalue,int dof){_assert_(this);
1690-
1691- if(type==PetscVecType){
1692- #ifdef _HAVE_PETSC_
1693- this->pvector->GetValue(pvalue,dof);
1694- #endif
1695- }
1696- else this->ivector->GetValue(pvalue,dof);
1697-
1698- }
1699- /*}}}*/
1700- /*FUNCTION GetSize{{{*/
1701- void GetSize(int* pM){_assert_(this);
1702-
1703- if(type==PetscVecType){
1704- #ifdef _HAVE_PETSC_
1705- this->pvector->GetSize(pM);
1706- #endif
1707- }
1708- else this->ivector->GetSize(pM);
1709-
1710- }
1711- /*}}}*/
1712- /*FUNCTION IsEmpty{{{*/
1713- bool IsEmpty(void){
1714- int M;
1715-
1716- _assert_(this);
1717- this->GetSize(&M);
1718-
1719- if(M==0)
1720- return true;
1721- else
1722- return false;
1723- }
1724- /*}}}*/
1725- /*FUNCTION GetLocalSize{{{*/
1726- void GetLocalSize(int* pM){_assert_(this);
1727-
1728- if(type==PetscVecType){
1729- #ifdef _HAVE_PETSC_
1730- this->pvector->GetLocalSize(pM);
1731- #endif
1732- }
1733- else this->ivector->GetLocalSize(pM);
1734-
1735- }
1736- /*}}}*/
1737- /*FUNCTION Duplicate{{{*/
1738- Vector<doubletype>* Duplicate(void){_assert_(this);
1739-
1740- Vector<doubletype>* output=NULL;
1741-
1742- output=new Vector<doubletype>();
1743-
1744- if(type==PetscVecType){
1745- #ifdef _HAVE_PETSC_
1746- output->pvector=this->pvector->Duplicate();
1747- #endif
1748- }
1749- else output->ivector=this->ivector->Duplicate();
1750-
1751- return output;
1752-
1753- }
1754- /*}}}*/
1755- /*FUNCTION Set{{{*/
1756- void Set(doubletype value){_assert_(this);
1757-
1758- if(type==PetscVecType){
1759- #ifdef _HAVE_PETSC_
1760- this->pvector->Set(value);
1761- #endif
1762- }
1763- else this->ivector->Set(value);
1764-
1765- }
1766- /*}}}*/
1767- /*FUNCTION AXPY{{{*/
1768- void AXPY(Vector* X, doubletype a){_assert_(this);
1769-
1770- if(type==PetscVecType){
1771- #ifdef _HAVE_PETSC_
1772- this->pvector->AXPY(X->pvector,a);
1773- #endif
1774- }
1775- else this->ivector->AXPY(X->ivector,a);
1776-
1777- }
1778- /*}}}*/
1779- /*FUNCTION AYPX{{{*/
1780- void AYPX(Vector* X, doubletype a){_assert_(this);
1781-
1782- if(type==PetscVecType){
1783- #ifdef _HAVE_PETSC_
1784- this->pvector->AYPX(X->pvector,a);
1785- #endif
1786- }
1787- else this->ivector->AYPX(X->ivector,a);
1788- }
1789- /*}}}*/
1790- /*FUNCTION ToMPISerial{{{*/
1791- doubletype* ToMPISerial(void){
1792-
1793- doubletype* vec_serial=NULL;
1794-
1795- _assert_(this);
1796- if(type==PetscVecType){
1797- #ifdef _HAVE_PETSC_
1798- vec_serial=this->pvector->ToMPISerial();
1799- #endif
1800- }
1801- else vec_serial=this->ivector->ToMPISerial();
1802-
1803- return vec_serial;
1804-
1805- }
1806- /*}}}*/
1807- /*FUNCTION Copy{{{*/
1808- void Copy(Vector* to){_assert_(this);
1809-
1810- if(type==PetscVecType){
1811- #ifdef _HAVE_PETSC_
1812- this->pvector->Copy(to->pvector);
1813- #endif
1814- }
1815- else this->ivector->Copy(to->ivector);
1816- }
1817- /*}}}*/
1818- /*FUNCTION Norm{{{*/
1819- doubletype Norm(NormMode norm_type){_assert_(this);
1820-
1821- doubletype norm=0;
1822-
1823- if(type==PetscVecType){
1824- #ifdef _HAVE_PETSC_
1825- norm=this->pvector->Norm(norm_type);
1826- #endif
1827- }
1828- else norm=this->ivector->Norm(norm_type);
1829- return norm;
1830- }
1831- /*}}}*/
1832- /*FUNCTION Scale{{{*/
1833- void Scale(doubletype scale_factor){_assert_(this);
1834-
1835- if(type==PetscVecType){
1836- #ifdef _HAVE_PETSC_
1837- this->pvector->Scale(scale_factor);
1838- #endif
1839- }
1840- else this->ivector->Scale(scale_factor);
1841- }
1842- /*}}}*/
1843- /*FUNCTION Dot{{{*/
1844- doubletype Dot(Vector* vector){_assert_(this);
1845-
1846- doubletype dot;
1847-
1848- if(type==PetscVecType){
1849- #ifdef _HAVE_PETSC_
1850- dot=this->pvector->Dot(vector->pvector);
1851- #endif
1852- }
1853- else dot=this->ivector->Dot(vector->ivector);
1854- return dot;
1855- }
1856- /*}}}*/
1857- /*FUNCTION PointwiseDivide{{{*/
1858- void PointwiseDivide(Vector* x,Vector* y){_assert_(this);
1859-
1860- if(type==PetscVecType){
1861- #ifdef _HAVE_PETSC_
1862- this->pvector->PointwiseDivide(x->pvector,y->pvector);
1863- #endif
1864- }
1865- else this->ivector->PointwiseDivide(x->ivector,y->ivector);
1866- }
1867- /*}}}*/
1868-};
1869-#endif //#ifndef _VECTOR_H_
1870Index: ../trunk-jpl/src/c/classes/matrix/Matrix.h
1871===================================================================
1872--- ../trunk-jpl/src/c/classes/matrix/Matrix.h (revision 14791)
1873+++ ../trunk-jpl/src/c/classes/matrix/Matrix.h (revision 14792)
1874@@ -1,328 +0,0 @@
1875-/*!\file: Matrix.h
1876- * \brief wrapper to matrix objects. The goal is to control which API (PETSc,Scalpack, Plapack?)
1877- * implements our underlying matrix format.
1878- */
1879-
1880-#ifndef _MATRIX_H_
1881-#define _MATRIX_H_
1882-
1883-/*Headers:*/
1884-/*{{{*/
1885-#ifdef HAVE_CONFIG_H
1886- #include <config.h>
1887-#else
1888-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
1889-#endif
1890-#include <cstring>
1891-#include "../../toolkits/toolkits.h"
1892-#include "../../EnumDefinitions/EnumDefinitions.h"
1893-/*}}}*/
1894-
1895-enum matrixtype { PetscMatType, IssmMatType };
1896-
1897-template <class doubletype> class Vector;
1898-
1899-template <class doubletype>
1900-class Matrix{
1901-
1902- public:
1903-
1904- int type;
1905- #ifdef _HAVE_PETSC_
1906- PetscMat *pmatrix;
1907- #endif
1908- IssmMat<doubletype> *imatrix;
1909-
1910- /*Matrix constructors, destructors*/
1911- /*FUNCTION Matrix(){{{*/
1912- Matrix(){
1913- InitCheckAndSetType();
1914- }
1915- /*}}}*/
1916- /*FUNCTION Matrix(int M,int N){{{*/
1917- Matrix(int M,int N){
1918-
1919- InitCheckAndSetType();
1920-
1921- if(type==PetscMatType){
1922- #ifdef _HAVE_PETSC_
1923- this->pmatrix=new PetscMat(M,N);
1924- #endif
1925- }
1926- else{
1927- this->imatrix=new IssmMat<doubletype>(M,N);
1928- }
1929-
1930- }
1931- /*}}}*/
1932- /*FUNCTION Matrix(int m,int n,int M,int N,int* d_nnz,int* o_nnz){{{*/
1933- Matrix(int m,int n,int M,int N,int* d_nnz,int* o_nnz){
1934-
1935- InitCheckAndSetType();
1936-
1937- if(type==PetscMatType){
1938- #ifdef _HAVE_PETSC_
1939- this->pmatrix=new PetscMat(m,n,M,N,d_nnz,o_nnz);
1940- #endif
1941- }
1942- else{
1943- this->imatrix=new IssmMat<doubletype>(m,n,M,N,d_nnz,o_nnz);
1944- }
1945-
1946- }
1947- /*}}}*/
1948- /*FUNCTION Matrix(int M,int N,IssmDouble sparsity){{{*/
1949- Matrix(int M,int N,double sparsity){
1950-
1951- InitCheckAndSetType();
1952-
1953- if(type==PetscMatType){
1954- #ifdef _HAVE_PETSC_
1955- this->pmatrix=new PetscMat(M,N,sparsity);
1956- #endif
1957- }
1958- else{
1959- this->imatrix=new IssmMat<doubletype>(M,N,sparsity);
1960- }
1961- }
1962- /*}}}*/
1963- /*FUNCTION Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity){{{*/
1964- Matrix(IssmPDouble* serial_mat,int M,int N,IssmPDouble sparsity){
1965-
1966- InitCheckAndSetType();
1967-
1968- if(type==PetscMatType){
1969- #ifdef _HAVE_PETSC_
1970- this->pmatrix=new PetscMat(serial_mat,M,N,sparsity);
1971- #endif
1972- }
1973- else{
1974- this->imatrix=new IssmMat<doubletype>(serial_mat,M,N,sparsity);
1975- }
1976-
1977- }
1978- /*}}}*/
1979- /*FUNCTION Matrix(int M,int N,int connectivity,int numberofdofspernode){{{*/
1980- Matrix(int M,int N,int connectivity,int numberofdofspernode){
1981-
1982- InitCheckAndSetType();
1983-
1984- if(type==PetscMatType){
1985- #ifdef _HAVE_PETSC_
1986- this->pmatrix=new PetscMat(M,N,connectivity,numberofdofspernode);
1987- #endif
1988- }
1989- else{
1990- this->imatrix=new IssmMat<doubletype>(M,N,connectivity,numberofdofspernode);
1991- }
1992-
1993- }
1994- /*}}}*/
1995- /*FUNCTION ~Matrix(){{{*/
1996- ~Matrix(){
1997-
1998- if(type==PetscMatType){
1999- #ifdef _HAVE_PETSC_
2000- delete this->pmatrix;
2001- #endif
2002- }
2003- else delete this->imatrix;
2004-
2005- }
2006- /*}}}*/
2007- /*FUNCTION InitCheckAndSetType(){{{*/
2008- void InitCheckAndSetType(void){
2009-
2010- char* toolkittype=NULL;
2011-
2012- #ifdef _HAVE_PETSC_
2013- pmatrix=NULL;
2014- #endif
2015- imatrix=NULL;
2016-
2017- /*retrieve toolkittype: */
2018- toolkittype=ToolkitOptions::GetToolkitType();
2019-
2020- /*set matrix type: */
2021- if (strcmp(toolkittype,"petsc")==0){
2022- #ifdef _HAVE_PETSC_
2023- type=PetscMatType;
2024- #else
2025- _error_("cannot create petsc matrix without PETSC compiled!");
2026- #endif
2027- }
2028- else if(strcmp(toolkittype,"issm")==0){
2029- /*let this choice stand:*/
2030- type=IssmMatType;
2031- }
2032- else _error_("unknow toolkit type ");
2033-
2034- /*Free ressources: */
2035- xDelete<char>(toolkittype);
2036- }
2037- /*}}}*/
2038-
2039- /*Matrix specific routines:*/
2040- /*FUNCTION Echo{{{*/
2041- void Echo(void){
2042- _assert_(this);
2043-
2044- if(type==PetscMatType){
2045- #ifdef _HAVE_PETSC_
2046- this->pmatrix->Echo();
2047- #endif
2048- }
2049- else{
2050- this->imatrix->Echo();
2051- }
2052-
2053- }
2054- /*}}}*/
2055- /*FUNCTION AllocationInfo{{{*/
2056- void AllocationInfo(void){
2057- _assert_(this);
2058- if(type==PetscMatType){
2059- #ifdef _HAVE_PETSC_
2060- this->pmatrix->AllocationInfo();
2061- #endif
2062- }
2063- else{
2064- this->imatrix->AllocationInfo();
2065- }
2066- }/*}}}*/
2067- /*FUNCTION Assemble{{{*/
2068- void Assemble(void){
2069-
2070- if(type==PetscMatType){
2071- #ifdef _HAVE_PETSC_
2072- this->pmatrix->Assemble();
2073- #endif
2074- }
2075- else{
2076- this->imatrix->Assemble();
2077- }
2078- }
2079- /*}}}*/
2080- /*FUNCTION Norm{{{*/
2081- IssmDouble Norm(NormMode norm_type){
2082-
2083- IssmDouble norm=0;
2084-
2085- if(type==PetscMatType){
2086- #ifdef _HAVE_PETSC_
2087- norm=this->pmatrix->Norm(norm_type);
2088- #endif
2089- }
2090- else{
2091- norm=this->imatrix->Norm(norm_type);
2092- }
2093-
2094- return norm;
2095- }
2096- /*}}}*/
2097- /*FUNCTION GetSize{{{*/
2098- void GetSize(int* pM,int* pN){
2099-
2100- if(type==PetscMatType){
2101- #ifdef _HAVE_PETSC_
2102- this->pmatrix->GetSize(pM,pN);
2103- #endif
2104- }
2105- else{
2106- this->imatrix->GetSize(pM,pN);
2107- }
2108-
2109- }
2110- /*}}}*/
2111- /*FUNCTION GetLocalSize{{{*/
2112- void GetLocalSize(int* pM,int* pN){
2113-
2114- if(type==PetscMatType){
2115- #ifdef _HAVE_PETSC_
2116- this->pmatrix->GetLocalSize(pM,pN);
2117- #endif
2118- }
2119- else{
2120- this->imatrix->GetLocalSize(pM,pN);
2121- }
2122-
2123- }
2124- /*}}}*/
2125- /*FUNCTION MatMult{{{*/
2126- void MatMult(Vector<doubletype>* X,Vector<doubletype>* AX){
2127-
2128- if(type==PetscMatType){
2129- #ifdef _HAVE_PETSC_
2130- this->pmatrix->MatMult(X->pvector,AX->pvector);
2131- #endif
2132- }
2133- else{
2134- this->imatrix->MatMult(X->ivector,AX->ivector);
2135- }
2136-
2137- }
2138- /*}}}*/
2139- /*FUNCTION Duplicate{{{*/
2140- Matrix<doubletype>* Duplicate(void){
2141-
2142- Matrix<doubletype>* output=new Matrix<doubletype>();
2143-
2144- if(type==PetscMatType){
2145- #ifdef _HAVE_PETSC_
2146- output->pmatrix=this->pmatrix->Duplicate();
2147- #endif
2148- }
2149- else{
2150- output->imatrix=this->imatrix->Duplicate();
2151- }
2152-
2153- return output;
2154- }
2155- /*}}}*/
2156- /*FUNCTION ToSerial{{{*/
2157- doubletype* ToSerial(void){
2158-
2159- doubletype* output=NULL;
2160-
2161- if(type==PetscMatType){
2162- #ifdef _HAVE_PETSC_
2163- output=this->pmatrix->ToSerial();
2164- #endif
2165- }
2166- else{
2167- output=this->imatrix->ToSerial();
2168- }
2169-
2170- return output;
2171- }
2172- /*}}}*/
2173- /*FUNCTION SetValues{{{*/
2174- void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){
2175-
2176- if(type==PetscMatType){
2177- #ifdef _HAVE_PETSC_
2178- this->pmatrix->SetValues(m,idxm,n,idxn,values,mode);
2179- #endif
2180- }
2181- else{
2182- this->imatrix->SetValues(m,idxm,n,idxn,values,mode);
2183- }
2184- }
2185- /*}}}*/
2186- /*FUNCTION Convert{{{*/
2187- void Convert(MatrixType newtype){
2188-
2189- if(type==PetscMatType){
2190- #ifdef _HAVE_PETSC_
2191- this->pmatrix->Convert(newtype);
2192- #endif
2193- }
2194- else{
2195- this->imatrix->Convert(newtype);
2196- }
2197-
2198- }
2199- /*}}}*/
2200-};
2201-
2202-#endif //#ifndef _MATRIX_H_
2203Index: ../trunk-jpl/src/c/classes/matrix/matrixobjects.h
2204===================================================================
2205--- ../trunk-jpl/src/c/classes/matrix/matrixobjects.h (revision 14791)
2206+++ ../trunk-jpl/src/c/classes/matrix/matrixobjects.h (revision 14792)
2207@@ -8,7 +8,5 @@
2208 /*Numerics:*/
2209 #include "./ElementMatrix.h"
2210 #include "./ElementVector.h"
2211-#include "./Vector.h"
2212-#include "./Matrix.h"
2213
2214 #endif
2215Index: ../trunk-jpl/src/c/classes/toolkits/Solver.h
2216===================================================================
2217--- ../trunk-jpl/src/c/classes/toolkits/Solver.h (revision 0)
2218+++ ../trunk-jpl/src/c/classes/toolkits/Solver.h (revision 14792)
2219@@ -0,0 +1,99 @@
2220+/*!\file: Solver.h
2221+ */
2222+
2223+#ifndef _SOLVER_CLASS_H_
2224+#define _SOLVER_CLASS_H_
2225+
2226+/*Headers:*/
2227+/*{{{*/
2228+#ifdef HAVE_CONFIG_H
2229+ #include <config.h>
2230+#else
2231+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
2232+#endif
2233+#include "../../toolkits/toolkits.h"
2234+/*}}}*/
2235+
2236+#include "./Matrix.h"
2237+#include "./Vector.h"
2238+class Parameters;
2239+
2240+template <class doubletype>
2241+class Solver{
2242+
2243+ private:
2244+ Matrix<doubletype>* Kff;
2245+ Vector<doubletype>* pf;
2246+ Vector<doubletype>* uf0;
2247+ Vector<doubletype>* df;
2248+ Parameters* parameters;
2249+
2250+ public:
2251+
2252+ /*Constructors, destructors:*/
2253+ /*Solver(){{{*/
2254+ Solver(){
2255+ }
2256+ /*}}}*/
2257+ /*Solver(Matrix<doubletype>* Kff, Vector<doubletype>* pf, Vector<doubletype>* uf0,Vector<doubletype>* df, Parameters* parameters):{{{*/
2258+ Solver(Matrix<doubletype>* Kff_in, Vector<doubletype>* pf_in, Vector<doubletype>* uf0_in,Vector<doubletype>* df_in, Parameters* parameters_in){
2259+
2260+ /*In debugging mode, check that stiffness matrix and load vectors are not NULL (they can be empty)*/
2261+ _assert_(Kff_in);
2262+ _assert_(pf_in);
2263+
2264+ /*initialize fields: */
2265+ this->Kff=Kff_in;
2266+ this->pf=pf_in;
2267+ this->uf0=uf0_in;
2268+ this->df=df_in;
2269+ this->parameters=parameters_in;
2270+ }
2271+ /*}}}*/
2272+ /*~Solver(){{{*/
2273+ ~Solver(){
2274+ }
2275+ /*}}}*/
2276+
2277+ /*Methods: */
2278+ /*Vector<doubletype Solve(void): {{{*/
2279+ Vector<doubletype>* Solve(){
2280+
2281+ /*output: */
2282+ Vector<doubletype>* uf=NULL;
2283+
2284+ /*Initialize vector: */
2285+ uf=new Vector<doubletype>();
2286+
2287+ /*According to matrix type, use specific solvers: */
2288+ switch(Kff->type){
2289+ #ifdef _HAVE_PETSC_
2290+ case PetscMatType:{
2291+ PetscVec* uf0_vector = NULL;
2292+ PetscVec* df_vector = NULL;
2293+ if(uf0) uf0_vector = uf0->pvector;
2294+ if(df) df_vector = df->pvector;
2295+ PetscSolve(&uf->pvector,Kff->pmatrix,pf->pvector,uf0_vector,df_vector,parameters);
2296+ break;
2297+ }
2298+ #endif
2299+ case IssmMatType:{
2300+ IssmSolve(&uf->ivector,Kff->imatrix,pf->ivector,parameters);
2301+ break;
2302+ }
2303+ default:
2304+ _error_("Matrix type: " << Kff->type << " not supported yet!");
2305+ }
2306+
2307+ /*allocate output pointer: */
2308+ return uf;
2309+
2310+ }
2311+ /*}}}*/
2312+
2313+};
2314+
2315+
2316+
2317+
2318+#endif //#ifndef _SOLVER_H_
2319Index: ../trunk-jpl/src/c/classes/toolkits/toolkitobjects.h
2320===================================================================
2321--- ../trunk-jpl/src/c/classes/toolkits/toolkitobjects.h (revision 0)
2322+++ ../trunk-jpl/src/c/classes/toolkits/toolkitobjects.h (revision 14792)
2323@@ -0,0 +1,13 @@
2324+
2325+/*!\file: toolkitobjects.h
2326+ * \brief wrappers to toolkits
2327+ */
2328+
2329+#ifndef _TOOLKIT_OBJECTS_H_
2330+#define _TOOLKIT_OBJECTS_H_
2331+
2332+#include "./Vector.h"
2333+#include "./Matrix.h"
2334+#include "./Solver.h"
2335+
2336+#endif
2337Index: ../trunk-jpl/src/c/classes/toolkits/Vector.h
2338===================================================================
2339--- ../trunk-jpl/src/c/classes/toolkits/Vector.h (revision 0)
2340+++ ../trunk-jpl/src/c/classes/toolkits/Vector.h (revision 14792)
2341@@ -0,0 +1,357 @@
2342+/*!\file: Vector.h
2343+ * \brief wrapper to vector objects. The goal is to control which API (PETSc,Scalpack, Plapack?)
2344+ * implements our underlying vector format.
2345+ */
2346+
2347+#ifndef _VECTOR_H_
2348+#define _VECTOR_H_
2349+
2350+/*Headers:*/
2351+/*{{{*/
2352+#ifdef HAVE_CONFIG_H
2353+ #include <config.h>
2354+#else
2355+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
2356+#endif
2357+#include <cstring>
2358+#include "../../toolkits/toolkits.h"
2359+#include "../../EnumDefinitions/EnumDefinitions.h"
2360+/*}}}*/
2361+
2362+enum vectortype { PetscVecType, IssmVecType };
2363+
2364+template <class doubletype>
2365+class Vector{
2366+
2367+ public:
2368+
2369+ int type;
2370+ #ifdef _HAVE_PETSC_
2371+ PetscVec* pvector;
2372+ #endif
2373+ IssmVec<doubletype>* ivector;
2374+
2375+ /*Vector constructors, destructors */
2376+ Vector(){ /*{{{*/
2377+
2378+ InitCheckAndSetType();
2379+ }
2380+ /*}}}*/
2381+ Vector(int M,bool fromlocalsize=false){ /*{{{*/
2382+
2383+ InitCheckAndSetType();
2384+
2385+ if(type==PetscVecType){
2386+ #ifdef _HAVE_PETSC_
2387+ this->pvector=new PetscVec(M,fromlocalsize);
2388+ #endif
2389+ }
2390+ else this->ivector=new IssmVec<doubletype>(M,fromlocalsize);
2391+
2392+ }
2393+ /*}}}*/
2394+ Vector(int m,int M){ /*{{{*/
2395+
2396+ InitCheckAndSetType();
2397+
2398+ if(type==PetscVecType){
2399+ #ifdef _HAVE_PETSC_
2400+ this->pvector=new PetscVec(m,M);
2401+ #endif
2402+ }
2403+ else this->ivector=new IssmVec<doubletype>(m,M);
2404+ }
2405+ /*}}}*/
2406+ Vector(doubletype* serial_vec,int M){ /*{{{*/
2407+
2408+ InitCheckAndSetType();
2409+
2410+ if(type==PetscVecType){
2411+ #ifdef _HAVE_PETSC_
2412+ this->pvector=new PetscVec(serial_vec,M);
2413+ #endif
2414+ }
2415+ else this->ivector=new IssmVec<doubletype>(serial_vec,M);
2416+ }
2417+ /*}}}*/
2418+ ~Vector(){ /*{{{*/
2419+
2420+ if(type==PetscVecType){
2421+ #ifdef _HAVE_PETSC_
2422+ delete this->pvector;
2423+ #endif
2424+ }
2425+ else delete this->ivector;
2426+ }
2427+ /*}}}*/
2428+ #ifdef _HAVE_PETSC_
2429+ Vector(Vec petsc_vector){ /*{{{*/
2430+
2431+ this->type=PetscVecType;
2432+ this->ivector=NULL;
2433+ this->pvector=new PetscVec(petsc_vector);
2434+
2435+ }
2436+ /*}}}*/
2437+ #endif
2438+ void InitCheckAndSetType(void){ /*{{{*/
2439+
2440+ char* toolkittype=NULL;
2441+
2442+ #ifdef _HAVE_PETSC_
2443+ pvector=NULL;
2444+ #endif
2445+ ivector=NULL;
2446+
2447+ /*retrieve toolkittype: */
2448+ toolkittype=ToolkitOptions::GetToolkitType();
2449+
2450+ /*set vector type: */
2451+ if (strcmp(toolkittype,"petsc")==0){
2452+ #ifdef _HAVE_PETSC_
2453+ type=PetscVecType;
2454+ #else
2455+ _error_("cannot create petsc vector without PETSC compiled!");
2456+ #endif
2457+ }
2458+ else if(strcmp(toolkittype,"issm")==0){
2459+ /*let this choice stand:*/
2460+ type=IssmVecType;
2461+ }
2462+ else _error_("unknow toolkit type ");
2463+
2464+ /*Free ressources: */
2465+ xDelete<char>(toolkittype);
2466+ }
2467+ /*}}}*/
2468+
2469+ /*Vector specific routines*/
2470+ /*FUNCTION Echo{{{*/
2471+ void Echo(void){_assert_(this);
2472+
2473+ if(type==PetscVecType){
2474+ #ifdef _HAVE_PETSC_
2475+ this->pvector->Echo();
2476+ #endif
2477+ }
2478+ else this->ivector->Echo();
2479+
2480+ }
2481+ /*}}}*/
2482+ /*FUNCTION Assemble{{{*/
2483+ void Assemble(void){_assert_(this);
2484+
2485+ if(type==PetscVecType){
2486+ #ifdef _HAVE_PETSC_
2487+ this->pvector->Assemble();
2488+ #endif
2489+ }
2490+ else this->ivector->Assemble();
2491+
2492+ }
2493+ /*}}}*/
2494+ /*FUNCTION SetValues{{{*/
2495+ void SetValues(int ssize, int* list, doubletype* values, InsMode mode){ _assert_(this);
2496+ if(type==PetscVecType){
2497+ #ifdef _HAVE_PETSC_
2498+ this->pvector->SetValues(ssize,list,values,mode);
2499+ #endif
2500+ }
2501+ else this->ivector->SetValues(ssize,list,values,mode);
2502+
2503+ }
2504+ /*}}}*/
2505+ /*FUNCTION SetValue{{{*/
2506+ void SetValue(int dof, doubletype value, InsMode mode){_assert_(this);
2507+
2508+ if(type==PetscVecType){
2509+ #ifdef _HAVE_PETSC_
2510+ this->pvector->SetValue(dof,value,mode);
2511+ #endif
2512+ }
2513+ else this->ivector->SetValue(dof,value,mode);
2514+
2515+ }
2516+ /*}}}*/
2517+ /*FUNCTION GetValue{{{*/
2518+ void GetValue(doubletype* pvalue,int dof){_assert_(this);
2519+
2520+ if(type==PetscVecType){
2521+ #ifdef _HAVE_PETSC_
2522+ this->pvector->GetValue(pvalue,dof);
2523+ #endif
2524+ }
2525+ else this->ivector->GetValue(pvalue,dof);
2526+
2527+ }
2528+ /*}}}*/
2529+ /*FUNCTION GetSize{{{*/
2530+ void GetSize(int* pM){_assert_(this);
2531+
2532+ if(type==PetscVecType){
2533+ #ifdef _HAVE_PETSC_
2534+ this->pvector->GetSize(pM);
2535+ #endif
2536+ }
2537+ else this->ivector->GetSize(pM);
2538+
2539+ }
2540+ /*}}}*/
2541+ /*FUNCTION IsEmpty{{{*/
2542+ bool IsEmpty(void){
2543+ int M;
2544+
2545+ _assert_(this);
2546+ this->GetSize(&M);
2547+
2548+ if(M==0)
2549+ return true;
2550+ else
2551+ return false;
2552+ }
2553+ /*}}}*/
2554+ /*FUNCTION GetLocalSize{{{*/
2555+ void GetLocalSize(int* pM){_assert_(this);
2556+
2557+ if(type==PetscVecType){
2558+ #ifdef _HAVE_PETSC_
2559+ this->pvector->GetLocalSize(pM);
2560+ #endif
2561+ }
2562+ else this->ivector->GetLocalSize(pM);
2563+
2564+ }
2565+ /*}}}*/
2566+ /*FUNCTION Duplicate{{{*/
2567+ Vector<doubletype>* Duplicate(void){_assert_(this);
2568+
2569+ Vector<doubletype>* output=NULL;
2570+
2571+ output=new Vector<doubletype>();
2572+
2573+ if(type==PetscVecType){
2574+ #ifdef _HAVE_PETSC_
2575+ output->pvector=this->pvector->Duplicate();
2576+ #endif
2577+ }
2578+ else output->ivector=this->ivector->Duplicate();
2579+
2580+ return output;
2581+
2582+ }
2583+ /*}}}*/
2584+ /*FUNCTION Set{{{*/
2585+ void Set(doubletype value){_assert_(this);
2586+
2587+ if(type==PetscVecType){
2588+ #ifdef _HAVE_PETSC_
2589+ this->pvector->Set(value);
2590+ #endif
2591+ }
2592+ else this->ivector->Set(value);
2593+
2594+ }
2595+ /*}}}*/
2596+ /*FUNCTION AXPY{{{*/
2597+ void AXPY(Vector* X, doubletype a){_assert_(this);
2598+
2599+ if(type==PetscVecType){
2600+ #ifdef _HAVE_PETSC_
2601+ this->pvector->AXPY(X->pvector,a);
2602+ #endif
2603+ }
2604+ else this->ivector->AXPY(X->ivector,a);
2605+
2606+ }
2607+ /*}}}*/
2608+ /*FUNCTION AYPX{{{*/
2609+ void AYPX(Vector* X, doubletype a){_assert_(this);
2610+
2611+ if(type==PetscVecType){
2612+ #ifdef _HAVE_PETSC_
2613+ this->pvector->AYPX(X->pvector,a);
2614+ #endif
2615+ }
2616+ else this->ivector->AYPX(X->ivector,a);
2617+ }
2618+ /*}}}*/
2619+ /*FUNCTION ToMPISerial{{{*/
2620+ doubletype* ToMPISerial(void){
2621+
2622+ doubletype* vec_serial=NULL;
2623+
2624+ _assert_(this);
2625+ if(type==PetscVecType){
2626+ #ifdef _HAVE_PETSC_
2627+ vec_serial=this->pvector->ToMPISerial();
2628+ #endif
2629+ }
2630+ else vec_serial=this->ivector->ToMPISerial();
2631+
2632+ return vec_serial;
2633+
2634+ }
2635+ /*}}}*/
2636+ /*FUNCTION Copy{{{*/
2637+ void Copy(Vector* to){_assert_(this);
2638+
2639+ if(type==PetscVecType){
2640+ #ifdef _HAVE_PETSC_
2641+ this->pvector->Copy(to->pvector);
2642+ #endif
2643+ }
2644+ else this->ivector->Copy(to->ivector);
2645+ }
2646+ /*}}}*/
2647+ /*FUNCTION Norm{{{*/
2648+ doubletype Norm(NormMode norm_type){_assert_(this);
2649+
2650+ doubletype norm=0;
2651+
2652+ if(type==PetscVecType){
2653+ #ifdef _HAVE_PETSC_
2654+ norm=this->pvector->Norm(norm_type);
2655+ #endif
2656+ }
2657+ else norm=this->ivector->Norm(norm_type);
2658+ return norm;
2659+ }
2660+ /*}}}*/
2661+ /*FUNCTION Scale{{{*/
2662+ void Scale(doubletype scale_factor){_assert_(this);
2663+
2664+ if(type==PetscVecType){
2665+ #ifdef _HAVE_PETSC_
2666+ this->pvector->Scale(scale_factor);
2667+ #endif
2668+ }
2669+ else this->ivector->Scale(scale_factor);
2670+ }
2671+ /*}}}*/
2672+ /*FUNCTION Dot{{{*/
2673+ doubletype Dot(Vector* vector){_assert_(this);
2674+
2675+ doubletype dot;
2676+
2677+ if(type==PetscVecType){
2678+ #ifdef _HAVE_PETSC_
2679+ dot=this->pvector->Dot(vector->pvector);
2680+ #endif
2681+ }
2682+ else dot=this->ivector->Dot(vector->ivector);
2683+ return dot;
2684+ }
2685+ /*}}}*/
2686+ /*FUNCTION PointwiseDivide{{{*/
2687+ void PointwiseDivide(Vector* x,Vector* y){_assert_(this);
2688+
2689+ if(type==PetscVecType){
2690+ #ifdef _HAVE_PETSC_
2691+ this->pvector->PointwiseDivide(x->pvector,y->pvector);
2692+ #endif
2693+ }
2694+ else this->ivector->PointwiseDivide(x->ivector,y->ivector);
2695+ }
2696+ /*}}}*/
2697+};
2698+#endif //#ifndef _VECTOR_H_
2699Index: ../trunk-jpl/src/c/classes/toolkits/Matrix.h
2700===================================================================
2701--- ../trunk-jpl/src/c/classes/toolkits/Matrix.h (revision 0)
2702+++ ../trunk-jpl/src/c/classes/toolkits/Matrix.h (revision 14792)
2703@@ -0,0 +1,329 @@
2704+/*!\file: Matrix.h
2705+ * \brief wrapper to matrix objects. The goal is to control which API (PETSc,Scalpack, Plapack?)
2706+ * implements our underlying matrix format.
2707+ */
2708+
2709+#ifndef _MATRIX_H_
2710+#define _MATRIX_H_
2711+
2712+/*Headers:*/
2713+/*{{{*/
2714+#ifdef HAVE_CONFIG_H
2715+ #include <config.h>
2716+#else
2717+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
2718+#endif
2719+#include <cstring>
2720+#include "../../toolkits/toolkits.h"
2721+#include "../../EnumDefinitions/EnumDefinitions.h"
2722+/*}}}*/
2723+
2724+enum matrixtype { PetscMatType, IssmMatType };
2725+
2726+template <class doubletype> class Vector;
2727+
2728+template <class doubletype>
2729+class Matrix{
2730+
2731+ public:
2732+
2733+ int type;
2734+ #ifdef _HAVE_PETSC_
2735+ PetscMat *pmatrix;
2736+ #endif
2737+ IssmMat<doubletype> *imatrix;
2738+
2739+ /*Matrix constructors, destructors*/
2740+ /*FUNCTION Matrix(){{{*/
2741+ Matrix(){
2742+ InitCheckAndSetType();
2743+ }
2744+ /*}}}*/
2745+ /*FUNCTION Matrix(int M,int N){{{*/
2746+ Matrix(int M,int N){
2747+
2748+ InitCheckAndSetType();
2749+
2750+ if(type==PetscMatType){
2751+ #ifdef _HAVE_PETSC_
2752+ this->pmatrix=new PetscMat(M,N);
2753+ #endif
2754+ }
2755+ else{
2756+ this->imatrix=new IssmMat<doubletype>(M,N);
2757+ }
2758+
2759+ }
2760+ /*}}}*/
2761+ /*FUNCTION Matrix(int m,int n,int M,int N,int* d_nnz,int* o_nnz){{{*/
2762+ Matrix(int m,int n,int M,int N,int* d_nnz,int* o_nnz){
2763+
2764+ InitCheckAndSetType();
2765+
2766+ if(type==PetscMatType){
2767+ #ifdef _HAVE_PETSC_
2768+ this->pmatrix=new PetscMat(m,n,M,N,d_nnz,o_nnz);
2769+ #endif
2770+ }
2771+ else{
2772+ this->imatrix=new IssmMat<doubletype>(m,n,M,N,d_nnz,o_nnz);
2773+ }
2774+
2775+ }
2776+ /*}}}*/
2777+ /*FUNCTION Matrix(int M,int N,IssmDouble sparsity){{{*/
2778+ Matrix(int M,int N,double sparsity){
2779+
2780+ InitCheckAndSetType();
2781+
2782+ if(type==PetscMatType){
2783+ #ifdef _HAVE_PETSC_
2784+ this->pmatrix=new PetscMat(M,N,sparsity);
2785+ #endif
2786+ }
2787+ else{
2788+ this->imatrix=new IssmMat<doubletype>(M,N,sparsity);
2789+ }
2790+ }
2791+ /*}}}*/
2792+ /*FUNCTION Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity){{{*/
2793+ Matrix(IssmPDouble* serial_mat,int M,int N,IssmPDouble sparsity){
2794+
2795+ InitCheckAndSetType();
2796+
2797+ if(type==PetscMatType){
2798+ #ifdef _HAVE_PETSC_
2799+ this->pmatrix=new PetscMat(serial_mat,M,N,sparsity);
2800+ #endif
2801+ }
2802+ else{
2803+ this->imatrix=new IssmMat<doubletype>(serial_mat,M,N,sparsity);
2804+ }
2805+
2806+ }
2807+ /*}}}*/
2808+ /*FUNCTION Matrix(int M,int N,int connectivity,int numberofdofspernode){{{*/
2809+ Matrix(int M,int N,int connectivity,int numberofdofspernode){
2810+
2811+ InitCheckAndSetType();
2812+
2813+ if(type==PetscMatType){
2814+ #ifdef _HAVE_PETSC_
2815+ this->pmatrix=new PetscMat(M,N,connectivity,numberofdofspernode);
2816+ #endif
2817+ }
2818+ else{
2819+ this->imatrix=new IssmMat<doubletype>(M,N,connectivity,numberofdofspernode);
2820+ }
2821+
2822+ }
2823+ /*}}}*/
2824+ /*FUNCTION ~Matrix(){{{*/
2825+ ~Matrix(){
2826+
2827+ if(type==PetscMatType){
2828+ #ifdef _HAVE_PETSC_
2829+ delete this->pmatrix;
2830+ #endif
2831+ }
2832+ else delete this->imatrix;
2833+
2834+ }
2835+ /*}}}*/
2836+ /*FUNCTION InitCheckAndSetType(){{{*/
2837+ void InitCheckAndSetType(void){
2838+
2839+ char* toolkittype=NULL;
2840+
2841+ #ifdef _HAVE_PETSC_
2842+ pmatrix=NULL;
2843+ #endif
2844+ imatrix=NULL;
2845+
2846+ /*retrieve toolkittype: */
2847+ toolkittype=ToolkitOptions::GetToolkitType();
2848+
2849+ /*set matrix type: */
2850+ if (strcmp(toolkittype,"petsc")==0){
2851+ #ifdef _HAVE_PETSC_
2852+ type=PetscMatType;
2853+ #else
2854+ _error_("cannot create petsc matrix without PETSC compiled!");
2855+ #endif
2856+ }
2857+ else if(strcmp(toolkittype,"issm")==0){
2858+ /*let this choice stand:*/
2859+ type=IssmMatType;
2860+ }
2861+ else _error_("unknow toolkit type ");
2862+
2863+ /*Free ressources: */
2864+ xDelete<char>(toolkittype);
2865+ }
2866+ /*}}}*/
2867+
2868+ /*Matrix specific routines:*/
2869+ /*FUNCTION Echo{{{*/
2870+ void Echo(void){
2871+ _assert_(this);
2872+
2873+ if(type==PetscMatType){
2874+ #ifdef _HAVE_PETSC_
2875+ this->pmatrix->Echo();
2876+ #endif
2877+ }
2878+ else{
2879+ this->imatrix->Echo();
2880+ }
2881+
2882+ }
2883+ /*}}}*/
2884+ /*FUNCTION AllocationInfo{{{*/
2885+ void AllocationInfo(void){
2886+ _assert_(this);
2887+ if(type==PetscMatType){
2888+ #ifdef _HAVE_PETSC_
2889+ this->pmatrix->AllocationInfo();
2890+ #endif
2891+ }
2892+ else{
2893+ this->imatrix->AllocationInfo();
2894+ }
2895+ }/*}}}*/
2896+ /*FUNCTION Assemble{{{*/
2897+ void Assemble(void){
2898+
2899+ if(type==PetscMatType){
2900+ #ifdef _HAVE_PETSC_
2901+ this->pmatrix->Assemble();
2902+ #endif
2903+ }
2904+ else{
2905+ this->imatrix->Assemble();
2906+ }
2907+ }
2908+ /*}}}*/
2909+ /*FUNCTION Norm{{{*/
2910+ IssmDouble Norm(NormMode norm_type){
2911+
2912+ IssmDouble norm=0;
2913+
2914+ if(type==PetscMatType){
2915+ #ifdef _HAVE_PETSC_
2916+ norm=this->pmatrix->Norm(norm_type);
2917+ #endif
2918+ }
2919+ else{
2920+ norm=this->imatrix->Norm(norm_type);
2921+ }
2922+
2923+ return norm;
2924+ }
2925+ /*}}}*/
2926+ /*FUNCTION GetSize{{{*/
2927+ void GetSize(int* pM,int* pN){
2928+
2929+ if(type==PetscMatType){
2930+ #ifdef _HAVE_PETSC_
2931+ this->pmatrix->GetSize(pM,pN);
2932+ #endif
2933+ }
2934+ else{
2935+ this->imatrix->GetSize(pM,pN);
2936+ }
2937+
2938+ }
2939+ /*}}}*/
2940+ /*FUNCTION GetLocalSize{{{*/
2941+ void GetLocalSize(int* pM,int* pN){
2942+
2943+ if(type==PetscMatType){
2944+ #ifdef _HAVE_PETSC_
2945+ this->pmatrix->GetLocalSize(pM,pN);
2946+ #endif
2947+ }
2948+ else{
2949+ this->imatrix->GetLocalSize(pM,pN);
2950+ }
2951+
2952+ }
2953+ /*}}}*/
2954+ /*FUNCTION MatMult{{{*/
2955+ void MatMult(Vector<doubletype>* X,Vector<doubletype>* AX){
2956+
2957+ if(type==PetscMatType){
2958+ #ifdef _HAVE_PETSC_
2959+ this->pmatrix->MatMult(X->pvector,AX->pvector);
2960+ #endif
2961+ }
2962+ else{
2963+ this->imatrix->MatMult(X->ivector,AX->ivector);
2964+ }
2965+
2966+ }
2967+ /*}}}*/
2968+ /*FUNCTION Duplicate{{{*/
2969+ Matrix<doubletype>* Duplicate(void){
2970+
2971+ Matrix<doubletype>* output=new Matrix<doubletype>();
2972+
2973+ if(type==PetscMatType){
2974+ #ifdef _HAVE_PETSC_
2975+ output->pmatrix=this->pmatrix->Duplicate();
2976+ #endif
2977+ }
2978+ else{
2979+ output->imatrix=this->imatrix->Duplicate();
2980+ }
2981+
2982+ return output;
2983+ }
2984+ /*}}}*/
2985+ /*FUNCTION ToSerial{{{*/
2986+ doubletype* ToSerial(void){
2987+
2988+ doubletype* output=NULL;
2989+
2990+ if(type==PetscMatType){
2991+ #ifdef _HAVE_PETSC_
2992+ output=this->pmatrix->ToSerial();
2993+ #endif
2994+ }
2995+ else{
2996+ output=this->imatrix->ToSerial();
2997+ }
2998+
2999+ return output;
3000+ }
3001+ /*}}}*/
3002+ /*FUNCTION SetValues{{{*/
3003+ void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){
3004+
3005+ if(type==PetscMatType){
3006+ #ifdef _HAVE_PETSC_
3007+ this->pmatrix->SetValues(m,idxm,n,idxn,values,mode);
3008+ #endif
3009+ }
3010+ else{
3011+ this->imatrix->SetValues(m,idxm,n,idxn,values,mode);
3012+ }
3013+ }
3014+ /*}}}*/
3015+ /*FUNCTION Convert{{{*/
3016+ void Convert(MatrixType newtype){
3017+
3018+ if(type==PetscMatType){
3019+ #ifdef _HAVE_PETSC_
3020+ this->pmatrix->Convert(newtype);
3021+ #endif
3022+ }
3023+ else{
3024+ this->imatrix->Convert(newtype);
3025+ }
3026+
3027+ }
3028+ /*}}}*/
3029+
3030+};
3031+
3032+#endif //#ifndef _MATRIX_H_
Note: See TracBrowser for help on using the repository browser.