![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
source:
issm/oecreview/Archive/14312-15392/ISSM-14791-14792.diff
| Last change on this file was 15393, checked in by , 12 years ago | |
|---|---|
| File size: 81.7 KB | |
-
../trunk-jpl/src/c/modules/Solverx/SolverxPetsc.cpp
1 /*!\file SolverxPetsc2 * \brief Petsc implementation of solver3 */4 5 #include "./Solverx.h"6 #include "../../shared/shared.h"7 #include "../../include/include.h"8 #include "../../io/io.h"9 10 #ifdef HAVE_CONFIG_H11 #include <config.h>12 #else13 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"14 #endif15 16 void SolverxPetsc(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters){17 18 PetscVec* uf=new PetscVec();19 20 Vec uf0_vector = NULL;21 Vec df_vector = NULL;22 23 if(uf0) uf0_vector = uf0->vector;24 if(df) df_vector = df->vector;25 26 SolverxPetsc(&uf->vector, Kff->matrix, pf->vector, uf0_vector, df_vector, parameters);27 28 *puf=uf;29 30 }31 void SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){32 33 /*Output: */34 Vec uf = NULL;35 36 /*Intermediary: */37 int local_m,local_n,global_m,global_n;38 39 /*Solver */40 KSP ksp = NULL;41 PC pc = NULL;42 int iteration_number;43 int solver_type;44 bool fromlocalsize = true;45 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)46 PetscTruth flag,flg;47 #else48 PetscBool flag,flg;49 #endif50 51 /*Stokes: */52 IS isv=NULL;53 IS isp=NULL;54 55 #if _PETSC_MAJOR_ >= 356 char ksp_type[50];57 #endif58 59 /*Display message*/60 if(VerboseModule()) _pprintLine_(" Solving");61 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)62 if(VerboseSolver())PetscOptionsPrint(stdout);63 #else64 if(VerboseSolver())PetscOptionsView(PETSC_VIEWER_STDOUT_WORLD);65 #endif66 67 /*First, check that f-set is not NULL, i.e. model is fully constrained:*/68 _assert_(Kff);69 MatGetSize(Kff,&global_m,&global_n); _assert_(global_m==global_n);70 if(!global_n){71 *puf=NewVec(0,IssmComm::GetComm()); return;72 }73 74 /*Initial guess */75 /*Now, check that we are not giving an initial guess to the solver, if we are running a direct solver: */76 #if _PETSC_MAJOR_ >= 377 PetscOptionsGetString(PETSC_NULL,"-ksp_type",ksp_type,49,&flg);78 if (strcmp(ksp_type,"preonly")==0)uf0=NULL;79 #endif80 81 /*If initial guess for the solution exists, use it to create uf, otherwise,82 * duplicate the right hand side so that the solution vector has the same structure*/83 if(uf0){84 VecDuplicate(uf0,&uf); VecCopy(uf0,uf);85 }86 else{87 MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,IssmComm::GetComm(),fromlocalsize);88 }89 90 /*Process petsc options to see if we are using special types of external solvers*/91 PetscOptionsDetermineSolverType(&solver_type);92 93 /*Check the solver is available*/94 if(solver_type==MUMPSPACKAGE_LU || solver_type==MUMPSPACKAGE_CHOL){95 #if _PETSC_MAJOR_ >=396 #ifndef _HAVE_MUMPS_97 _error_("requested MUMPS solver, which was not compiled into ISSM!\n");98 #endif99 #endif100 }101 102 /*Prepare solver*/103 KSPCreate(IssmComm::GetComm(),&ksp);104 KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN);105 KSPSetFromOptions(ksp);106 107 #if _PETSC_MAJOR_==3108 /*Specific solver?: */109 KSPGetPC(ksp,&pc);110 if (solver_type==MUMPSPACKAGE_LU){111 #if _PETSC_MINOR_==1112 PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS);113 #else114 PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);115 #endif116 }117 118 /*Stokes: */119 if (solver_type==StokesSolverEnum){120 /*Make indices out of doftypes: */121 if(!df)_error_("need doftypes for Stokes solver!\n");122 DofTypesToIndexSet(&isv,&isp,df,StokesSolverEnum);123 124 /*Set field splits: */125 KSPGetPC(ksp,&pc);126 #if _PETSC_MINOR_==1127 PCFieldSplitSetIS(pc,isv);128 PCFieldSplitSetIS(pc,isp);129 #else130 PCFieldSplitSetIS(pc,PETSC_NULL,isv);131 PCFieldSplitSetIS(pc,PETSC_NULL,isp);132 #endif133 134 }135 #endif136 137 /*If there is an initial guess for the solution, use it138 * except if we are using the MUMPS direct solver139 * where any initial guess will crash Petsc*/140 if (uf0){141 if((solver_type!=MUMPSPACKAGE_LU) && (solver_type!=MUMPSPACKAGE_CHOL) && (solver_type!=SPOOLESPACKAGE_LU) && (solver_type!=SPOOLESPACKAGE_CHOL) && (solver_type!=SUPERLUDISTPACKAGE)){142 KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);143 }144 }145 146 /*Solve: */147 if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);148 KSPSolve(ksp,pf,uf);149 150 /*Check convergence*/151 KSPGetIterationNumber(ksp,&iteration_number);152 if (iteration_number<0) _error_("Solver diverged at iteration number: " << -iteration_number);153 154 /*Free resources:*/155 KSPFree(&ksp);156 157 /*Assign output pointers:*/158 *puf=uf;159 } -
../trunk-jpl/src/c/modules/Solverx/DofTypesToIndexSet.cpp
1 /*!\file Solverx2 * \brief solver3 */4 5 #include "./Solverx.h"6 #include "../../shared/shared.h"7 #include "../../include/include.h"8 9 #ifdef HAVE_CONFIG_H10 #include <config.h>11 #else12 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"13 #endif14 15 void DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum){16 17 /*output: */18 IS isv=NULL;19 IS isp=NULL;20 21 int start,end;22 IssmDouble* df_local=NULL;23 int df_local_size;24 int i;25 26 int* pressure_indices=NULL;27 int* velocity_indices=NULL;28 int pressure_num=0;29 int velocity_num=0;30 int pressure_count=0;31 int velocity_count=0;32 33 if(typeenum==StokesSolverEnum){34 35 /*Ok, recover doftypes vector values and indices: */36 VecGetOwnershipRange(df,&start,&end);37 VecGetLocalSize(df,&df_local_size);38 VecGetArray(df,&df_local);39 40 pressure_num=0;41 velocity_num=0;42 for(i=0;i<df_local_size;i++){43 if (df_local[i]==PressureEnum)pressure_num++;44 else velocity_num++;45 }46 47 /*Allocate indices: */48 if(pressure_num)pressure_indices=xNew<int>(pressure_num);49 if(velocity_num)velocity_indices=xNew<int>(velocity_num);50 51 pressure_count=0;52 velocity_count=0;53 for(i=0;i<df_local_size;i++){54 if (df_local[i]==PressureEnum){55 pressure_indices[pressure_count]=start+i;56 pressure_count++;57 }58 if (df_local[i]==VelocityEnum){59 velocity_indices[velocity_count]=start+i;60 velocity_count++;61 }62 }63 VecRestoreArray(df,&df_local);64 65 /*Create indices sets: */66 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)67 ISCreateGeneral(IssmComm::GetComm(),pressure_num,pressure_indices,&isp);68 ISCreateGeneral(IssmComm::GetComm(),velocity_num,velocity_indices,&isv);69 #else70 ISCreateGeneral(IssmComm::GetComm(),pressure_num,pressure_indices,PETSC_COPY_VALUES,&isp);71 ISCreateGeneral(IssmComm::GetComm(),velocity_num,velocity_indices,PETSC_COPY_VALUES,&isv);72 #endif73 }74 75 /*Free ressources:*/76 xDelete<int>(pressure_indices);77 xDelete<int>(velocity_indices);78 79 /*Assign output pointers:*/80 *pisv=isv;81 *pisp=isp;82 } -
../trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp
1 /*!\file SolverxSeq2 * \brief implementation of sequential solver using the GSL librarie3 */4 5 #ifdef HAVE_CONFIG_H6 #include <config.h>7 #else8 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"9 #endif10 #include <cstring>11 12 #include "./Solverx.h"13 #include "../../shared/shared.h"14 #include "../../classes/classes.h"15 #include "../../include/include.h"16 #include "../../io/io.h"17 18 #ifdef _HAVE_GSL_19 #include <gsl/gsl_linalg.h>20 #endif21 22 void SolverxSeq(IssmVec<IssmDouble>** pout,IssmMat<IssmDouble>* Kffin, IssmVec<IssmDouble>* pfin, Parameters* parameters){/*{{{*/23 24 /*First off, we assume that the type of IssmVec is IssmSeqVec and the type of IssmMat is IssmDenseMat. So downcast: */25 IssmSeqVec<IssmDouble>* pf=(IssmSeqVec<IssmDouble>*)pfin->vector;26 IssmDenseMat<IssmDouble>* Kff=(IssmDenseMat<IssmDouble>*)Kffin->matrix;27 28 #ifdef _HAVE_GSL_29 /*Intermediary: */30 int M,N,N2;31 IssmSeqVec<IssmDouble> *uf = NULL;32 33 Kff->GetSize(&M,&N);34 pf->GetSize(&N2);35 36 if(N!=N2)_error_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !");37 if(M!=N)_error_("Stiffness matrix should be square!");38 #ifdef _HAVE_ADOLC_39 ensureContiguousLocations(N);40 #endif41 IssmDouble *x = xNew<IssmDouble>(N);42 #ifdef _HAVE_ADOLC_43 SolverxSeq(x,Kff->matrix,pf->vector,N,parameters);44 #else45 SolverxSeq(x,Kff->matrix,pf->vector,N);46 #endif47 48 uf=new IssmSeqVec<IssmDouble>(x,N);49 xDelete(x);50 51 /*Assign output pointers:*/52 IssmVec<IssmDouble>* out=new IssmVec<IssmDouble>();53 out->vector=(IssmAbsVec<IssmDouble>*)uf;54 *pout=out;55 56 #else57 _error_("GSL support not compiled in!");58 #endif59 60 }/*}}}*/61 void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/62 63 /*Allocate output*/64 double* X = xNew<double>(n);65 66 /*Solve*/67 SolverxSeq(X,A,B,n);68 69 /*Assign output pointer*/70 *pX=X;71 }72 /*}}}*/73 74 #ifdef _HAVE_ADOLC_75 int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y){ /*{{{*/76 SolverxSeq(y,x, x+m*m, m); // x is where the matrix starts, x+m*m is where the right-hand side starts77 return 0;78 } /*}}}*/79 int EDF_fos_forward_for_solverx(int n, IssmPDouble *inVal, IssmPDouble *inDeriv, int m, IssmPDouble *outVal, IssmPDouble *outDeriv) { /*{{{*/80 #ifdef _HAVE_GSL_81 // for (int i=0; i<m*m; ++i) std::cout << "EDF_fos_forward_for_solverx A["<< i << "]=" << inVal[i] << std::endl;82 // for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx b["<< i << "]=" << inVal[i+m*m] << std::endl;83 // the matrix will be modified by LU decomposition. Use gsl_A copy84 double* Acopy = xNew<double>(m*m);85 xMemCpy(Acopy,inVal,m*m);86 /*Initialize gsl matrices and vectors: */87 gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);88 gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m89 gsl_permutation *perm_p = gsl_permutation_alloc (m);90 int signPerm;91 // factorize92 gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);93 gsl_vector *gsl_x_p = gsl_vector_alloc (m);94 // solve for the value95 gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);96 /*Copy result*/97 xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);98 gsl_vector_free(gsl_x_p);99 // for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx x["<< i << "]=" << outVal[i] << std::endl;100 // solve for the derivatives acc. to A * dx = r with r=db - dA * x101 // compute the RHS102 double* r=xNew<double>(m);103 for (int i=0; i<m; i++) {104 r[i]=inDeriv[m*m+i]; // this is db[i]105 for (int j=0;j<m; j++) {106 r[i]-=inDeriv[i*m+j]*outVal[j]; // this is dA[i][j]*x[j]107 }108 }109 gsl_vector_view gsl_r=gsl_vector_view_array(r,m);110 gsl_vector *gsl_dx_p = gsl_vector_alloc(m);111 gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);112 xMemCpy(outDeriv,gsl_vector_ptr(gsl_dx_p,0),m);113 gsl_vector_free(gsl_dx_p);114 xDelete(r);115 gsl_permutation_free(perm_p);116 xDelete(Acopy);117 #endif118 return 0;119 } /*}}}*/120 int EDF_fov_forward_for_solverx(int n, IssmPDouble *inVal, int directionCount, IssmPDouble **inDeriv, int m, IssmPDouble *outVal, IssmPDouble **outDeriv) { /*{{{*/121 #ifdef _HAVE_GSL_122 // the matrix will be modified by LU decomposition. Use gsl_A copy123 double* Acopy = xNew<double>(m*m);124 xMemCpy(Acopy,inVal,m*m);125 /*Initialize gsl matrices and vectors: */126 gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);127 gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m128 gsl_permutation *perm_p = gsl_permutation_alloc (m);129 int signPerm;130 // factorize131 gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);132 gsl_vector *gsl_x_p = gsl_vector_alloc (m);133 // solve for the value134 gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);135 /*Copy result*/136 xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);137 gsl_vector_free(gsl_x_p);138 // solve for the derivatives acc. to A * dx = r with r=db - dA * x139 double* r=xNew<double>(m);140 gsl_vector *gsl_dx_p = gsl_vector_alloc(m);141 for (int dir=0;dir<directionCount;++dir) {142 // compute the RHS143 for (int i=0; i<m; i++) {144 r[i]=inDeriv[m*m+i][dir]; // this is db[i]145 for (int j=0;j<m; j++) {146 r[i]-=inDeriv[i*m+j][dir]*outVal[j]; // this is dA[i][j]*x[j]147 }148 }149 gsl_vector_view gsl_r=gsl_vector_view_array(r,m);150 gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);151 // reuse r152 xMemCpy(r,gsl_vector_ptr(gsl_dx_p,0),m);153 for (int i=0; i<m; i++) {154 outDeriv[i][dir]=r[i];155 }156 }157 gsl_vector_free(gsl_dx_p);158 xDelete(r);159 gsl_permutation_free(perm_p);160 xDelete(Acopy);161 #endif162 return 0;163 }164 /*}}}*/165 int EDF_fos_reverse_for_solverx(int m, double *dp_U, int n, double *dp_Z, double* dp_x, double* dp_y) { /*{{{*/166 // copy to transpose the matrix167 double* transposed=xNew<double>(m*m);168 for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j];169 gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m);170 // the adjoint of the solution is our right-hand side171 gsl_vector_view x_bar=gsl_vector_view_array(dp_U,m);172 // the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute:173 gsl_vector_view b_bar=gsl_vector_view_array(dp_Z+m*m,m);174 gsl_permutation *perm_p = gsl_permutation_alloc (m);175 int permSign;176 gsl_linalg_LU_decomp (&aTransposed.matrix, perm_p, &permSign);177 gsl_linalg_LU_solve (&aTransposed.matrix, perm_p, &x_bar.vector, &b_bar.vector);178 // now do the adjoint of the matrix179 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];180 gsl_permutation_free(perm_p);181 xDelete(transposed);182 return 0;183 }184 /*}}}*/185 int EDF_fov_reverse_for_solverx(int m, int p, double **dpp_U, int n, double **dpp_Z, double* dp_x, double* dp_y) { /*{{{*/186 // copy to transpose the matrix187 double* transposed=xNew<double>(m*m);188 for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j];189 gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m);190 gsl_permutation *perm_p = gsl_permutation_alloc (m);191 int permSign;192 gsl_linalg_LU_decomp (&aTransposed.matrix, perm_p, &permSign);193 for (int weightsRowIndex=0;weightsRowIndex<p;++weightsRowIndex) {194 // the adjoint of the solution is our right-hand side195 gsl_vector_view x_bar=gsl_vector_view_array(dpp_U[weightsRowIndex],m);196 // the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute:197 gsl_vector_view b_bar=gsl_vector_view_array(dpp_Z[weightsRowIndex]+m*m,m);198 gsl_linalg_LU_solve (&aTransposed.matrix, perm_p, &x_bar.vector, &b_bar.vector);199 // now do the adjoint of the matrix200 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];201 }202 gsl_permutation_free(perm_p);203 xDelete(transposed);204 return 0;205 }206 /*}}}*/207 void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/208 // pack inputs to conform to the EDF-prescribed interface209 // ensure a contiguous block of locations:210 ensureContiguousLocations(n*(n+1));211 IssmDouble* adoubleEDFin=xNew<IssmDouble>(n*(n+1)); // packed inputs, i.e. matrix and right hand side212 for(int i=0; i<n*n;i++)adoubleEDFin[i] =A[i]; // pack matrix213 for(int i=0; i<n; i++)adoubleEDFin[i+n*n]=B[i]; // pack the right hand side214 IssmPDouble* pdoubleEDFin=xNew<IssmPDouble>(n*(n+1)); // provide space to transfer inputs during call_ext_fct215 IssmPDouble* pdoubleEDFout=xNew<IssmPDouble>(n); // provide space to transfer outputs during call_ext_fct216 // call the wrapped solver through the registry entry we retrieve from parameters217 call_ext_fct(dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p,218 n*(n+1), pdoubleEDFin, adoubleEDFin,219 n, pdoubleEDFout,X);220 // for(int i=0; i<n; i++) {ADOLC_DUMP_MACRO(X[i]);}221 xDelete(adoubleEDFin);222 xDelete(pdoubleEDFin);223 xDelete(pdoubleEDFout);224 }225 /*}}}*/226 #endif227 void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/228 #ifdef _HAVE_GSL_229 /*GSL Matrices and vectors: */230 int s;231 gsl_matrix_view a;232 gsl_vector_view b,x;233 gsl_permutation *p = NULL;234 // for (int i=0; i<n*n; ++i) std::cout << "SolverxSeq A["<< i << "]=" << A[i] << std::endl;235 // for (int i=0; i<n; ++i) std::cout << "SolverxSeq b["<< i << "]=" << B[i] << std::endl;236 /*A will be modified by LU decomposition. Use copy*/237 double* Acopy = xNew<double>(n*n);238 xMemCpy(Acopy,A,n*n);239 240 /*Initialize gsl matrices and vectors: */241 a = gsl_matrix_view_array (Acopy,n,n);242 b = gsl_vector_view_array (B,n);243 x = gsl_vector_view_array (X,n);244 245 /*Run LU and solve: */246 p = gsl_permutation_alloc (n);247 gsl_linalg_LU_decomp (&a.matrix, p, &s);248 gsl_linalg_LU_solve (&a.matrix, p, &b.vector, &x.vector);249 250 /*Clean up and assign output pointer*/251 xDelete(Acopy);252 gsl_permutation_free(p);253 #endif254 }255 /*}}}*/ -
../trunk-jpl/src/c/modules/Solverx/Solverx.cpp
2 2 * \brief solver 3 3 */ 4 4 5 #include "./Solverx.h"6 #include "../../shared/shared.h"7 #include "../../include/include.h"8 #include "../../io/io.h"9 10 5 #ifdef HAVE_CONFIG_H 11 6 #include <config.h> 12 7 #else 13 8 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" 14 9 #endif 15 10 11 #include "./Solverx.h" 12 #include "../../shared/shared.h" 13 16 14 void Solverx(Vector<IssmDouble>** puf, Matrix<IssmDouble>* Kff, Vector<IssmDouble>* pf, Vector<IssmDouble>* uf0,Vector<IssmDouble>* df, Parameters* parameters){ 17 15 16 /*intermediary: */ 17 Solver<IssmDouble> *solver=NULL; 18 18 19 /*output: */ 19 20 Vector<IssmDouble> *uf=NULL; 20 21 21 /*In debugging mode, check that stiffness matrix and load vectors are not NULL (they can be empty)*/22 _assert_(Kff);23 _assert_(pf);24 25 22 if(VerboseModule()) _pprintLine_(" Solving matrix system"); 26 23 27 /*Initialize vector: */28 uf=new Vector<IssmDouble>();24 /*Initialize solver: */ 25 solver=new Solver<IssmDouble>(Kff,pf,uf0,df,parameters); 29 26 30 /*According to matrix type, use specific solvers: */ 31 switch(Kff->type){ 32 #ifdef _HAVE_PETSC_ 33 case PetscMatType:{ 34 PetscVec* uf0_vector = NULL; 35 PetscVec* df_vector = NULL; 36 if(uf0) uf0_vector = uf0->pvector; 37 if(df) df_vector = df->pvector; 38 SolverxPetsc(&uf->pvector,Kff->pmatrix,pf->pvector,uf0_vector,df_vector,parameters); 39 break;} 40 #endif 41 case IssmMatType:{ 42 SolverxSeq(&uf->ivector,Kff->imatrix,pf->ivector,parameters); 43 break;} 44 default: 45 _error_("Matrix type: " << Kff->type << " not supported yet!"); 46 } 27 /*Solve:*/ 28 uf=solver->Solve(); 47 29 48 30 /*Assign output pointers:*/ 49 31 *puf=uf; -
../trunk-jpl/src/c/modules/Solverx/Solverx.h
11 11 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" 12 12 #endif 13 13 14 #include "../../classes/objects/objects.h" 15 #include "../../toolkits/toolkits.h" 14 #include "../../classes/toolkits/toolkitobjects.h" 16 15 17 16 /* local prototypes: */ 18 17 void Solverx(Vector<IssmDouble>** puf, Matrix<IssmDouble>* Kff, Vector<IssmDouble>* pf, Vector<IssmDouble>* uf0,Vector<IssmDouble>* df, Parameters* parameters); 19 18 20 #ifdef _HAVE_PETSC_21 void SolverxPetsc(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters);22 void SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters);23 void DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum);24 #endif25 26 void SolverxSeq(IssmVec<IssmDouble>** puf,IssmMat<IssmDouble>* Kff, IssmVec<IssmDouble>* pf,Parameters* parameters);27 void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n);28 void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n);29 30 #if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_)31 void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters);32 // call back functions:33 ADOLC_ext_fct EDF_for_solverx;34 ADOLC_ext_fct_fos_forward EDF_fos_forward_for_solverx;35 ADOLC_ext_fct_fos_reverse EDF_fos_reverse_for_solverx;36 ADOLC_ext_fct_fov_forward EDF_fov_forward_for_solverx;37 ADOLC_ext_fct_fov_reverse EDF_fov_reverse_for_solverx;38 #endif39 40 19 #endif /* _SOLVERX_H */ -
../trunk-jpl/src/c/modules/Solverx/CMakeLists.txt
4 4 include_directories(AFTER $ENV{ISSM_DIR}/src/c/modules/Solverx) 5 5 # }}} 6 6 # CORE_SOURCES {{{ 7 set(CORE_SOURCES $ENV{ISSM_DIR}/src/c/modules/Solverx/Solverx.cpp 8 $ENV{ISSM_DIR}/src/c/modules/Solverx/SolverxSeq.cpp PARENT_SCOPE) 7 set(CORE_SOURCES $ENV{ISSM_DIR}/src/c/modules/Solverx/Solverx.cpp PARENT_SCOPE) 9 8 # }}} -
../trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h
21 21 #include "../../shared/MemOps/xMemCpy.h" 22 22 #include "../../shared/Alloc/alloc.h" 23 23 #include "../../include/macros.h" 24 #include "../../io/io.h" 24 25 #ifdef _HAVE_MPI_ 25 26 #include "../mpi/mpiincludes.h" 26 27 #endif … … 302 303 /*}}}*/ 303 304 /*FUNCTION AXPY{{{*/ 304 305 void AXPY(IssmAbsVec<doubletype>* Xin, doubletype a){ 305 printf("AXPY not implemented yet!"); 306 exit(1); 306 307 int i; 308 309 /*Assume X is of the correct type, and downcast: */ 310 IssmMpiVec* X=NULL; 311 312 X=(IssmMpiVec<doubletype>*)Xin; 313 314 /*y=a*x+y where this->vector is y*/ 315 for(i=0;i<this->m;i++)this->vector[i]=a*X->vector[i]+this->vector[i]; 316 317 307 318 } 308 319 /*}}}*/ 309 320 /*FUNCTION AYPX{{{*/ 310 321 void AYPX(IssmAbsVec<doubletype>* Xin, doubletype a){ 311 printf("AYPX not implemented yet!"); 312 exit(1); 322 int i; 323 324 /*Assume X is of the correct type, and downcast: */ 325 IssmMpiVec* X=NULL; 326 327 X=(IssmMpiVec<doubletype>*)Xin; 328 329 /*y=x+a*y where this->vector is y*/ 330 for(i=0;i<this->m;i++)this->vector[i]=X->vector[i]+a*this->vector[i]; 331 313 332 } 314 333 /*}}}*/ 315 334 /*FUNCTION ToMPISerial{{{*/ 316 335 doubletype* ToMPISerial(void){ 336 337 /*communicator info: */ 338 MPI_Comm comm; 339 int num_procs; 340 341 /*MPI_Allgatherv info: */ 342 int lower_row,upper_row; 343 int* recvcounts=NULL; 344 int* displs=NULL; 317 345 318 printf("IssmMpiVec ToMPISerial not implemented yet!");319 exit(1);346 /*output: */ 347 doubletype* buffer=NULL; 320 348 349 /*initialize comm info: */ 350 comm=IssmComm::GetComm(); 351 num_procs=IssmComm::GetSize(); 352 353 /*Allocate: */ 354 buffer=xNew<doubletype>(M); 355 recvcounts=xNew<int>(num_procs); 356 displs=xNew<int>(num_procs); 357 358 /*recvcounts:*/ 359 MPI_Allgather(&this->m,1,MPI_INT,recvcounts,1,MPI_INT,comm); 360 361 /*get lower_row: */ 362 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,this->m,comm); 363 364 /*displs: */ 365 MPI_Allgather(&lower_row,1,MPI_INT,displs,1,MPI_INT,comm); 366 367 /*All gather:*/ 368 MPI_Allgatherv(this->vector, this->m, MPI_DOUBLE, buffer, recvcounts, displs, MPI_DOUBLE,comm); 369 370 /*free ressources: */ 371 xDelete<int>(recvcounts); 372 xDelete<int>(displs); 373 374 /*return: */ 375 return buffer; 376 321 377 } 322 378 /*}}}*/ 323 379 /*FUNCTION Copy{{{*/ … … 337 393 /*}}}*/ 338 394 /*FUNCTION Norm{{{*/ 339 395 doubletype Norm(NormMode mode){ 396 397 doubletype local_norm; 398 doubletype norm; 399 int i; 340 400 341 printf("IssmMpiVec Norm not implemented yet!"); 342 exit(1); 401 switch(mode){ 402 case NORM_INF: 403 //local_norm=0; for(i=0;i<this->m;i++)local_norm=max(local_norm,fabs(this->vector[i])); 404 local_norm=0; for(i=0;i<this->m;i++)local_norm=max(local_norm,this->vector[i]); 405 MPI_Reduce(&local_norm, &norm, 1, MPI_DOUBLE, MPI_MAX, 0, IssmComm::GetComm()); 406 return norm; 407 break; 408 case NORM_TWO: 409 local_norm=0; 410 for(i=0;i<this->m;i++)local_norm+=pow(this->vector[i],2); 411 MPI_Reduce(&local_norm, &norm, 1, MPI_DOUBLE, MPI_SUM, 0, IssmComm::GetComm()); 412 return sqrt(norm); 413 break; 414 default: 415 _error_("unknown norm !"); 416 break; 417 } 343 418 } 344 419 /*}}}*/ 345 420 /*FUNCTION Scale{{{*/ … … 391 466 } 392 467 /*}}}*/ 393 468 }; 394 #endif //#ifndef _ISSM_MPI_VEC_H_ 469 #endif //#ifndef _ISSM_MPI_VEC_H_ -
../trunk-jpl/src/c/toolkits/issm/IssmSolver.cpp
1 /*!\file SolverxSeq 2 * \brief implementation of sequential solver using the GSL librarie 3 */ 4 5 #ifdef HAVE_CONFIG_H 6 #include <config.h> 7 #else 8 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" 9 #endif 10 #include <cstring> 11 12 #include "./IssmSolver.h" 13 #include "../../shared/shared.h" 14 #include "../../classes/classes.h" 15 #include "../../include/include.h" 16 #include "../../io/io.h" 17 18 #ifdef _HAVE_GSL_ 19 #include <gsl/gsl_linalg.h> 20 #endif 21 22 void IssmSolve(IssmVec<IssmDouble>** pout,IssmMat<IssmDouble>* Kffin, IssmVec<IssmDouble>* pfin, Parameters* parameters){/*{{{*/ 23 24 /*First off, we assume that the type of IssmVec is IssmSeqVec and the type of IssmMat is IssmDenseMat. So downcast: */ 25 IssmSeqVec<IssmDouble>* pf=(IssmSeqVec<IssmDouble>*)pfin->vector; 26 IssmDenseMat<IssmDouble>* Kff=(IssmDenseMat<IssmDouble>*)Kffin->matrix; 27 28 #ifdef _HAVE_GSL_ 29 /*Intermediary: */ 30 int M,N,N2; 31 IssmSeqVec<IssmDouble> *uf = NULL; 32 33 Kff->GetSize(&M,&N); 34 pf->GetSize(&N2); 35 36 if(N!=N2)_error_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !"); 37 if(M!=N)_error_("Stiffness matrix should be square!"); 38 #ifdef _HAVE_ADOLC_ 39 ensureContiguousLocations(N); 40 #endif 41 IssmDouble *x = xNew<IssmDouble>(N); 42 #ifdef _HAVE_ADOLC_ 43 SolverxSeq(x,Kff->matrix,pf->vector,N,parameters); 44 #else 45 SolverxSeq(x,Kff->matrix,pf->vector,N); 46 #endif 47 48 uf=new IssmSeqVec<IssmDouble>(x,N); 49 xDelete(x); 50 51 /*Assign output pointers:*/ 52 IssmVec<IssmDouble>* out=new IssmVec<IssmDouble>(); 53 out->vector=(IssmAbsVec<IssmDouble>*)uf; 54 *pout=out; 55 56 #else 57 _error_("GSL support not compiled in!"); 58 #endif 59 60 }/*}}}*/ 61 void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/ 62 63 /*Allocate output*/ 64 double* X = xNew<double>(n); 65 66 /*Solve*/ 67 SolverxSeq(X,A,B,n); 68 69 /*Assign output pointer*/ 70 *pX=X; 71 } 72 /*}}}*/ 73 void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/ 74 #ifdef _HAVE_GSL_ 75 /*GSL Matrices and vectors: */ 76 int s; 77 gsl_matrix_view a; 78 gsl_vector_view b,x; 79 gsl_permutation *p = NULL; 80 // for (int i=0; i<n*n; ++i) std::cout << "SolverxSeq A["<< i << "]=" << A[i] << std::endl; 81 // for (int i=0; i<n; ++i) std::cout << "SolverxSeq b["<< i << "]=" << B[i] << std::endl; 82 /*A will be modified by LU decomposition. Use copy*/ 83 double* Acopy = xNew<double>(n*n); 84 xMemCpy(Acopy,A,n*n); 85 86 /*Initialize gsl matrices and vectors: */ 87 a = gsl_matrix_view_array (Acopy,n,n); 88 b = gsl_vector_view_array (B,n); 89 x = gsl_vector_view_array (X,n); 90 91 /*Run LU and solve: */ 92 p = gsl_permutation_alloc (n); 93 gsl_linalg_LU_decomp (&a.matrix, p, &s); 94 gsl_linalg_LU_solve (&a.matrix, p, &b.vector, &x.vector); 95 96 /*Clean up and assign output pointer*/ 97 xDelete(Acopy); 98 gsl_permutation_free(p); 99 #endif 100 } 101 /*}}}*/ 102 103 #ifdef _HAVE_ADOLC_ 104 int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y){ /*{{{*/ 105 SolverxSeq(y,x, x+m*m, m); // x is where the matrix starts, x+m*m is where the right-hand side starts 106 return 0; 107 } /*}}}*/ 108 int EDF_fos_forward_for_solverx(int n, IssmPDouble *inVal, IssmPDouble *inDeriv, int m, IssmPDouble *outVal, IssmPDouble *outDeriv) { /*{{{*/ 109 #ifdef _HAVE_GSL_ 110 // for (int i=0; i<m*m; ++i) std::cout << "EDF_fos_forward_for_solverx A["<< i << "]=" << inVal[i] << std::endl; 111 // for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx b["<< i << "]=" << inVal[i+m*m] << std::endl; 112 // the matrix will be modified by LU decomposition. Use gsl_A copy 113 double* Acopy = xNew<double>(m*m); 114 xMemCpy(Acopy,inVal,m*m); 115 /*Initialize gsl matrices and vectors: */ 116 gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m); 117 gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m 118 gsl_permutation *perm_p = gsl_permutation_alloc (m); 119 int signPerm; 120 // factorize 121 gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm); 122 gsl_vector *gsl_x_p = gsl_vector_alloc (m); 123 // solve for the value 124 gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p); 125 /*Copy result*/ 126 xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m); 127 gsl_vector_free(gsl_x_p); 128 // for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx x["<< i << "]=" << outVal[i] << std::endl; 129 // solve for the derivatives acc. to A * dx = r with r=db - dA * x 130 // compute the RHS 131 double* r=xNew<double>(m); 132 for (int i=0; i<m; i++) { 133 r[i]=inDeriv[m*m+i]; // this is db[i] 134 for (int j=0;j<m; j++) { 135 r[i]-=inDeriv[i*m+j]*outVal[j]; // this is dA[i][j]*x[j] 136 } 137 } 138 gsl_vector_view gsl_r=gsl_vector_view_array(r,m); 139 gsl_vector *gsl_dx_p = gsl_vector_alloc(m); 140 gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p); 141 xMemCpy(outDeriv,gsl_vector_ptr(gsl_dx_p,0),m); 142 gsl_vector_free(gsl_dx_p); 143 xDelete(r); 144 gsl_permutation_free(perm_p); 145 xDelete(Acopy); 146 #endif 147 return 0; 148 } /*}}}*/ 149 int EDF_fov_forward_for_solverx(int n, IssmPDouble *inVal, int directionCount, IssmPDouble **inDeriv, int m, IssmPDouble *outVal, IssmPDouble **outDeriv) { /*{{{*/ 150 #ifdef _HAVE_GSL_ 151 // the matrix will be modified by LU decomposition. Use gsl_A copy 152 double* Acopy = xNew<double>(m*m); 153 xMemCpy(Acopy,inVal,m*m); 154 /*Initialize gsl matrices and vectors: */ 155 gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m); 156 gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m 157 gsl_permutation *perm_p = gsl_permutation_alloc (m); 158 int signPerm; 159 // factorize 160 gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm); 161 gsl_vector *gsl_x_p = gsl_vector_alloc (m); 162 // solve for the value 163 gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p); 164 /*Copy result*/ 165 xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m); 166 gsl_vector_free(gsl_x_p); 167 // solve for the derivatives acc. to A * dx = r with r=db - dA * x 168 double* r=xNew<double>(m); 169 gsl_vector *gsl_dx_p = gsl_vector_alloc(m); 170 for (int dir=0;dir<directionCount;++dir) { 171 // compute the RHS 172 for (int i=0; i<m; i++) { 173 r[i]=inDeriv[m*m+i][dir]; // this is db[i] 174 for (int j=0;j<m; j++) { 175 r[i]-=inDeriv[i*m+j][dir]*outVal[j]; // this is dA[i][j]*x[j] 176 } 177 } 178 gsl_vector_view gsl_r=gsl_vector_view_array(r,m); 179 gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p); 180 // reuse r 181 xMemCpy(r,gsl_vector_ptr(gsl_dx_p,0),m); 182 for (int i=0; i<m; i++) { 183 outDeriv[i][dir]=r[i]; 184 } 185 } 186 gsl_vector_free(gsl_dx_p); 187 xDelete(r); 188 gsl_permutation_free(perm_p); 189 xDelete(Acopy); 190 #endif 191 return 0; 192 } 193 /*}}}*/ 194 int EDF_fos_reverse_for_solverx(int m, double *dp_U, int n, double *dp_Z, double* dp_x, double* dp_y) { /*{{{*/ 195 // copy to transpose the matrix 196 double* transposed=xNew<double>(m*m); 197 for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j]; 198 gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m); 199 // the adjoint of the solution is our right-hand side 200 gsl_vector_view x_bar=gsl_vector_view_array(dp_U,m); 201 // the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute: 202 gsl_vector_view b_bar=gsl_vector_view_array(dp_Z+m*m,m); 203 gsl_permutation *perm_p = gsl_permutation_alloc (m); 204 int permSign; 205 gsl_linalg_LU_decomp (&aTransposed.matrix, perm_p, &permSign); 206 gsl_linalg_LU_solve (&aTransposed.matrix, perm_p, &x_bar.vector, &b_bar.vector); 207 // now do the adjoint of the matrix 208 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]; 209 gsl_permutation_free(perm_p); 210 xDelete(transposed); 211 return 0; 212 } 213 /*}}}*/ 214 int EDF_fov_reverse_for_solverx(int m, int p, double **dpp_U, int n, double **dpp_Z, double* dp_x, double* dp_y) { /*{{{*/ 215 // copy to transpose the matrix 216 double* transposed=xNew<double>(m*m); 217 for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j]; 218 gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m); 219 gsl_permutation *perm_p = gsl_permutation_alloc (m); 220 int permSign; 221 gsl_linalg_LU_decomp (&aTransposed.matrix, perm_p, &permSign); 222 for (int weightsRowIndex=0;weightsRowIndex<p;++weightsRowIndex) { 223 // the adjoint of the solution is our right-hand side 224 gsl_vector_view x_bar=gsl_vector_view_array(dpp_U[weightsRowIndex],m); 225 // the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute: 226 gsl_vector_view b_bar=gsl_vector_view_array(dpp_Z[weightsRowIndex]+m*m,m); 227 gsl_linalg_LU_solve (&aTransposed.matrix, perm_p, &x_bar.vector, &b_bar.vector); 228 // now do the adjoint of the matrix 229 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]; 230 } 231 gsl_permutation_free(perm_p); 232 xDelete(transposed); 233 return 0; 234 } 235 /*}}}*/ 236 void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/ 237 // pack inputs to conform to the EDF-prescribed interface 238 // ensure a contiguous block of locations: 239 ensureContiguousLocations(n*(n+1)); 240 IssmDouble* adoubleEDFin=xNew<IssmDouble>(n*(n+1)); // packed inputs, i.e. matrix and right hand side 241 for(int i=0; i<n*n;i++)adoubleEDFin[i] =A[i]; // pack matrix 242 for(int i=0; i<n; i++)adoubleEDFin[i+n*n]=B[i]; // pack the right hand side 243 IssmPDouble* pdoubleEDFin=xNew<IssmPDouble>(n*(n+1)); // provide space to transfer inputs during call_ext_fct 244 IssmPDouble* pdoubleEDFout=xNew<IssmPDouble>(n); // provide space to transfer outputs during call_ext_fct 245 // call the wrapped solver through the registry entry we retrieve from parameters 246 call_ext_fct(dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p, 247 n*(n+1), pdoubleEDFin, adoubleEDFin, 248 n, pdoubleEDFout,X); 249 // for(int i=0; i<n; i++) {ADOLC_DUMP_MACRO(X[i]);} 250 xDelete(adoubleEDFin); 251 xDelete(pdoubleEDFin); 252 xDelete(pdoubleEDFout); 253 } 254 /*}}}*/ 255 #endif -
../trunk-jpl/src/c/toolkits/issm/issmtoolkit.h
17 17 #include "./IssmMat.h" 18 18 #include "./IssmSeqVec.h" 19 19 #include "./IssmVec.h" 20 #include "./IssmSolver.h" 20 21 21 22 #ifdef _HAVE_MPI_ 22 23 #include "./IssmMpiDenseMat.h" -
../trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h
271 271 /*}}}*/ 272 272 /*FUNCTION Norm{{{*/ 273 273 doubletype Norm(NormMode mode){ 274 _error_("not supported yet!"); 274 275 276 doubletype norm,local_norm; 277 doubletype absolute; 278 int i,j; 279 280 switch(mode){ 281 case NORM_INF: 282 local_norm=0; 283 for(i=0;i<this->M;i++){ 284 absolute=0; 285 for(j=0;j<this->N;j++){ 286 absolute+=fabs(this->matrix[N*i+j]); 287 } 288 local_norm=max(local_norm,absolute); 289 } 290 MPI_Reduce(&local_norm, &norm, 1, MPI_DOUBLE, MPI_MAX, 0, IssmComm::GetComm()); 291 MPI_Bcast(&norm,1,MPI_DOUBLE,0,IssmComm::GetComm()); 292 return norm; 293 break; 294 default: 295 _error_("unknown norm !"); 296 break; 297 } 275 298 } 276 299 /*}}}*/ 277 300 /*FUNCTION GetSize{{{*/ 278 301 void GetSize(int* pM,int* pN){ 279 _error_("not supported yet!"); 302 *pM=M; 303 *pN=N; 280 304 } 281 305 /*}}}*/ 282 306 /*FUNCTION GetLocalSize{{{*/ 283 307 void GetLocalSize(int* pM,int* pN){ 284 _error_("not supported yet!"); 308 *pM=m; 309 *pN=N; 285 310 } 286 311 /*}}}*/ 287 312 /*FUNCTION MatMult{{{*/ 288 313 void MatMult(IssmAbsVec<doubletype>* Xin,IssmAbsVec<doubletype>* AXin){ 289 _error_("not supported yet!"); 314 315 316 int i,j; 317 doubletype *X_serial = NULL; 318 319 320 /*A check on the types: */ 321 if(IssmVecTypeFromToolkitOptions()!=MpiEnum)_error_("MatMult operation only possible with 'mpi' vectors"); 322 323 /*Now that we are sure, cast vectors: */ 324 IssmMpiVec<doubletype>* X=(IssmMpiVec<doubletype>*)Xin; 325 IssmMpiVec<doubletype>* AX=(IssmMpiVec<doubletype>*)AXin; 326 327 /*Serialize input Xin: */ 328 X_serial=X->ToMPISerial(); 329 330 /*Every cpu has a serial version of the input vector. Use it to do the Matrix-Vector multiply 331 *locally and plug it into AXin: */ 332 for(i=0;i<this->m;i++){ 333 for(j=0;j<this->N;j++){ 334 AX->vector[i]+=this->matrix[i*N+j]*X_serial[j]; 335 } 336 } 337 338 /*Free ressources: */ 339 xDelete<doubletype>(X_serial); 290 340 } 291 341 /*}}}*/ 292 342 /*FUNCTION Duplicate{{{*/ 293 343 IssmMpiDenseMat<doubletype>* Duplicate(void){ 294 _error_("not supported yet!"); 344 345 IssmMpiDenseMat<doubletype>* dup=new IssmMpiDenseMat<doubletype>(this->matrix,this->M,this->N,0); 346 return dup; 347 295 348 } 296 349 /*}}}*/ 297 350 /*FUNCTION ToSerial{{{*/ … … 313 366 /*}}}*/ 314 367 /*FUNCTION Convert{{{*/ 315 368 void Convert(MatrixType type){ 316 /*do nothing: */369 _error_("not supported yet!"); 317 370 } 318 371 /*}}}*/ 319 372 }; -
../trunk-jpl/src/c/toolkits/issm/IssmSolver.h
1 /*!\file: IssmSolver.h 2 * \brief main hook up from Solver toolkit object (in src/c/classes/toolkits) to the ISSM toolkit 3 */ 4 5 #ifndef _ISSM_SOLVER_H_ 6 #define _ISSM_SOLVER_H_ 7 8 /*Headers:*/ 9 /*{{{*/ 10 #ifdef HAVE_CONFIG_H 11 #include <config.h> 12 #else 13 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" 14 #endif 15 16 #include "../../include/types.h" 17 18 /*}}}*/ 19 20 template <class doubletype> class IssmVec; 21 template <class doubletype> class IssmMat; 22 class Parameters; 23 24 void IssmSolve(IssmVec<IssmDouble>** puf,IssmMat<IssmDouble>* Kff, IssmVec<IssmDouble>* pf,Parameters* parameters); 25 void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n); 26 void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n); 27 28 #if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_) 29 void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters); 30 // call back functions: 31 ADOLC_ext_fct EDF_for_solverx; 32 ADOLC_ext_fct_fos_forward EDF_fos_forward_for_solverx; 33 ADOLC_ext_fct_fos_reverse EDF_fos_reverse_for_solverx; 34 ADOLC_ext_fct_fov_forward EDF_fov_forward_for_solverx; 35 ADOLC_ext_fct_fov_reverse EDF_fov_reverse_for_solverx; 36 #endif 37 38 #endif -
../trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.h
1 /*!\file: PetscMat.h 2 */ 3 4 #ifndef _PETSC_SOLVER_H_ 5 #define _PETSC_SOLVER_H_ 6 7 /*Headers:*/ 8 /*{{{*/ 9 #ifdef HAVE_CONFIG_H 10 #include <config.h> 11 #else 12 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" 13 #endif 14 15 #include "../petscincludes.h" 16 class Parameters; 17 18 /*}}}*/ 19 20 void PetscSolve(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters); 21 void SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters); 22 void DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum); 23 24 #endif //#ifndef _PETSCSOLVER_H_ -
../trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.cpp
1 /*!\file SolverxPetsc 2 * \brief Petsc implementation of solver 3 */ 4 5 6 #ifdef HAVE_CONFIG_H 7 #include <config.h> 8 #else 9 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" 10 #endif 11 12 #include "./PetscSolver.h" 13 #include "../../../shared/Numerics/Verbosity.h" 14 #include "../../../shared/Alloc/alloc.h" 15 #include "../../../shared/Exceptions/exceptions.h" 16 #include "../../../classes/IssmComm.h" 17 #include "../../../EnumDefinitions/EnumDefinitions.h" 18 19 void PetscSolve(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters){ /*{{{*/ 20 21 PetscVec* uf=new PetscVec(); 22 23 Vec uf0_vector = NULL; 24 Vec df_vector = NULL; 25 26 if(uf0) uf0_vector = uf0->vector; 27 if(df) df_vector = df->vector; 28 29 SolverxPetsc(&uf->vector, Kff->matrix, pf->vector, uf0_vector, df_vector, parameters); 30 31 *puf=uf; 32 33 } 34 /*}}}*/ 35 void SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){ /*{{{*/ 36 37 /*Output: */ 38 Vec uf = NULL; 39 40 /*Intermediary: */ 41 int local_m,local_n,global_m,global_n; 42 43 /*Solver */ 44 KSP ksp = NULL; 45 PC pc = NULL; 46 int iteration_number; 47 int solver_type; 48 bool fromlocalsize = true; 49 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2) 50 PetscTruth flag,flg; 51 #else 52 PetscBool flag,flg; 53 #endif 54 55 /*Stokes: */ 56 IS isv=NULL; 57 IS isp=NULL; 58 59 #if _PETSC_MAJOR_ >= 3 60 char ksp_type[50]; 61 #endif 62 63 /*Display message*/ 64 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2) 65 if(VerboseSolver())PetscOptionsPrint(stdout); 66 #else 67 if(VerboseSolver())PetscOptionsView(PETSC_VIEWER_STDOUT_WORLD); 68 #endif 69 70 /*First, check that f-set is not NULL, i.e. model is fully constrained:*/ 71 _assert_(Kff); 72 MatGetSize(Kff,&global_m,&global_n); _assert_(global_m==global_n); 73 if(!global_n){ 74 *puf=NewVec(0,IssmComm::GetComm()); return; 75 } 76 77 /*Initial guess */ 78 /*Now, check that we are not giving an initial guess to the solver, if we are running a direct solver: */ 79 #if _PETSC_MAJOR_ >= 3 80 PetscOptionsGetString(PETSC_NULL,"-ksp_type",ksp_type,49,&flg); 81 if (strcmp(ksp_type,"preonly")==0)uf0=NULL; 82 #endif 83 84 /*If initial guess for the solution exists, use it to create uf, otherwise, 85 * duplicate the right hand side so that the solution vector has the same structure*/ 86 if(uf0){ 87 VecDuplicate(uf0,&uf); VecCopy(uf0,uf); 88 } 89 else{ 90 MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,IssmComm::GetComm(),fromlocalsize); 91 } 92 93 /*Process petsc options to see if we are using special types of external solvers*/ 94 PetscOptionsDetermineSolverType(&solver_type); 95 96 /*Check the solver is available*/ 97 if(solver_type==MUMPSPACKAGE_LU || solver_type==MUMPSPACKAGE_CHOL){ 98 #if _PETSC_MAJOR_ >=3 99 #ifndef _HAVE_MUMPS_ 100 _error_("requested MUMPS solver, which was not compiled into ISSM!\n"); 101 #endif 102 #endif 103 } 104 105 /*Prepare solver*/ 106 KSPCreate(IssmComm::GetComm(),&ksp); 107 KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN); 108 KSPSetFromOptions(ksp); 109 110 #if _PETSC_MAJOR_==3 111 /*Specific solver?: */ 112 KSPGetPC(ksp,&pc); 113 if (solver_type==MUMPSPACKAGE_LU){ 114 #if _PETSC_MINOR_==1 115 PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS); 116 #else 117 PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS); 118 #endif 119 } 120 121 /*Stokes: */ 122 if (solver_type==StokesSolverEnum){ 123 /*Make indices out of doftypes: */ 124 if(!df)_error_("need doftypes for Stokes solver!\n"); 125 DofTypesToIndexSet(&isv,&isp,df,StokesSolverEnum); 126 127 /*Set field splits: */ 128 KSPGetPC(ksp,&pc); 129 #if _PETSC_MINOR_==1 130 PCFieldSplitSetIS(pc,isv); 131 PCFieldSplitSetIS(pc,isp); 132 #else 133 PCFieldSplitSetIS(pc,PETSC_NULL,isv); 134 PCFieldSplitSetIS(pc,PETSC_NULL,isp); 135 #endif 136 137 } 138 #endif 139 140 /*If there is an initial guess for the solution, use it 141 * except if we are using the MUMPS direct solver 142 * where any initial guess will crash Petsc*/ 143 if (uf0){ 144 if((solver_type!=MUMPSPACKAGE_LU) && (solver_type!=MUMPSPACKAGE_CHOL) && (solver_type!=SPOOLESPACKAGE_LU) && (solver_type!=SPOOLESPACKAGE_CHOL) && (solver_type!=SUPERLUDISTPACKAGE)){ 145 KSPSetInitialGuessNonzero(ksp,PETSC_TRUE); 146 } 147 } 148 149 /*Solve: */ 150 if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD); 151 KSPSolve(ksp,pf,uf); 152 153 /*Check convergence*/ 154 KSPGetIterationNumber(ksp,&iteration_number); 155 if (iteration_number<0) _error_("Solver diverged at iteration number: " << -iteration_number); 156 157 /*Free resources:*/ 158 KSPFree(&ksp); 159 160 /*Assign output pointers:*/ 161 *puf=uf; 162 } 163 /*}}}*/ 164 void DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum){ /*{{{*/ 165 166 /*output: */ 167 IS isv=NULL; 168 IS isp=NULL; 169 170 int start,end; 171 IssmDouble* df_local=NULL; 172 int df_local_size; 173 int i; 174 175 int* pressure_indices=NULL; 176 int* velocity_indices=NULL; 177 int pressure_num=0; 178 int velocity_num=0; 179 int pressure_count=0; 180 int velocity_count=0; 181 182 if(typeenum==StokesSolverEnum){ 183 184 /*Ok, recover doftypes vector values and indices: */ 185 VecGetOwnershipRange(df,&start,&end); 186 VecGetLocalSize(df,&df_local_size); 187 VecGetArray(df,&df_local); 188 189 pressure_num=0; 190 velocity_num=0; 191 for(i=0;i<df_local_size;i++){ 192 if (df_local[i]==PressureEnum)pressure_num++; 193 else velocity_num++; 194 } 195 196 /*Allocate indices: */ 197 if(pressure_num)pressure_indices=xNew<int>(pressure_num); 198 if(velocity_num)velocity_indices=xNew<int>(velocity_num); 199 200 pressure_count=0; 201 velocity_count=0; 202 for(i=0;i<df_local_size;i++){ 203 if (df_local[i]==PressureEnum){ 204 pressure_indices[pressure_count]=start+i; 205 pressure_count++; 206 } 207 if (df_local[i]==VelocityEnum){ 208 velocity_indices[velocity_count]=start+i; 209 velocity_count++; 210 } 211 } 212 VecRestoreArray(df,&df_local); 213 214 /*Create indices sets: */ 215 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2) 216 ISCreateGeneral(IssmComm::GetComm(),pressure_num,pressure_indices,&isp); 217 ISCreateGeneral(IssmComm::GetComm(),velocity_num,velocity_indices,&isv); 218 #else 219 ISCreateGeneral(IssmComm::GetComm(),pressure_num,pressure_indices,PETSC_COPY_VALUES,&isp); 220 ISCreateGeneral(IssmComm::GetComm(),velocity_num,velocity_indices,PETSC_COPY_VALUES,&isv); 221 #endif 222 } 223 224 /*Free ressources:*/ 225 xDelete<int>(pressure_indices); 226 xDelete<int>(velocity_indices); 227 228 /*Assign output pointers:*/ 229 *pisv=isv; 230 *pisp=isp; 231 } 232 /*}}}*/ -
../trunk-jpl/src/c/toolkits/petsc/objects/petscobjects.h
7 7 8 8 #include "./PetscMat.h" 9 9 #include "./PetscVec.h" 10 #include "./PetscSolver.h" 10 11 11 12 #endif -
../trunk-jpl/src/c/Makefile.am
118 118 ./classes/matrix/ElementMatrix.cpp\ 119 119 ./classes/matrix/ElementVector.h\ 120 120 ./classes/matrix/ElementVector.cpp\ 121 ./classes/matrix/Matrix.h\ 122 ./classes/matrix/Vector.h\ 121 ./classes/toolkits/toolkitobjects.h\ 122 ./classes/toolkits/Matrix.h\ 123 ./classes/toolkits/Vector.h\ 124 ./classes/toolkits/Solver.h\ 123 125 ./classes/objects/Params/Param.h\ 124 126 ./classes/objects/Params/GenericParam.h\ 125 127 ./classes/objects/Params/BoolParam.cpp\ … … 228 230 ./toolkits/issm/IssmMat.h\ 229 231 ./toolkits/issm/IssmSeqVec.h\ 230 232 ./toolkits/issm/IssmVec.h\ 233 ./toolkits/issm/IssmSolver.h\ 234 ./toolkits/issm/IssmSolver.cpp\ 231 235 ./toolkits/triangle/triangleincludes.h\ 232 236 ./toolkitsenums.h\ 233 237 ./toolkits.h\ … … 319 323 ./modules/ResetCoordinateSystemx/ResetCoordinateSystemx.cpp\ 320 324 ./modules/Solverx/Solverx.cpp\ 321 325 ./modules/Solverx/Solverx.h\ 322 ./modules/Solverx/SolverxSeq.cpp\323 326 ./modules/VecMergex/VecMergex.cpp\ 324 327 ./modules/VecMergex/VecMergex.h\ 325 328 ./modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp\ … … 759 762 ./toolkits/petsc/objects/PetscMat.cpp\ 760 763 ./toolkits/petsc/objects/PetscVec.h\ 761 764 ./toolkits/petsc/objects/PetscVec.cpp\ 762 ./toolkits/petsc/ petscincludes.h\763 ./ modules/Solverx/SolverxPetsc.cpp\764 ./ modules/Solverx/DofTypesToIndexSet.cpp765 ./toolkits/petsc/objects/PetscSolver.cpp\ 766 ./toolkits/petsc/objects/PetscSolver.h\ 767 ./toolkits/petsc/petscincludes.h 765 768 766 769 #}}} 767 770 #Mpi sources {{{ -
../trunk-jpl/src/c/classes/classes.h
11 11 /*matrix: */ 12 12 #include "./matrix/matrixobjects.h" 13 13 14 /*toolkit objects: */ 15 #include "./toolkits/toolkitobjects.h" 16 14 17 /*bamg: */ 15 18 #include "./bamg/bamgobjects.h" 16 19 -
../trunk-jpl/src/c/classes/matrix/Vector.h
1 /*!\file: Vector.h2 * \brief wrapper to vector objects. The goal is to control which API (PETSc,Scalpack, Plapack?)3 * implements our underlying vector format.4 */5 6 #ifndef _VECTOR_H_7 #define _VECTOR_H_8 9 /*Headers:*/10 /*{{{*/11 #ifdef HAVE_CONFIG_H12 #include <config.h>13 #else14 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"15 #endif16 #include <cstring>17 #include "../../toolkits/toolkits.h"18 #include "../../EnumDefinitions/EnumDefinitions.h"19 /*}}}*/20 21 enum vectortype { PetscVecType, IssmVecType };22 23 template <class doubletype>24 class Vector{25 26 public:27 28 int type;29 #ifdef _HAVE_PETSC_30 PetscVec* pvector;31 #endif32 IssmVec<doubletype>* ivector;33 34 /*Vector constructors, destructors */35 Vector(){ /*{{{*/36 37 InitCheckAndSetType();38 }39 /*}}}*/40 Vector(int M,bool fromlocalsize=false){ /*{{{*/41 42 InitCheckAndSetType();43 44 if(type==PetscVecType){45 #ifdef _HAVE_PETSC_46 this->pvector=new PetscVec(M,fromlocalsize);47 #endif48 }49 else this->ivector=new IssmVec<doubletype>(M,fromlocalsize);50 51 }52 /*}}}*/53 Vector(int m,int M){ /*{{{*/54 55 InitCheckAndSetType();56 57 if(type==PetscVecType){58 #ifdef _HAVE_PETSC_59 this->pvector=new PetscVec(m,M);60 #endif61 }62 else this->ivector=new IssmVec<doubletype>(m,M);63 }64 /*}}}*/65 Vector(doubletype* serial_vec,int M){ /*{{{*/66 67 InitCheckAndSetType();68 69 if(type==PetscVecType){70 #ifdef _HAVE_PETSC_71 this->pvector=new PetscVec(serial_vec,M);72 #endif73 }74 else this->ivector=new IssmVec<doubletype>(serial_vec,M);75 }76 /*}}}*/77 ~Vector(){ /*{{{*/78 79 if(type==PetscVecType){80 #ifdef _HAVE_PETSC_81 delete this->pvector;82 #endif83 }84 else delete this->ivector;85 }86 /*}}}*/87 #ifdef _HAVE_PETSC_88 Vector(Vec petsc_vector){ /*{{{*/89 90 this->type=PetscVecType;91 this->ivector=NULL;92 this->pvector=new PetscVec(petsc_vector);93 94 }95 /*}}}*/96 #endif97 void InitCheckAndSetType(void){ /*{{{*/98 99 char* toolkittype=NULL;100 101 #ifdef _HAVE_PETSC_102 pvector=NULL;103 #endif104 ivector=NULL;105 106 /*retrieve toolkittype: */107 toolkittype=ToolkitOptions::GetToolkitType();108 109 /*set vector type: */110 if (strcmp(toolkittype,"petsc")==0){111 #ifdef _HAVE_PETSC_112 type=PetscVecType;113 #else114 _error_("cannot create petsc vector without PETSC compiled!");115 #endif116 }117 else if(strcmp(toolkittype,"issm")==0){118 /*let this choice stand:*/119 type=IssmVecType;120 }121 else _error_("unknow toolkit type ");122 123 /*Free ressources: */124 xDelete<char>(toolkittype);125 }126 /*}}}*/127 128 /*Vector specific routines*/129 /*FUNCTION Echo{{{*/130 void Echo(void){_assert_(this);131 132 if(type==PetscVecType){133 #ifdef _HAVE_PETSC_134 this->pvector->Echo();135 #endif136 }137 else this->ivector->Echo();138 139 }140 /*}}}*/141 /*FUNCTION Assemble{{{*/142 void Assemble(void){_assert_(this);143 144 if(type==PetscVecType){145 #ifdef _HAVE_PETSC_146 this->pvector->Assemble();147 #endif148 }149 else this->ivector->Assemble();150 151 }152 /*}}}*/153 /*FUNCTION SetValues{{{*/154 void SetValues(int ssize, int* list, doubletype* values, InsMode mode){ _assert_(this);155 if(type==PetscVecType){156 #ifdef _HAVE_PETSC_157 this->pvector->SetValues(ssize,list,values,mode);158 #endif159 }160 else this->ivector->SetValues(ssize,list,values,mode);161 162 }163 /*}}}*/164 /*FUNCTION SetValue{{{*/165 void SetValue(int dof, doubletype value, InsMode mode){_assert_(this);166 167 if(type==PetscVecType){168 #ifdef _HAVE_PETSC_169 this->pvector->SetValue(dof,value,mode);170 #endif171 }172 else this->ivector->SetValue(dof,value,mode);173 174 }175 /*}}}*/176 /*FUNCTION GetValue{{{*/177 void GetValue(doubletype* pvalue,int dof){_assert_(this);178 179 if(type==PetscVecType){180 #ifdef _HAVE_PETSC_181 this->pvector->GetValue(pvalue,dof);182 #endif183 }184 else this->ivector->GetValue(pvalue,dof);185 186 }187 /*}}}*/188 /*FUNCTION GetSize{{{*/189 void GetSize(int* pM){_assert_(this);190 191 if(type==PetscVecType){192 #ifdef _HAVE_PETSC_193 this->pvector->GetSize(pM);194 #endif195 }196 else this->ivector->GetSize(pM);197 198 }199 /*}}}*/200 /*FUNCTION IsEmpty{{{*/201 bool IsEmpty(void){202 int M;203 204 _assert_(this);205 this->GetSize(&M);206 207 if(M==0)208 return true;209 else210 return false;211 }212 /*}}}*/213 /*FUNCTION GetLocalSize{{{*/214 void GetLocalSize(int* pM){_assert_(this);215 216 if(type==PetscVecType){217 #ifdef _HAVE_PETSC_218 this->pvector->GetLocalSize(pM);219 #endif220 }221 else this->ivector->GetLocalSize(pM);222 223 }224 /*}}}*/225 /*FUNCTION Duplicate{{{*/226 Vector<doubletype>* Duplicate(void){_assert_(this);227 228 Vector<doubletype>* output=NULL;229 230 output=new Vector<doubletype>();231 232 if(type==PetscVecType){233 #ifdef _HAVE_PETSC_234 output->pvector=this->pvector->Duplicate();235 #endif236 }237 else output->ivector=this->ivector->Duplicate();238 239 return output;240 241 }242 /*}}}*/243 /*FUNCTION Set{{{*/244 void Set(doubletype value){_assert_(this);245 246 if(type==PetscVecType){247 #ifdef _HAVE_PETSC_248 this->pvector->Set(value);249 #endif250 }251 else this->ivector->Set(value);252 253 }254 /*}}}*/255 /*FUNCTION AXPY{{{*/256 void AXPY(Vector* X, doubletype a){_assert_(this);257 258 if(type==PetscVecType){259 #ifdef _HAVE_PETSC_260 this->pvector->AXPY(X->pvector,a);261 #endif262 }263 else this->ivector->AXPY(X->ivector,a);264 265 }266 /*}}}*/267 /*FUNCTION AYPX{{{*/268 void AYPX(Vector* X, doubletype a){_assert_(this);269 270 if(type==PetscVecType){271 #ifdef _HAVE_PETSC_272 this->pvector->AYPX(X->pvector,a);273 #endif274 }275 else this->ivector->AYPX(X->ivector,a);276 }277 /*}}}*/278 /*FUNCTION ToMPISerial{{{*/279 doubletype* ToMPISerial(void){280 281 doubletype* vec_serial=NULL;282 283 _assert_(this);284 if(type==PetscVecType){285 #ifdef _HAVE_PETSC_286 vec_serial=this->pvector->ToMPISerial();287 #endif288 }289 else vec_serial=this->ivector->ToMPISerial();290 291 return vec_serial;292 293 }294 /*}}}*/295 /*FUNCTION Copy{{{*/296 void Copy(Vector* to){_assert_(this);297 298 if(type==PetscVecType){299 #ifdef _HAVE_PETSC_300 this->pvector->Copy(to->pvector);301 #endif302 }303 else this->ivector->Copy(to->ivector);304 }305 /*}}}*/306 /*FUNCTION Norm{{{*/307 doubletype Norm(NormMode norm_type){_assert_(this);308 309 doubletype norm=0;310 311 if(type==PetscVecType){312 #ifdef _HAVE_PETSC_313 norm=this->pvector->Norm(norm_type);314 #endif315 }316 else norm=this->ivector->Norm(norm_type);317 return norm;318 }319 /*}}}*/320 /*FUNCTION Scale{{{*/321 void Scale(doubletype scale_factor){_assert_(this);322 323 if(type==PetscVecType){324 #ifdef _HAVE_PETSC_325 this->pvector->Scale(scale_factor);326 #endif327 }328 else this->ivector->Scale(scale_factor);329 }330 /*}}}*/331 /*FUNCTION Dot{{{*/332 doubletype Dot(Vector* vector){_assert_(this);333 334 doubletype dot;335 336 if(type==PetscVecType){337 #ifdef _HAVE_PETSC_338 dot=this->pvector->Dot(vector->pvector);339 #endif340 }341 else dot=this->ivector->Dot(vector->ivector);342 return dot;343 }344 /*}}}*/345 /*FUNCTION PointwiseDivide{{{*/346 void PointwiseDivide(Vector* x,Vector* y){_assert_(this);347 348 if(type==PetscVecType){349 #ifdef _HAVE_PETSC_350 this->pvector->PointwiseDivide(x->pvector,y->pvector);351 #endif352 }353 else this->ivector->PointwiseDivide(x->ivector,y->ivector);354 }355 /*}}}*/356 };357 #endif //#ifndef _VECTOR_H_ -
../trunk-jpl/src/c/classes/matrix/Matrix.h
1 /*!\file: Matrix.h2 * \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_H12 #include <config.h>13 #else14 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"15 #endif16 #include <cstring>17 #include "../../toolkits/toolkits.h"18 #include "../../EnumDefinitions/EnumDefinitions.h"19 /*}}}*/20 21 enum matrixtype { PetscMatType, IssmMatType };22 23 template <class doubletype> class Vector;24 25 template <class doubletype>26 class Matrix{27 28 public:29 30 int type;31 #ifdef _HAVE_PETSC_32 PetscMat *pmatrix;33 #endif34 IssmMat<doubletype> *imatrix;35 36 /*Matrix constructors, destructors*/37 /*FUNCTION Matrix(){{{*/38 Matrix(){39 InitCheckAndSetType();40 }41 /*}}}*/42 /*FUNCTION Matrix(int M,int N){{{*/43 Matrix(int M,int N){44 45 InitCheckAndSetType();46 47 if(type==PetscMatType){48 #ifdef _HAVE_PETSC_49 this->pmatrix=new PetscMat(M,N);50 #endif51 }52 else{53 this->imatrix=new IssmMat<doubletype>(M,N);54 }55 56 }57 /*}}}*/58 /*FUNCTION Matrix(int m,int n,int M,int N,int* d_nnz,int* o_nnz){{{*/59 Matrix(int m,int n,int M,int N,int* d_nnz,int* o_nnz){60 61 InitCheckAndSetType();62 63 if(type==PetscMatType){64 #ifdef _HAVE_PETSC_65 this->pmatrix=new PetscMat(m,n,M,N,d_nnz,o_nnz);66 #endif67 }68 else{69 this->imatrix=new IssmMat<doubletype>(m,n,M,N,d_nnz,o_nnz);70 }71 72 }73 /*}}}*/74 /*FUNCTION Matrix(int M,int N,IssmDouble sparsity){{{*/75 Matrix(int M,int N,double sparsity){76 77 InitCheckAndSetType();78 79 if(type==PetscMatType){80 #ifdef _HAVE_PETSC_81 this->pmatrix=new PetscMat(M,N,sparsity);82 #endif83 }84 else{85 this->imatrix=new IssmMat<doubletype>(M,N,sparsity);86 }87 }88 /*}}}*/89 /*FUNCTION Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity){{{*/90 Matrix(IssmPDouble* serial_mat,int M,int N,IssmPDouble sparsity){91 92 InitCheckAndSetType();93 94 if(type==PetscMatType){95 #ifdef _HAVE_PETSC_96 this->pmatrix=new PetscMat(serial_mat,M,N,sparsity);97 #endif98 }99 else{100 this->imatrix=new IssmMat<doubletype>(serial_mat,M,N,sparsity);101 }102 103 }104 /*}}}*/105 /*FUNCTION Matrix(int M,int N,int connectivity,int numberofdofspernode){{{*/106 Matrix(int M,int N,int connectivity,int numberofdofspernode){107 108 InitCheckAndSetType();109 110 if(type==PetscMatType){111 #ifdef _HAVE_PETSC_112 this->pmatrix=new PetscMat(M,N,connectivity,numberofdofspernode);113 #endif114 }115 else{116 this->imatrix=new IssmMat<doubletype>(M,N,connectivity,numberofdofspernode);117 }118 119 }120 /*}}}*/121 /*FUNCTION ~Matrix(){{{*/122 ~Matrix(){123 124 if(type==PetscMatType){125 #ifdef _HAVE_PETSC_126 delete this->pmatrix;127 #endif128 }129 else delete this->imatrix;130 131 }132 /*}}}*/133 /*FUNCTION InitCheckAndSetType(){{{*/134 void InitCheckAndSetType(void){135 136 char* toolkittype=NULL;137 138 #ifdef _HAVE_PETSC_139 pmatrix=NULL;140 #endif141 imatrix=NULL;142 143 /*retrieve toolkittype: */144 toolkittype=ToolkitOptions::GetToolkitType();145 146 /*set matrix type: */147 if (strcmp(toolkittype,"petsc")==0){148 #ifdef _HAVE_PETSC_149 type=PetscMatType;150 #else151 _error_("cannot create petsc matrix without PETSC compiled!");152 #endif153 }154 else if(strcmp(toolkittype,"issm")==0){155 /*let this choice stand:*/156 type=IssmMatType;157 }158 else _error_("unknow toolkit type ");159 160 /*Free ressources: */161 xDelete<char>(toolkittype);162 }163 /*}}}*/164 165 /*Matrix specific routines:*/166 /*FUNCTION Echo{{{*/167 void Echo(void){168 _assert_(this);169 170 if(type==PetscMatType){171 #ifdef _HAVE_PETSC_172 this->pmatrix->Echo();173 #endif174 }175 else{176 this->imatrix->Echo();177 }178 179 }180 /*}}}*/181 /*FUNCTION AllocationInfo{{{*/182 void AllocationInfo(void){183 _assert_(this);184 if(type==PetscMatType){185 #ifdef _HAVE_PETSC_186 this->pmatrix->AllocationInfo();187 #endif188 }189 else{190 this->imatrix->AllocationInfo();191 }192 }/*}}}*/193 /*FUNCTION Assemble{{{*/194 void Assemble(void){195 196 if(type==PetscMatType){197 #ifdef _HAVE_PETSC_198 this->pmatrix->Assemble();199 #endif200 }201 else{202 this->imatrix->Assemble();203 }204 }205 /*}}}*/206 /*FUNCTION Norm{{{*/207 IssmDouble Norm(NormMode norm_type){208 209 IssmDouble norm=0;210 211 if(type==PetscMatType){212 #ifdef _HAVE_PETSC_213 norm=this->pmatrix->Norm(norm_type);214 #endif215 }216 else{217 norm=this->imatrix->Norm(norm_type);218 }219 220 return norm;221 }222 /*}}}*/223 /*FUNCTION GetSize{{{*/224 void GetSize(int* pM,int* pN){225 226 if(type==PetscMatType){227 #ifdef _HAVE_PETSC_228 this->pmatrix->GetSize(pM,pN);229 #endif230 }231 else{232 this->imatrix->GetSize(pM,pN);233 }234 235 }236 /*}}}*/237 /*FUNCTION GetLocalSize{{{*/238 void GetLocalSize(int* pM,int* pN){239 240 if(type==PetscMatType){241 #ifdef _HAVE_PETSC_242 this->pmatrix->GetLocalSize(pM,pN);243 #endif244 }245 else{246 this->imatrix->GetLocalSize(pM,pN);247 }248 249 }250 /*}}}*/251 /*FUNCTION MatMult{{{*/252 void MatMult(Vector<doubletype>* X,Vector<doubletype>* AX){253 254 if(type==PetscMatType){255 #ifdef _HAVE_PETSC_256 this->pmatrix->MatMult(X->pvector,AX->pvector);257 #endif258 }259 else{260 this->imatrix->MatMult(X->ivector,AX->ivector);261 }262 263 }264 /*}}}*/265 /*FUNCTION Duplicate{{{*/266 Matrix<doubletype>* Duplicate(void){267 268 Matrix<doubletype>* output=new Matrix<doubletype>();269 270 if(type==PetscMatType){271 #ifdef _HAVE_PETSC_272 output->pmatrix=this->pmatrix->Duplicate();273 #endif274 }275 else{276 output->imatrix=this->imatrix->Duplicate();277 }278 279 return output;280 }281 /*}}}*/282 /*FUNCTION ToSerial{{{*/283 doubletype* ToSerial(void){284 285 doubletype* output=NULL;286 287 if(type==PetscMatType){288 #ifdef _HAVE_PETSC_289 output=this->pmatrix->ToSerial();290 #endif291 }292 else{293 output=this->imatrix->ToSerial();294 }295 296 return output;297 }298 /*}}}*/299 /*FUNCTION SetValues{{{*/300 void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){301 302 if(type==PetscMatType){303 #ifdef _HAVE_PETSC_304 this->pmatrix->SetValues(m,idxm,n,idxn,values,mode);305 #endif306 }307 else{308 this->imatrix->SetValues(m,idxm,n,idxn,values,mode);309 }310 }311 /*}}}*/312 /*FUNCTION Convert{{{*/313 void Convert(MatrixType newtype){314 315 if(type==PetscMatType){316 #ifdef _HAVE_PETSC_317 this->pmatrix->Convert(newtype);318 #endif319 }320 else{321 this->imatrix->Convert(newtype);322 }323 324 }325 /*}}}*/326 };327 328 #endif //#ifndef _MATRIX_H_ -
../trunk-jpl/src/c/classes/matrix/matrixobjects.h
8 8 /*Numerics:*/ 9 9 #include "./ElementMatrix.h" 10 10 #include "./ElementVector.h" 11 #include "./Vector.h"12 #include "./Matrix.h"13 11 14 12 #endif -
../trunk-jpl/src/c/classes/toolkits/Solver.h
1 /*!\file: Solver.h 2 */ 3 4 #ifndef _SOLVER_CLASS_H_ 5 #define _SOLVER_CLASS_H_ 6 7 /*Headers:*/ 8 /*{{{*/ 9 #ifdef HAVE_CONFIG_H 10 #include <config.h> 11 #else 12 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" 13 #endif 14 #include "../../toolkits/toolkits.h" 15 /*}}}*/ 16 17 #include "./Matrix.h" 18 #include "./Vector.h" 19 class Parameters; 20 21 template <class doubletype> 22 class Solver{ 23 24 private: 25 Matrix<doubletype>* Kff; 26 Vector<doubletype>* pf; 27 Vector<doubletype>* uf0; 28 Vector<doubletype>* df; 29 Parameters* parameters; 30 31 public: 32 33 /*Constructors, destructors:*/ 34 /*Solver(){{{*/ 35 Solver(){ 36 } 37 /*}}}*/ 38 /*Solver(Matrix<doubletype>* Kff, Vector<doubletype>* pf, Vector<doubletype>* uf0,Vector<doubletype>* df, Parameters* parameters):{{{*/ 39 Solver(Matrix<doubletype>* Kff_in, Vector<doubletype>* pf_in, Vector<doubletype>* uf0_in,Vector<doubletype>* df_in, Parameters* parameters_in){ 40 41 /*In debugging mode, check that stiffness matrix and load vectors are not NULL (they can be empty)*/ 42 _assert_(Kff_in); 43 _assert_(pf_in); 44 45 /*initialize fields: */ 46 this->Kff=Kff_in; 47 this->pf=pf_in; 48 this->uf0=uf0_in; 49 this->df=df_in; 50 this->parameters=parameters_in; 51 } 52 /*}}}*/ 53 /*~Solver(){{{*/ 54 ~Solver(){ 55 } 56 /*}}}*/ 57 58 /*Methods: */ 59 /*Vector<doubletype Solve(void): {{{*/ 60 Vector<doubletype>* Solve(){ 61 62 /*output: */ 63 Vector<doubletype>* uf=NULL; 64 65 /*Initialize vector: */ 66 uf=new Vector<doubletype>(); 67 68 /*According to matrix type, use specific solvers: */ 69 switch(Kff->type){ 70 #ifdef _HAVE_PETSC_ 71 case PetscMatType:{ 72 PetscVec* uf0_vector = NULL; 73 PetscVec* df_vector = NULL; 74 if(uf0) uf0_vector = uf0->pvector; 75 if(df) df_vector = df->pvector; 76 PetscSolve(&uf->pvector,Kff->pmatrix,pf->pvector,uf0_vector,df_vector,parameters); 77 break; 78 } 79 #endif 80 case IssmMatType:{ 81 IssmSolve(&uf->ivector,Kff->imatrix,pf->ivector,parameters); 82 break; 83 } 84 default: 85 _error_("Matrix type: " << Kff->type << " not supported yet!"); 86 } 87 88 /*allocate output pointer: */ 89 return uf; 90 91 } 92 /*}}}*/ 93 94 }; 95 96 97 98 99 #endif //#ifndef _SOLVER_H_ -
../trunk-jpl/src/c/classes/toolkits/toolkitobjects.h
1 2 /*!\file: toolkitobjects.h 3 * \brief wrappers to toolkits 4 */ 5 6 #ifndef _TOOLKIT_OBJECTS_H_ 7 #define _TOOLKIT_OBJECTS_H_ 8 9 #include "./Vector.h" 10 #include "./Matrix.h" 11 #include "./Solver.h" 12 13 #endif -
../trunk-jpl/src/c/classes/toolkits/Vector.h
1 /*!\file: Vector.h 2 * \brief wrapper to vector objects. The goal is to control which API (PETSc,Scalpack, Plapack?) 3 * implements our underlying vector format. 4 */ 5 6 #ifndef _VECTOR_H_ 7 #define _VECTOR_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 "../../toolkits/toolkits.h" 18 #include "../../EnumDefinitions/EnumDefinitions.h" 19 /*}}}*/ 20 21 enum vectortype { PetscVecType, IssmVecType }; 22 23 template <class doubletype> 24 class Vector{ 25 26 public: 27 28 int type; 29 #ifdef _HAVE_PETSC_ 30 PetscVec* pvector; 31 #endif 32 IssmVec<doubletype>* ivector; 33 34 /*Vector constructors, destructors */ 35 Vector(){ /*{{{*/ 36 37 InitCheckAndSetType(); 38 } 39 /*}}}*/ 40 Vector(int M,bool fromlocalsize=false){ /*{{{*/ 41 42 InitCheckAndSetType(); 43 44 if(type==PetscVecType){ 45 #ifdef _HAVE_PETSC_ 46 this->pvector=new PetscVec(M,fromlocalsize); 47 #endif 48 } 49 else this->ivector=new IssmVec<doubletype>(M,fromlocalsize); 50 51 } 52 /*}}}*/ 53 Vector(int m,int M){ /*{{{*/ 54 55 InitCheckAndSetType(); 56 57 if(type==PetscVecType){ 58 #ifdef _HAVE_PETSC_ 59 this->pvector=new PetscVec(m,M); 60 #endif 61 } 62 else this->ivector=new IssmVec<doubletype>(m,M); 63 } 64 /*}}}*/ 65 Vector(doubletype* serial_vec,int M){ /*{{{*/ 66 67 InitCheckAndSetType(); 68 69 if(type==PetscVecType){ 70 #ifdef _HAVE_PETSC_ 71 this->pvector=new PetscVec(serial_vec,M); 72 #endif 73 } 74 else this->ivector=new IssmVec<doubletype>(serial_vec,M); 75 } 76 /*}}}*/ 77 ~Vector(){ /*{{{*/ 78 79 if(type==PetscVecType){ 80 #ifdef _HAVE_PETSC_ 81 delete this->pvector; 82 #endif 83 } 84 else delete this->ivector; 85 } 86 /*}}}*/ 87 #ifdef _HAVE_PETSC_ 88 Vector(Vec petsc_vector){ /*{{{*/ 89 90 this->type=PetscVecType; 91 this->ivector=NULL; 92 this->pvector=new PetscVec(petsc_vector); 93 94 } 95 /*}}}*/ 96 #endif 97 void InitCheckAndSetType(void){ /*{{{*/ 98 99 char* toolkittype=NULL; 100 101 #ifdef _HAVE_PETSC_ 102 pvector=NULL; 103 #endif 104 ivector=NULL; 105 106 /*retrieve toolkittype: */ 107 toolkittype=ToolkitOptions::GetToolkitType(); 108 109 /*set vector type: */ 110 if (strcmp(toolkittype,"petsc")==0){ 111 #ifdef _HAVE_PETSC_ 112 type=PetscVecType; 113 #else 114 _error_("cannot create petsc vector without PETSC compiled!"); 115 #endif 116 } 117 else if(strcmp(toolkittype,"issm")==0){ 118 /*let this choice stand:*/ 119 type=IssmVecType; 120 } 121 else _error_("unknow toolkit type "); 122 123 /*Free ressources: */ 124 xDelete<char>(toolkittype); 125 } 126 /*}}}*/ 127 128 /*Vector specific routines*/ 129 /*FUNCTION Echo{{{*/ 130 void Echo(void){_assert_(this); 131 132 if(type==PetscVecType){ 133 #ifdef _HAVE_PETSC_ 134 this->pvector->Echo(); 135 #endif 136 } 137 else this->ivector->Echo(); 138 139 } 140 /*}}}*/ 141 /*FUNCTION Assemble{{{*/ 142 void Assemble(void){_assert_(this); 143 144 if(type==PetscVecType){ 145 #ifdef _HAVE_PETSC_ 146 this->pvector->Assemble(); 147 #endif 148 } 149 else this->ivector->Assemble(); 150 151 } 152 /*}}}*/ 153 /*FUNCTION SetValues{{{*/ 154 void SetValues(int ssize, int* list, doubletype* values, InsMode mode){ _assert_(this); 155 if(type==PetscVecType){ 156 #ifdef _HAVE_PETSC_ 157 this->pvector->SetValues(ssize,list,values,mode); 158 #endif 159 } 160 else this->ivector->SetValues(ssize,list,values,mode); 161 162 } 163 /*}}}*/ 164 /*FUNCTION SetValue{{{*/ 165 void SetValue(int dof, doubletype value, InsMode mode){_assert_(this); 166 167 if(type==PetscVecType){ 168 #ifdef _HAVE_PETSC_ 169 this->pvector->SetValue(dof,value,mode); 170 #endif 171 } 172 else this->ivector->SetValue(dof,value,mode); 173 174 } 175 /*}}}*/ 176 /*FUNCTION GetValue{{{*/ 177 void GetValue(doubletype* pvalue,int dof){_assert_(this); 178 179 if(type==PetscVecType){ 180 #ifdef _HAVE_PETSC_ 181 this->pvector->GetValue(pvalue,dof); 182 #endif 183 } 184 else this->ivector->GetValue(pvalue,dof); 185 186 } 187 /*}}}*/ 188 /*FUNCTION GetSize{{{*/ 189 void GetSize(int* pM){_assert_(this); 190 191 if(type==PetscVecType){ 192 #ifdef _HAVE_PETSC_ 193 this->pvector->GetSize(pM); 194 #endif 195 } 196 else this->ivector->GetSize(pM); 197 198 } 199 /*}}}*/ 200 /*FUNCTION IsEmpty{{{*/ 201 bool IsEmpty(void){ 202 int M; 203 204 _assert_(this); 205 this->GetSize(&M); 206 207 if(M==0) 208 return true; 209 else 210 return false; 211 } 212 /*}}}*/ 213 /*FUNCTION GetLocalSize{{{*/ 214 void GetLocalSize(int* pM){_assert_(this); 215 216 if(type==PetscVecType){ 217 #ifdef _HAVE_PETSC_ 218 this->pvector->GetLocalSize(pM); 219 #endif 220 } 221 else this->ivector->GetLocalSize(pM); 222 223 } 224 /*}}}*/ 225 /*FUNCTION Duplicate{{{*/ 226 Vector<doubletype>* Duplicate(void){_assert_(this); 227 228 Vector<doubletype>* output=NULL; 229 230 output=new Vector<doubletype>(); 231 232 if(type==PetscVecType){ 233 #ifdef _HAVE_PETSC_ 234 output->pvector=this->pvector->Duplicate(); 235 #endif 236 } 237 else output->ivector=this->ivector->Duplicate(); 238 239 return output; 240 241 } 242 /*}}}*/ 243 /*FUNCTION Set{{{*/ 244 void Set(doubletype value){_assert_(this); 245 246 if(type==PetscVecType){ 247 #ifdef _HAVE_PETSC_ 248 this->pvector->Set(value); 249 #endif 250 } 251 else this->ivector->Set(value); 252 253 } 254 /*}}}*/ 255 /*FUNCTION AXPY{{{*/ 256 void AXPY(Vector* X, doubletype a){_assert_(this); 257 258 if(type==PetscVecType){ 259 #ifdef _HAVE_PETSC_ 260 this->pvector->AXPY(X->pvector,a); 261 #endif 262 } 263 else this->ivector->AXPY(X->ivector,a); 264 265 } 266 /*}}}*/ 267 /*FUNCTION AYPX{{{*/ 268 void AYPX(Vector* X, doubletype a){_assert_(this); 269 270 if(type==PetscVecType){ 271 #ifdef _HAVE_PETSC_ 272 this->pvector->AYPX(X->pvector,a); 273 #endif 274 } 275 else this->ivector->AYPX(X->ivector,a); 276 } 277 /*}}}*/ 278 /*FUNCTION ToMPISerial{{{*/ 279 doubletype* ToMPISerial(void){ 280 281 doubletype* vec_serial=NULL; 282 283 _assert_(this); 284 if(type==PetscVecType){ 285 #ifdef _HAVE_PETSC_ 286 vec_serial=this->pvector->ToMPISerial(); 287 #endif 288 } 289 else vec_serial=this->ivector->ToMPISerial(); 290 291 return vec_serial; 292 293 } 294 /*}}}*/ 295 /*FUNCTION Copy{{{*/ 296 void Copy(Vector* to){_assert_(this); 297 298 if(type==PetscVecType){ 299 #ifdef _HAVE_PETSC_ 300 this->pvector->Copy(to->pvector); 301 #endif 302 } 303 else this->ivector->Copy(to->ivector); 304 } 305 /*}}}*/ 306 /*FUNCTION Norm{{{*/ 307 doubletype Norm(NormMode norm_type){_assert_(this); 308 309 doubletype norm=0; 310 311 if(type==PetscVecType){ 312 #ifdef _HAVE_PETSC_ 313 norm=this->pvector->Norm(norm_type); 314 #endif 315 } 316 else norm=this->ivector->Norm(norm_type); 317 return norm; 318 } 319 /*}}}*/ 320 /*FUNCTION Scale{{{*/ 321 void Scale(doubletype scale_factor){_assert_(this); 322 323 if(type==PetscVecType){ 324 #ifdef _HAVE_PETSC_ 325 this->pvector->Scale(scale_factor); 326 #endif 327 } 328 else this->ivector->Scale(scale_factor); 329 } 330 /*}}}*/ 331 /*FUNCTION Dot{{{*/ 332 doubletype Dot(Vector* vector){_assert_(this); 333 334 doubletype dot; 335 336 if(type==PetscVecType){ 337 #ifdef _HAVE_PETSC_ 338 dot=this->pvector->Dot(vector->pvector); 339 #endif 340 } 341 else dot=this->ivector->Dot(vector->ivector); 342 return dot; 343 } 344 /*}}}*/ 345 /*FUNCTION PointwiseDivide{{{*/ 346 void PointwiseDivide(Vector* x,Vector* y){_assert_(this); 347 348 if(type==PetscVecType){ 349 #ifdef _HAVE_PETSC_ 350 this->pvector->PointwiseDivide(x->pvector,y->pvector); 351 #endif 352 } 353 else this->ivector->PointwiseDivide(x->ivector,y->ivector); 354 } 355 /*}}}*/ 356 }; 357 #endif //#ifndef _VECTOR_H_ -
../trunk-jpl/src/c/classes/toolkits/Matrix.h
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 "../../toolkits/toolkits.h" 18 #include "../../EnumDefinitions/EnumDefinitions.h" 19 /*}}}*/ 20 21 enum matrixtype { PetscMatType, IssmMatType }; 22 23 template <class doubletype> class Vector; 24 25 template <class doubletype> 26 class Matrix{ 27 28 public: 29 30 int type; 31 #ifdef _HAVE_PETSC_ 32 PetscMat *pmatrix; 33 #endif 34 IssmMat<doubletype> *imatrix; 35 36 /*Matrix constructors, destructors*/ 37 /*FUNCTION Matrix(){{{*/ 38 Matrix(){ 39 InitCheckAndSetType(); 40 } 41 /*}}}*/ 42 /*FUNCTION Matrix(int M,int N){{{*/ 43 Matrix(int M,int N){ 44 45 InitCheckAndSetType(); 46 47 if(type==PetscMatType){ 48 #ifdef _HAVE_PETSC_ 49 this->pmatrix=new PetscMat(M,N); 50 #endif 51 } 52 else{ 53 this->imatrix=new IssmMat<doubletype>(M,N); 54 } 55 56 } 57 /*}}}*/ 58 /*FUNCTION Matrix(int m,int n,int M,int N,int* d_nnz,int* o_nnz){{{*/ 59 Matrix(int m,int n,int M,int N,int* d_nnz,int* o_nnz){ 60 61 InitCheckAndSetType(); 62 63 if(type==PetscMatType){ 64 #ifdef _HAVE_PETSC_ 65 this->pmatrix=new PetscMat(m,n,M,N,d_nnz,o_nnz); 66 #endif 67 } 68 else{ 69 this->imatrix=new IssmMat<doubletype>(m,n,M,N,d_nnz,o_nnz); 70 } 71 72 } 73 /*}}}*/ 74 /*FUNCTION Matrix(int M,int N,IssmDouble sparsity){{{*/ 75 Matrix(int M,int N,double sparsity){ 76 77 InitCheckAndSetType(); 78 79 if(type==PetscMatType){ 80 #ifdef _HAVE_PETSC_ 81 this->pmatrix=new PetscMat(M,N,sparsity); 82 #endif 83 } 84 else{ 85 this->imatrix=new IssmMat<doubletype>(M,N,sparsity); 86 } 87 } 88 /*}}}*/ 89 /*FUNCTION Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity){{{*/ 90 Matrix(IssmPDouble* serial_mat,int M,int N,IssmPDouble sparsity){ 91 92 InitCheckAndSetType(); 93 94 if(type==PetscMatType){ 95 #ifdef _HAVE_PETSC_ 96 this->pmatrix=new PetscMat(serial_mat,M,N,sparsity); 97 #endif 98 } 99 else{ 100 this->imatrix=new IssmMat<doubletype>(serial_mat,M,N,sparsity); 101 } 102 103 } 104 /*}}}*/ 105 /*FUNCTION Matrix(int M,int N,int connectivity,int numberofdofspernode){{{*/ 106 Matrix(int M,int N,int connectivity,int numberofdofspernode){ 107 108 InitCheckAndSetType(); 109 110 if(type==PetscMatType){ 111 #ifdef _HAVE_PETSC_ 112 this->pmatrix=new PetscMat(M,N,connectivity,numberofdofspernode); 113 #endif 114 } 115 else{ 116 this->imatrix=new IssmMat<doubletype>(M,N,connectivity,numberofdofspernode); 117 } 118 119 } 120 /*}}}*/ 121 /*FUNCTION ~Matrix(){{{*/ 122 ~Matrix(){ 123 124 if(type==PetscMatType){ 125 #ifdef _HAVE_PETSC_ 126 delete this->pmatrix; 127 #endif 128 } 129 else delete this->imatrix; 130 131 } 132 /*}}}*/ 133 /*FUNCTION InitCheckAndSetType(){{{*/ 134 void InitCheckAndSetType(void){ 135 136 char* toolkittype=NULL; 137 138 #ifdef _HAVE_PETSC_ 139 pmatrix=NULL; 140 #endif 141 imatrix=NULL; 142 143 /*retrieve toolkittype: */ 144 toolkittype=ToolkitOptions::GetToolkitType(); 145 146 /*set matrix type: */ 147 if (strcmp(toolkittype,"petsc")==0){ 148 #ifdef _HAVE_PETSC_ 149 type=PetscMatType; 150 #else 151 _error_("cannot create petsc matrix without PETSC compiled!"); 152 #endif 153 } 154 else if(strcmp(toolkittype,"issm")==0){ 155 /*let this choice stand:*/ 156 type=IssmMatType; 157 } 158 else _error_("unknow toolkit type "); 159 160 /*Free ressources: */ 161 xDelete<char>(toolkittype); 162 } 163 /*}}}*/ 164 165 /*Matrix specific routines:*/ 166 /*FUNCTION Echo{{{*/ 167 void Echo(void){ 168 _assert_(this); 169 170 if(type==PetscMatType){ 171 #ifdef _HAVE_PETSC_ 172 this->pmatrix->Echo(); 173 #endif 174 } 175 else{ 176 this->imatrix->Echo(); 177 } 178 179 } 180 /*}}}*/ 181 /*FUNCTION AllocationInfo{{{*/ 182 void AllocationInfo(void){ 183 _assert_(this); 184 if(type==PetscMatType){ 185 #ifdef _HAVE_PETSC_ 186 this->pmatrix->AllocationInfo(); 187 #endif 188 } 189 else{ 190 this->imatrix->AllocationInfo(); 191 } 192 }/*}}}*/ 193 /*FUNCTION Assemble{{{*/ 194 void Assemble(void){ 195 196 if(type==PetscMatType){ 197 #ifdef _HAVE_PETSC_ 198 this->pmatrix->Assemble(); 199 #endif 200 } 201 else{ 202 this->imatrix->Assemble(); 203 } 204 } 205 /*}}}*/ 206 /*FUNCTION Norm{{{*/ 207 IssmDouble Norm(NormMode norm_type){ 208 209 IssmDouble norm=0; 210 211 if(type==PetscMatType){ 212 #ifdef _HAVE_PETSC_ 213 norm=this->pmatrix->Norm(norm_type); 214 #endif 215 } 216 else{ 217 norm=this->imatrix->Norm(norm_type); 218 } 219 220 return norm; 221 } 222 /*}}}*/ 223 /*FUNCTION GetSize{{{*/ 224 void GetSize(int* pM,int* pN){ 225 226 if(type==PetscMatType){ 227 #ifdef _HAVE_PETSC_ 228 this->pmatrix->GetSize(pM,pN); 229 #endif 230 } 231 else{ 232 this->imatrix->GetSize(pM,pN); 233 } 234 235 } 236 /*}}}*/ 237 /*FUNCTION GetLocalSize{{{*/ 238 void GetLocalSize(int* pM,int* pN){ 239 240 if(type==PetscMatType){ 241 #ifdef _HAVE_PETSC_ 242 this->pmatrix->GetLocalSize(pM,pN); 243 #endif 244 } 245 else{ 246 this->imatrix->GetLocalSize(pM,pN); 247 } 248 249 } 250 /*}}}*/ 251 /*FUNCTION MatMult{{{*/ 252 void MatMult(Vector<doubletype>* X,Vector<doubletype>* AX){ 253 254 if(type==PetscMatType){ 255 #ifdef _HAVE_PETSC_ 256 this->pmatrix->MatMult(X->pvector,AX->pvector); 257 #endif 258 } 259 else{ 260 this->imatrix->MatMult(X->ivector,AX->ivector); 261 } 262 263 } 264 /*}}}*/ 265 /*FUNCTION Duplicate{{{*/ 266 Matrix<doubletype>* Duplicate(void){ 267 268 Matrix<doubletype>* output=new Matrix<doubletype>(); 269 270 if(type==PetscMatType){ 271 #ifdef _HAVE_PETSC_ 272 output->pmatrix=this->pmatrix->Duplicate(); 273 #endif 274 } 275 else{ 276 output->imatrix=this->imatrix->Duplicate(); 277 } 278 279 return output; 280 } 281 /*}}}*/ 282 /*FUNCTION ToSerial{{{*/ 283 doubletype* ToSerial(void){ 284 285 doubletype* output=NULL; 286 287 if(type==PetscMatType){ 288 #ifdef _HAVE_PETSC_ 289 output=this->pmatrix->ToSerial(); 290 #endif 291 } 292 else{ 293 output=this->imatrix->ToSerial(); 294 } 295 296 return output; 297 } 298 /*}}}*/ 299 /*FUNCTION SetValues{{{*/ 300 void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){ 301 302 if(type==PetscMatType){ 303 #ifdef _HAVE_PETSC_ 304 this->pmatrix->SetValues(m,idxm,n,idxn,values,mode); 305 #endif 306 } 307 else{ 308 this->imatrix->SetValues(m,idxm,n,idxn,values,mode); 309 } 310 } 311 /*}}}*/ 312 /*FUNCTION Convert{{{*/ 313 void Convert(MatrixType newtype){ 314 315 if(type==PetscMatType){ 316 #ifdef _HAVE_PETSC_ 317 this->pmatrix->Convert(newtype); 318 #endif 319 } 320 else{ 321 this->imatrix->Convert(newtype); 322 } 323 324 } 325 /*}}}*/ 326 327 }; 328 329 #endif //#ifndef _MATRIX_H_
Note:
See TracBrowser
for help on using the repository browser.
