Index: ../trunk-jpl/src/c/modules/Solverx/SolverxPetsc.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/Solverx/SolverxPetsc.cpp	(revision 14791)
+++ ../trunk-jpl/src/c/modules/Solverx/SolverxPetsc.cpp	(revision 14792)
@@ -1,159 +0,0 @@
-/*!\file SolverxPetsc
- * \brief Petsc implementation of solver
- */
-
-#include "./Solverx.h"
-#include "../../shared/shared.h"
-#include "../../include/include.h"
-#include "../../io/io.h"
-
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-void	SolverxPetsc(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters){
-
-	PetscVec* uf=new PetscVec();
-
-	Vec uf0_vector = NULL;
-	Vec df_vector  = NULL;
-
-	if(uf0) uf0_vector = uf0->vector;
-	if(df)  df_vector  = df->vector;
-
-	SolverxPetsc(&uf->vector, Kff->matrix, pf->vector, uf0_vector, df_vector, parameters);
-
-	*puf=uf;
-
-}
-void	SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){
-
-	/*Output: */
-	Vec        uf               = NULL;
-
-	/*Intermediary: */
-	int        local_m,local_n,global_m,global_n;
-
-	/*Solver */
-	KSP        ksp              = NULL;
-	PC         pc               = NULL;
-	int        iteration_number;
-	int        solver_type;
-	bool       fromlocalsize    = true;
-	#if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
-	PetscTruth flag,flg;
-	#else
-	PetscBool flag,flg;
-	#endif
-
-	/*Stokes: */
-	IS         isv=NULL;
-	IS         isp=NULL;
-
-	#if _PETSC_MAJOR_ >= 3 
-	char ksp_type[50];
-	#endif
-
-	/*Display message*/
-	if(VerboseModule()) _pprintLine_("   Solving");
-	#if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
-	if(VerboseSolver())PetscOptionsPrint(stdout);
-	#else
-	if(VerboseSolver())PetscOptionsView(PETSC_VIEWER_STDOUT_WORLD);
-	#endif
-
-	/*First, check that f-set is not NULL, i.e. model is fully constrained:*/ 
-	_assert_(Kff);
-	MatGetSize(Kff,&global_m,&global_n); _assert_(global_m==global_n);
-	if(!global_n){
-		*puf=NewVec(0,IssmComm::GetComm()); return;
-	}
-
-	/*Initial guess */
-	/*Now, check that we are not giving an initial guess to the solver, if we are running a direct solver: */
-	#if _PETSC_MAJOR_ >= 3 
-	PetscOptionsGetString(PETSC_NULL,"-ksp_type",ksp_type,49,&flg);
-	if (strcmp(ksp_type,"preonly")==0)uf0=NULL;
-	#endif
-
-	/*If initial guess for the solution exists, use it to create uf, otherwise, 
-	 * duplicate the right hand side so that the solution vector has the same structure*/
-	if(uf0){
-		VecDuplicate(uf0,&uf); VecCopy(uf0,uf);
-	}
-	else{
-		MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,IssmComm::GetComm(),fromlocalsize);
-	}
-
-	/*Process petsc options to see if we are using special types of external solvers*/
-	PetscOptionsDetermineSolverType(&solver_type);
-
-	/*Check the solver is available*/
-	if(solver_type==MUMPSPACKAGE_LU || solver_type==MUMPSPACKAGE_CHOL){
-		#if _PETSC_MAJOR_ >=3
-			#ifndef _HAVE_MUMPS_
-			_error_("requested MUMPS solver, which was not compiled into ISSM!\n");
-			#endif
-		#endif
-	}
-
-	/*Prepare solver*/
-	KSPCreate(IssmComm::GetComm(),&ksp);
-	KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN);
-	KSPSetFromOptions(ksp);
-
-	#if _PETSC_MAJOR_==3
-	/*Specific solver?: */
-	KSPGetPC(ksp,&pc);
-	if (solver_type==MUMPSPACKAGE_LU){
-		#if _PETSC_MINOR_==1
-		PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS);
-		#else
-		PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
-		#endif
-	}
-
-	/*Stokes: */
-	if (solver_type==StokesSolverEnum){
-		/*Make indices out of doftypes: */
-		if(!df)_error_("need doftypes for Stokes solver!\n");
-		DofTypesToIndexSet(&isv,&isp,df,StokesSolverEnum);
-
-		/*Set field splits: */
-		KSPGetPC(ksp,&pc);
-		#if _PETSC_MINOR_==1
-		PCFieldSplitSetIS(pc,isv);
-		PCFieldSplitSetIS(pc,isp);
-		#else
-		PCFieldSplitSetIS(pc,PETSC_NULL,isv);
-		PCFieldSplitSetIS(pc,PETSC_NULL,isp);
-		#endif
-
-	}
-	#endif
-
-	/*If there is an initial guess for the solution, use it
-	 * except if we are using the MUMPS direct solver
-	 * where any initial guess will crash Petsc*/
-	if (uf0){
-		if((solver_type!=MUMPSPACKAGE_LU) && (solver_type!=MUMPSPACKAGE_CHOL) && (solver_type!=SPOOLESPACKAGE_LU) && (solver_type!=SPOOLESPACKAGE_CHOL) && (solver_type!=SUPERLUDISTPACKAGE)){
-			KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
-		}
-	}
-
-	/*Solve: */
-	if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
-	KSPSolve(ksp,pf,uf);
-
-	/*Check convergence*/
-	KSPGetIterationNumber(ksp,&iteration_number);
-	if (iteration_number<0) _error_("Solver diverged at iteration number: " << -iteration_number);
-
-	/*Free resources:*/
-	KSPFree(&ksp);
-
-	/*Assign output pointers:*/
-	*puf=uf;
-}
Index: ../trunk-jpl/src/c/modules/Solverx/DofTypesToIndexSet.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/Solverx/DofTypesToIndexSet.cpp	(revision 14791)
+++ ../trunk-jpl/src/c/modules/Solverx/DofTypesToIndexSet.cpp	(revision 14792)
@@ -1,82 +0,0 @@
-/*!\file Solverx
- * \brief solver
- */
-
-#include "./Solverx.h"
-#include "../../shared/shared.h"
-#include "../../include/include.h"
-
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-void DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum){
-
-	/*output: */
-	IS          isv=NULL;
-	IS          isp=NULL;
-
-	int         start,end;
-	IssmDouble*     df_local=NULL;
-	int         df_local_size;
-	int         i;
-
-	int*     pressure_indices=NULL;
-	int*     velocity_indices=NULL;
-	int      pressure_num=0;
-	int      velocity_num=0;
-	int      pressure_count=0;
-	int      velocity_count=0;
-
-	if(typeenum==StokesSolverEnum){
-
-		/*Ok, recover doftypes vector values and indices: */
-		VecGetOwnershipRange(df,&start,&end);
-		VecGetLocalSize(df,&df_local_size);
-		VecGetArray(df,&df_local);
-
-		pressure_num=0;
-		velocity_num=0;
-		for(i=0;i<df_local_size;i++){
-			if (df_local[i]==PressureEnum)pressure_num++;
-			else velocity_num++;
-		}
-
-		/*Allocate indices: */
-		if(pressure_num)pressure_indices=xNew<int>(pressure_num);
-		if(velocity_num)velocity_indices=xNew<int>(velocity_num);
-
-		pressure_count=0;
-		velocity_count=0;
-		for(i=0;i<df_local_size;i++){
-			if (df_local[i]==PressureEnum){
-				pressure_indices[pressure_count]=start+i;
-				pressure_count++;
-			}
-			if (df_local[i]==VelocityEnum){
-				velocity_indices[velocity_count]=start+i;
-				velocity_count++;
-			}
-		}
-		VecRestoreArray(df,&df_local);
-
-		/*Create indices sets: */
-		#if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
-		ISCreateGeneral(IssmComm::GetComm(),pressure_num,pressure_indices,&isp);
-		ISCreateGeneral(IssmComm::GetComm(),velocity_num,velocity_indices,&isv);
-		#else
-		ISCreateGeneral(IssmComm::GetComm(),pressure_num,pressure_indices,PETSC_COPY_VALUES,&isp);
-		ISCreateGeneral(IssmComm::GetComm(),velocity_num,velocity_indices,PETSC_COPY_VALUES,&isv);
-		#endif
-	}
-
-	/*Free ressources:*/
-	xDelete<int>(pressure_indices);
-	xDelete<int>(velocity_indices);
-
-	/*Assign output pointers:*/
-	*pisv=isv;
-	*pisp=isp;
-}
Index: ../trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp	(revision 14791)
+++ ../trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp	(revision 14792)
@@ -1,255 +0,0 @@
-/*!\file SolverxSeq
- * \brief implementation of sequential solver using the GSL librarie
- */
-
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-#include <cstring>
-
-#include "./Solverx.h"
-#include "../../shared/shared.h"
-#include "../../classes/classes.h"
-#include "../../include/include.h"
-#include "../../io/io.h"
-
-#ifdef _HAVE_GSL_
-#include <gsl/gsl_linalg.h>
-#endif
-
-void SolverxSeq(IssmVec<IssmDouble>** pout,IssmMat<IssmDouble>* Kffin, IssmVec<IssmDouble>* pfin, Parameters* parameters){/*{{{*/
-
-	/*First off, we assume that the type of IssmVec is IssmSeqVec and the type of IssmMat is IssmDenseMat. So downcast: */
-	IssmSeqVec<IssmDouble>* pf=(IssmSeqVec<IssmDouble>*)pfin->vector;
-	IssmDenseMat<IssmDouble>* Kff=(IssmDenseMat<IssmDouble>*)Kffin->matrix;
-
-#ifdef _HAVE_GSL_
-	/*Intermediary: */
-	int M,N,N2;
-	IssmSeqVec<IssmDouble> *uf = NULL;
-
-	Kff->GetSize(&M,&N);
-	pf->GetSize(&N2);
-
-	if(N!=N2)_error_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !");
-	if(M!=N)_error_("Stiffness matrix should be square!");
-#ifdef _HAVE_ADOLC_
-	ensureContiguousLocations(N);
-#endif
-	IssmDouble *x  = xNew<IssmDouble>(N);
-#ifdef _HAVE_ADOLC_
-	SolverxSeq(x,Kff->matrix,pf->vector,N,parameters);
-#else
-	SolverxSeq(x,Kff->matrix,pf->vector,N);
-#endif
-
-	uf=new IssmSeqVec<IssmDouble>(x,N);
-	xDelete(x);
-
-	/*Assign output pointers:*/
-	IssmVec<IssmDouble>* out=new IssmVec<IssmDouble>();
-	out->vector=(IssmAbsVec<IssmDouble>*)uf;
-	*pout=out;
-
-#else
-	_error_("GSL support not compiled in!");
-#endif
-
-}/*}}}*/
-void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
-
-	/*Allocate output*/
-	double* X = xNew<double>(n); 
-
-	/*Solve*/
-	SolverxSeq(X,A,B,n);
-
-	/*Assign output pointer*/
-	*pX=X;
-}
-/*}}}*/
-
-#ifdef _HAVE_ADOLC_
-int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y){ /*{{{*/
-	SolverxSeq(y,x, x+m*m, m); // x is where the matrix starts, x+m*m is where the right-hand side starts
-	return 0;
-} /*}}}*/
-int EDF_fos_forward_for_solverx(int n, IssmPDouble *inVal, IssmPDouble *inDeriv, int m, IssmPDouble *outVal, IssmPDouble *outDeriv) { /*{{{*/
-#ifdef _HAVE_GSL_
-	//  for (int i=0; i<m*m; ++i) std::cout << "EDF_fos_forward_for_solverx A["<< i << "]=" << inVal[i] << std::endl;
-	//  for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx b["<< i << "]=" << inVal[i+m*m] << std::endl;
-	// the matrix will be modified by LU decomposition. Use gsl_A copy
-	double* Acopy = xNew<double>(m*m);
-	xMemCpy(Acopy,inVal,m*m);
-	/*Initialize gsl matrices and vectors: */
-	gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
-	gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
-	gsl_permutation *perm_p = gsl_permutation_alloc (m);
-	int  signPerm;
-	// factorize
-	gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
-	gsl_vector *gsl_x_p = gsl_vector_alloc (m);
-	// solve for the value
-	gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
-	/*Copy result*/
-	xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
-	gsl_vector_free(gsl_x_p);
-	//  for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx x["<< i << "]=" << outVal[i] << std::endl;
-	// solve for the derivatives acc. to A * dx = r  with r=db - dA * x
-	// compute the RHS
-	double* r=xNew<double>(m);
-	for (int i=0; i<m; i++) {
-		r[i]=inDeriv[m*m+i]; // this is db[i]
-		for (int j=0;j<m; j++) {
-			r[i]-=inDeriv[i*m+j]*outVal[j]; // this is dA[i][j]*x[j]
-		}
-	}
-	gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
-	gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
-	gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
-	xMemCpy(outDeriv,gsl_vector_ptr(gsl_dx_p,0),m);
-	gsl_vector_free(gsl_dx_p);
-	xDelete(r);
-	gsl_permutation_free(perm_p);
-	xDelete(Acopy);
-#endif
-	return 0;
-} /*}}}*/
-int EDF_fov_forward_for_solverx(int n, IssmPDouble *inVal, int directionCount, IssmPDouble **inDeriv, int m, IssmPDouble *outVal, IssmPDouble **outDeriv) { /*{{{*/
-#ifdef _HAVE_GSL_
-	// the matrix will be modified by LU decomposition. Use gsl_A copy
-	double* Acopy = xNew<double>(m*m);
-	xMemCpy(Acopy,inVal,m*m);
-	/*Initialize gsl matrices and vectors: */
-	gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
-	gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
-	gsl_permutation *perm_p = gsl_permutation_alloc (m);
-	int  signPerm;
-	// factorize
-	gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
-	gsl_vector *gsl_x_p = gsl_vector_alloc (m);
-	// solve for the value
-	gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
-	/*Copy result*/
-	xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
-	gsl_vector_free(gsl_x_p);
-	// solve for the derivatives acc. to A * dx = r  with r=db - dA * x
-	double* r=xNew<double>(m);
-	gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
-	for (int dir=0;dir<directionCount;++dir) {
-		// compute the RHS
-		for (int i=0; i<m; i++) {
-			r[i]=inDeriv[m*m+i][dir]; // this is db[i]
-			for (int j=0;j<m; j++) {
-				r[i]-=inDeriv[i*m+j][dir]*outVal[j]; // this is dA[i][j]*x[j]
-			}
-		}
-		gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
-		gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
-		// reuse r
-		xMemCpy(r,gsl_vector_ptr(gsl_dx_p,0),m);
-		for (int i=0; i<m; i++) {
-			outDeriv[i][dir]=r[i];
-		}
-	}
-	gsl_vector_free(gsl_dx_p);
-	xDelete(r);
-	gsl_permutation_free(perm_p);
-	xDelete(Acopy);
-#endif
-	return 0;
-}
-/*}}}*/
-int EDF_fos_reverse_for_solverx(int m, double *dp_U, int n, double *dp_Z, double* dp_x, double* dp_y) { /*{{{*/
-	// copy to transpose the matrix
-	double* transposed=xNew<double>(m*m);
-	for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j];
-	gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m);
-	// the adjoint of the solution is our right-hand side
-	gsl_vector_view x_bar=gsl_vector_view_array(dp_U,m);
-	// the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute:
-	gsl_vector_view b_bar=gsl_vector_view_array(dp_Z+m*m,m);
-	gsl_permutation *perm_p = gsl_permutation_alloc (m);
-	int permSign;
-	gsl_linalg_LU_decomp (&aTransposed.matrix, perm_p, &permSign);
-	gsl_linalg_LU_solve (&aTransposed.matrix, perm_p, &x_bar.vector, &b_bar.vector);
-	// now do the adjoint of the matrix
-	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];
-	gsl_permutation_free(perm_p);
-	xDelete(transposed);
-	return 0;
-}
-/*}}}*/
-int EDF_fov_reverse_for_solverx(int m, int p, double **dpp_U, int n, double **dpp_Z, double* dp_x, double* dp_y) { /*{{{*/
-	// copy to transpose the matrix
-	double* transposed=xNew<double>(m*m);
-	for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j];
-	gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m);
-	gsl_permutation *perm_p = gsl_permutation_alloc (m);
-	int permSign;
-	gsl_linalg_LU_decomp (&aTransposed.matrix, perm_p, &permSign);
-	for (int weightsRowIndex=0;weightsRowIndex<p;++weightsRowIndex) {
-		// the adjoint of the solution is our right-hand side
-		gsl_vector_view x_bar=gsl_vector_view_array(dpp_U[weightsRowIndex],m);
-		// the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute:
-		gsl_vector_view b_bar=gsl_vector_view_array(dpp_Z[weightsRowIndex]+m*m,m);
-		gsl_linalg_LU_solve (&aTransposed.matrix, perm_p, &x_bar.vector, &b_bar.vector);
-		// now do the adjoint of the matrix
-		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];
-	}
-	gsl_permutation_free(perm_p);
-	xDelete(transposed);
-	return 0;
-}
-/*}}}*/
-void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
-	// pack inputs to conform to the EDF-prescribed interface
-        // ensure a contiguous block of locations:
-        ensureContiguousLocations(n*(n+1));
-        IssmDouble*  adoubleEDFin=xNew<IssmDouble>(n*(n+1));  // packed inputs, i.e. matrix and right hand side
-        for(int i=0; i<n*n;i++)adoubleEDFin[i]    =A[i];      // pack matrix
-        for(int i=0; i<n;  i++)adoubleEDFin[i+n*n]=B[i];      // pack the right hand side
-        IssmPDouble* pdoubleEDFin=xNew<IssmPDouble>(n*(n+1)); // provide space to transfer inputs during call_ext_fct
-	IssmPDouble* pdoubleEDFout=xNew<IssmPDouble>(n);           // provide space to transfer outputs during call_ext_fct
-	// call the wrapped solver through the registry entry we retrieve from parameters
-	call_ext_fct(dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p,
-	             n*(n+1), pdoubleEDFin, adoubleEDFin,
-	             n, pdoubleEDFout,X);
-	// for(int i=0; i<n;  i++) {ADOLC_DUMP_MACRO(X[i]);}
-	xDelete(adoubleEDFin);
-	xDelete(pdoubleEDFin);
-	xDelete(pdoubleEDFout);
-}
-/*}}}*/
-#endif
-void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
-#ifdef _HAVE_GSL_
-	/*GSL Matrices and vectors: */
-	int              s;
-	gsl_matrix_view  a;
-	gsl_vector_view  b,x;
-	gsl_permutation *p = NULL;
-//	for (int i=0; i<n*n; ++i) std::cout << "SolverxSeq A["<< i << "]=" << A[i] << std::endl;
-//	for (int i=0; i<n; ++i) std::cout << "SolverxSeq b["<< i << "]=" << B[i] << std::endl;
-	/*A will be modified by LU decomposition. Use copy*/
-	double* Acopy = xNew<double>(n*n);
-	xMemCpy(Acopy,A,n*n);
-
-	/*Initialize gsl matrices and vectors: */
-	a = gsl_matrix_view_array (Acopy,n,n);
-	b = gsl_vector_view_array (B,n);
-	x = gsl_vector_view_array (X,n);
-
-	/*Run LU and solve: */
-	p = gsl_permutation_alloc (n);
-	gsl_linalg_LU_decomp (&a.matrix, p, &s);
-	gsl_linalg_LU_solve (&a.matrix, p, &b.vector, &x.vector);
-
-	/*Clean up and assign output pointer*/
-	xDelete(Acopy);
-	gsl_permutation_free(p);
-#endif
-}
-/*}}}*/
Index: ../trunk-jpl/src/c/modules/Solverx/Solverx.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/Solverx/Solverx.cpp	(revision 14791)
+++ ../trunk-jpl/src/c/modules/Solverx/Solverx.cpp	(revision 14792)
@@ -2,48 +2,30 @@
  * \brief solver
  */
 
-#include "./Solverx.h"
-#include "../../shared/shared.h"
-#include "../../include/include.h"
-#include "../../io/io.h"
-
 #ifdef HAVE_CONFIG_H
 	#include <config.h>
 #else
 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
 #endif
 
+#include "./Solverx.h"
+#include "../../shared/shared.h"
+
 void	Solverx(Vector<IssmDouble>** puf, Matrix<IssmDouble>* Kff, Vector<IssmDouble>* pf, Vector<IssmDouble>* uf0,Vector<IssmDouble>* df, Parameters* parameters){
 
+	/*intermediary: */
+	Solver<IssmDouble> *solver=NULL;
+	
 	/*output: */
 	Vector<IssmDouble> *uf=NULL;
 
-	/*In debugging mode, check that stiffness matrix and load vectors are not NULL (they can be empty)*/
-	_assert_(Kff);
-	_assert_(pf);
-	
 	if(VerboseModule()) _pprintLine_("   Solving matrix system");
 
-	/*Initialize vector: */
-	uf=new Vector<IssmDouble>();
+	/*Initialize solver: */
+	solver=new Solver<IssmDouble>(Kff,pf,uf0,df,parameters);
 
-	/*According to matrix type, use specific solvers: */
-	switch(Kff->type){
-		#ifdef _HAVE_PETSC_
-		case PetscMatType:{
-			PetscVec* uf0_vector = NULL;
-			PetscVec* df_vector  = NULL;
-			if(uf0) uf0_vector = uf0->pvector;
-			if(df)  df_vector  = df->pvector;
-			SolverxPetsc(&uf->pvector,Kff->pmatrix,pf->pvector,uf0_vector,df_vector,parameters);
-			break;}
-		#endif
-		case IssmMatType:{
-			SolverxSeq(&uf->ivector,Kff->imatrix,pf->ivector,parameters);
-			break;}
-		default:
-			  _error_("Matrix type: " << Kff->type << " not supported yet!");
-	}
+	/*Solve:*/
+	uf=solver->Solve();
 
 	/*Assign output pointers:*/
 	*puf=uf;
Index: ../trunk-jpl/src/c/modules/Solverx/Solverx.h
===================================================================
--- ../trunk-jpl/src/c/modules/Solverx/Solverx.h	(revision 14791)
+++ ../trunk-jpl/src/c/modules/Solverx/Solverx.h	(revision 14792)
@@ -11,30 +11,9 @@
 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
 #endif
 
-#include "../../classes/objects/objects.h"
-#include "../../toolkits/toolkits.h"
+#include "../../classes/toolkits/toolkitobjects.h"
 
 /* local prototypes: */
 void	Solverx(Vector<IssmDouble>** puf, Matrix<IssmDouble>* Kff, Vector<IssmDouble>* pf, Vector<IssmDouble>* uf0,Vector<IssmDouble>* df, Parameters* parameters);
 
-#ifdef _HAVE_PETSC_
-void	SolverxPetsc(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters);
-void	SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters);
-void  DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum);
-#endif
-
-void SolverxSeq(IssmVec<IssmDouble>** puf,IssmMat<IssmDouble>* Kff, IssmVec<IssmDouble>* pf,Parameters* parameters);
-void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n);
-void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n);
-
-#if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_)
-void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters);
-// call back functions:
-ADOLC_ext_fct EDF_for_solverx;
-ADOLC_ext_fct_fos_forward EDF_fos_forward_for_solverx;
-ADOLC_ext_fct_fos_reverse EDF_fos_reverse_for_solverx;
-ADOLC_ext_fct_fov_forward EDF_fov_forward_for_solverx;
-ADOLC_ext_fct_fov_reverse EDF_fov_reverse_for_solverx;
-#endif
-
 #endif  /* _SOLVERX_H */
Index: ../trunk-jpl/src/c/modules/Solverx/CMakeLists.txt
===================================================================
--- ../trunk-jpl/src/c/modules/Solverx/CMakeLists.txt	(revision 14791)
+++ ../trunk-jpl/src/c/modules/Solverx/CMakeLists.txt	(revision 14792)
@@ -4,6 +4,5 @@
 include_directories(AFTER $ENV{ISSM_DIR}/src/c/modules/Solverx)
 # }}}
 # CORE_SOURCES {{{
-set(CORE_SOURCES $ENV{ISSM_DIR}/src/c/modules/Solverx/Solverx.cpp
-              $ENV{ISSM_DIR}/src/c/modules/Solverx/SolverxSeq.cpp PARENT_SCOPE)
+set(CORE_SOURCES $ENV{ISSM_DIR}/src/c/modules/Solverx/Solverx.cpp PARENT_SCOPE)
 # }}}
Index: ../trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h
===================================================================
--- ../trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h	(revision 14791)
+++ ../trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h	(revision 14792)
@@ -21,6 +21,7 @@
 #include "../../shared/MemOps/xMemCpy.h"
 #include "../../shared/Alloc/alloc.h"
 #include "../../include/macros.h"
+#include "../../io/io.h"
 #ifdef _HAVE_MPI_
 #include "../mpi/mpiincludes.h"
 #endif
@@ -302,22 +303,77 @@
 		/*}}}*/
 		/*FUNCTION AXPY{{{*/
 		void AXPY(IssmAbsVec<doubletype>* Xin, doubletype a){
-			printf("AXPY not implemented yet!");
-			exit(1);
+				
+			int i;
+
+			/*Assume X is of the correct type, and downcast: */
+			IssmMpiVec* X=NULL;
+
+			X=(IssmMpiVec<doubletype>*)Xin;
+
+			/*y=a*x+y where this->vector is y*/
+			for(i=0;i<this->m;i++)this->vector[i]=a*X->vector[i]+this->vector[i];
+
+
 		}
 		/*}}}*/
 		/*FUNCTION AYPX{{{*/
 		void AYPX(IssmAbsVec<doubletype>* Xin, doubletype a){
-			printf("AYPX not implemented yet!");
-			exit(1);
+			int i;
+
+			/*Assume X is of the correct type, and downcast: */
+			IssmMpiVec* X=NULL;
+
+			X=(IssmMpiVec<doubletype>*)Xin;
+
+			/*y=x+a*y where this->vector is y*/
+			for(i=0;i<this->m;i++)this->vector[i]=X->vector[i]+a*this->vector[i];
+
 		}
 		/*}}}*/
 		/*FUNCTION ToMPISerial{{{*/
 		doubletype* ToMPISerial(void){
+
+			/*communicator info: */
+			MPI_Comm comm;
+			int num_procs;
+
+			/*MPI_Allgatherv info: */
+			int  lower_row,upper_row;
+			int* recvcounts=NULL;
+			int* displs=NULL;
 			
-			printf("IssmMpiVec ToMPISerial not implemented yet!");
-			exit(1);
+			/*output: */
+			doubletype* buffer=NULL;
 
+			/*initialize comm info: */
+			comm=IssmComm::GetComm();
+			num_procs=IssmComm::GetSize();
+
+			/*Allocate: */
+			buffer=xNew<doubletype>(M);
+			recvcounts=xNew<int>(num_procs);
+			displs=xNew<int>(num_procs);
+
+			/*recvcounts:*/
+			MPI_Allgather(&this->m,1,MPI_INT,recvcounts,1,MPI_INT,comm);
+			
+			/*get lower_row: */
+			GetOwnershipBoundariesFromRange(&lower_row,&upper_row,this->m,comm);
+
+			/*displs: */
+			MPI_Allgather(&lower_row,1,MPI_INT,displs,1,MPI_INT,comm);
+
+			/*All gather:*/
+			MPI_Allgatherv(this->vector, this->m, MPI_DOUBLE, buffer, recvcounts, displs, MPI_DOUBLE,comm);
+
+			/*free ressources: */
+			xDelete<int>(recvcounts);
+			xDelete<int>(displs);
+
+			/*return: */
+			return buffer;
+
 		}
 		/*}}}*/
 		/*FUNCTION Copy{{{*/
@@ -337,9 +393,28 @@
 		/*}}}*/
 		/*FUNCTION Norm{{{*/
 		doubletype Norm(NormMode mode){
+		
+			doubletype local_norm;
+			doubletype norm;
+			int i;
 
-			printf("IssmMpiVec Norm not implemented yet!");
-			exit(1);
+			switch(mode){
+				case NORM_INF:
+					//local_norm=0; for(i=0;i<this->m;i++)local_norm=max(local_norm,fabs(this->vector[i]));
+					local_norm=0; for(i=0;i<this->m;i++)local_norm=max(local_norm,this->vector[i]);
+					MPI_Reduce(&local_norm, &norm, 1, MPI_DOUBLE, MPI_MAX, 0, IssmComm::GetComm());
+					return norm;
+					break;
+				case NORM_TWO:
+					local_norm=0; 
+					for(i=0;i<this->m;i++)local_norm+=pow(this->vector[i],2);
+					MPI_Reduce(&local_norm, &norm, 1, MPI_DOUBLE, MPI_SUM, 0, IssmComm::GetComm());
+					return sqrt(norm);
+					break;
+				default:
+					_error_("unknown norm !");
+					break;
+			}
 		}
 		/*}}}*/
 		/*FUNCTION Scale{{{*/
@@ -391,4 +466,4 @@
 		}
 		/*}}}*/
 };
-#endif //#ifndef _ISSM_MPI_VEC_H_
+#endif //#ifndef _ISSM_MPI_VEC_H_	
Index: ../trunk-jpl/src/c/toolkits/issm/IssmSolver.cpp
===================================================================
--- ../trunk-jpl/src/c/toolkits/issm/IssmSolver.cpp	(revision 0)
+++ ../trunk-jpl/src/c/toolkits/issm/IssmSolver.cpp	(revision 14792)
@@ -0,0 +1,255 @@
+/*!\file SolverxSeq
+ * \brief implementation of sequential solver using the GSL librarie
+ */
+
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+#include <cstring>
+
+#include "./IssmSolver.h"
+#include "../../shared/shared.h"
+#include "../../classes/classes.h"
+#include "../../include/include.h"
+#include "../../io/io.h"
+
+#ifdef _HAVE_GSL_
+#include <gsl/gsl_linalg.h>
+#endif
+
+void IssmSolve(IssmVec<IssmDouble>** pout,IssmMat<IssmDouble>* Kffin, IssmVec<IssmDouble>* pfin, Parameters* parameters){/*{{{*/
+
+	/*First off, we assume that the type of IssmVec is IssmSeqVec and the type of IssmMat is IssmDenseMat. So downcast: */
+	IssmSeqVec<IssmDouble>* pf=(IssmSeqVec<IssmDouble>*)pfin->vector;
+	IssmDenseMat<IssmDouble>* Kff=(IssmDenseMat<IssmDouble>*)Kffin->matrix;
+
+#ifdef _HAVE_GSL_
+	/*Intermediary: */
+	int M,N,N2;
+	IssmSeqVec<IssmDouble> *uf = NULL;
+
+	Kff->GetSize(&M,&N);
+	pf->GetSize(&N2);
+
+	if(N!=N2)_error_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !");
+	if(M!=N)_error_("Stiffness matrix should be square!");
+#ifdef _HAVE_ADOLC_
+	ensureContiguousLocations(N);
+#endif
+	IssmDouble *x  = xNew<IssmDouble>(N);
+#ifdef _HAVE_ADOLC_
+	SolverxSeq(x,Kff->matrix,pf->vector,N,parameters);
+#else
+	SolverxSeq(x,Kff->matrix,pf->vector,N);
+#endif
+
+	uf=new IssmSeqVec<IssmDouble>(x,N);
+	xDelete(x);
+
+	/*Assign output pointers:*/
+	IssmVec<IssmDouble>* out=new IssmVec<IssmDouble>();
+	out->vector=(IssmAbsVec<IssmDouble>*)uf;
+	*pout=out;
+
+#else
+	_error_("GSL support not compiled in!");
+#endif
+
+}/*}}}*/
+void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
+
+	/*Allocate output*/
+	double* X = xNew<double>(n); 
+
+	/*Solve*/
+	SolverxSeq(X,A,B,n);
+
+	/*Assign output pointer*/
+	*pX=X;
+}
+/*}}}*/
+void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
+#ifdef _HAVE_GSL_
+	/*GSL Matrices and vectors: */
+	int              s;
+	gsl_matrix_view  a;
+	gsl_vector_view  b,x;
+	gsl_permutation *p = NULL;
+//	for (int i=0; i<n*n; ++i) std::cout << "SolverxSeq A["<< i << "]=" << A[i] << std::endl;
+//	for (int i=0; i<n; ++i) std::cout << "SolverxSeq b["<< i << "]=" << B[i] << std::endl;
+	/*A will be modified by LU decomposition. Use copy*/
+	double* Acopy = xNew<double>(n*n);
+	xMemCpy(Acopy,A,n*n);
+
+	/*Initialize gsl matrices and vectors: */
+	a = gsl_matrix_view_array (Acopy,n,n);
+	b = gsl_vector_view_array (B,n);
+	x = gsl_vector_view_array (X,n);
+
+	/*Run LU and solve: */
+	p = gsl_permutation_alloc (n);
+	gsl_linalg_LU_decomp (&a.matrix, p, &s);
+	gsl_linalg_LU_solve (&a.matrix, p, &b.vector, &x.vector);
+
+	/*Clean up and assign output pointer*/
+	xDelete(Acopy);
+	gsl_permutation_free(p);
+#endif
+}
+/*}}}*/
+
+#ifdef _HAVE_ADOLC_
+int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y){ /*{{{*/
+	SolverxSeq(y,x, x+m*m, m); // x is where the matrix starts, x+m*m is where the right-hand side starts
+	return 0;
+} /*}}}*/
+int EDF_fos_forward_for_solverx(int n, IssmPDouble *inVal, IssmPDouble *inDeriv, int m, IssmPDouble *outVal, IssmPDouble *outDeriv) { /*{{{*/
+#ifdef _HAVE_GSL_
+	//  for (int i=0; i<m*m; ++i) std::cout << "EDF_fos_forward_for_solverx A["<< i << "]=" << inVal[i] << std::endl;
+	//  for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx b["<< i << "]=" << inVal[i+m*m] << std::endl;
+	// the matrix will be modified by LU decomposition. Use gsl_A copy
+	double* Acopy = xNew<double>(m*m);
+	xMemCpy(Acopy,inVal,m*m);
+	/*Initialize gsl matrices and vectors: */
+	gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
+	gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
+	gsl_permutation *perm_p = gsl_permutation_alloc (m);
+	int  signPerm;
+	// factorize
+	gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
+	gsl_vector *gsl_x_p = gsl_vector_alloc (m);
+	// solve for the value
+	gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
+	/*Copy result*/
+	xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
+	gsl_vector_free(gsl_x_p);
+	//  for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx x["<< i << "]=" << outVal[i] << std::endl;
+	// solve for the derivatives acc. to A * dx = r  with r=db - dA * x
+	// compute the RHS
+	double* r=xNew<double>(m);
+	for (int i=0; i<m; i++) {
+		r[i]=inDeriv[m*m+i]; // this is db[i]
+		for (int j=0;j<m; j++) {
+			r[i]-=inDeriv[i*m+j]*outVal[j]; // this is dA[i][j]*x[j]
+		}
+	}
+	gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
+	gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
+	gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
+	xMemCpy(outDeriv,gsl_vector_ptr(gsl_dx_p,0),m);
+	gsl_vector_free(gsl_dx_p);
+	xDelete(r);
+	gsl_permutation_free(perm_p);
+	xDelete(Acopy);
+#endif
+	return 0;
+} /*}}}*/
+int EDF_fov_forward_for_solverx(int n, IssmPDouble *inVal, int directionCount, IssmPDouble **inDeriv, int m, IssmPDouble *outVal, IssmPDouble **outDeriv) { /*{{{*/
+#ifdef _HAVE_GSL_
+	// the matrix will be modified by LU decomposition. Use gsl_A copy
+	double* Acopy = xNew<double>(m*m);
+	xMemCpy(Acopy,inVal,m*m);
+	/*Initialize gsl matrices and vectors: */
+	gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
+	gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
+	gsl_permutation *perm_p = gsl_permutation_alloc (m);
+	int  signPerm;
+	// factorize
+	gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
+	gsl_vector *gsl_x_p = gsl_vector_alloc (m);
+	// solve for the value
+	gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
+	/*Copy result*/
+	xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
+	gsl_vector_free(gsl_x_p);
+	// solve for the derivatives acc. to A * dx = r  with r=db - dA * x
+	double* r=xNew<double>(m);
+	gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
+	for (int dir=0;dir<directionCount;++dir) {
+		// compute the RHS
+		for (int i=0; i<m; i++) {
+			r[i]=inDeriv[m*m+i][dir]; // this is db[i]
+			for (int j=0;j<m; j++) {
+				r[i]-=inDeriv[i*m+j][dir]*outVal[j]; // this is dA[i][j]*x[j]
+			}
+		}
+		gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
+		gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
+		// reuse r
+		xMemCpy(r,gsl_vector_ptr(gsl_dx_p,0),m);
+		for (int i=0; i<m; i++) {
+			outDeriv[i][dir]=r[i];
+		}
+	}
+	gsl_vector_free(gsl_dx_p);
+	xDelete(r);
+	gsl_permutation_free(perm_p);
+	xDelete(Acopy);
+#endif
+	return 0;
+}
+/*}}}*/
+int EDF_fos_reverse_for_solverx(int m, double *dp_U, int n, double *dp_Z, double* dp_x, double* dp_y) { /*{{{*/
+	// copy to transpose the matrix
+	double* transposed=xNew<double>(m*m);
+	for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j];
+	gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m);
+	// the adjoint of the solution is our right-hand side
+	gsl_vector_view x_bar=gsl_vector_view_array(dp_U,m);
+	// the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute:
+	gsl_vector_view b_bar=gsl_vector_view_array(dp_Z+m*m,m);
+	gsl_permutation *perm_p = gsl_permutation_alloc (m);
+	int permSign;
+	gsl_linalg_LU_decomp (&aTransposed.matrix, perm_p, &permSign);
+	gsl_linalg_LU_solve (&aTransposed.matrix, perm_p, &x_bar.vector, &b_bar.vector);
+	// now do the adjoint of the matrix
+	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];
+	gsl_permutation_free(perm_p);
+	xDelete(transposed);
+	return 0;
+}
+/*}}}*/
+int EDF_fov_reverse_for_solverx(int m, int p, double **dpp_U, int n, double **dpp_Z, double* dp_x, double* dp_y) { /*{{{*/
+	// copy to transpose the matrix
+	double* transposed=xNew<double>(m*m);
+	for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j];
+	gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m);
+	gsl_permutation *perm_p = gsl_permutation_alloc (m);
+	int permSign;
+	gsl_linalg_LU_decomp (&aTransposed.matrix, perm_p, &permSign);
+	for (int weightsRowIndex=0;weightsRowIndex<p;++weightsRowIndex) {
+		// the adjoint of the solution is our right-hand side
+		gsl_vector_view x_bar=gsl_vector_view_array(dpp_U[weightsRowIndex],m);
+		// the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute:
+		gsl_vector_view b_bar=gsl_vector_view_array(dpp_Z[weightsRowIndex]+m*m,m);
+		gsl_linalg_LU_solve (&aTransposed.matrix, perm_p, &x_bar.vector, &b_bar.vector);
+		// now do the adjoint of the matrix
+		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];
+	}
+	gsl_permutation_free(perm_p);
+	xDelete(transposed);
+	return 0;
+}
+/*}}}*/
+void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
+	// pack inputs to conform to the EDF-prescribed interface
+        // ensure a contiguous block of locations:
+        ensureContiguousLocations(n*(n+1));
+        IssmDouble*  adoubleEDFin=xNew<IssmDouble>(n*(n+1));  // packed inputs, i.e. matrix and right hand side
+        for(int i=0; i<n*n;i++)adoubleEDFin[i]    =A[i];      // pack matrix
+        for(int i=0; i<n;  i++)adoubleEDFin[i+n*n]=B[i];      // pack the right hand side
+        IssmPDouble* pdoubleEDFin=xNew<IssmPDouble>(n*(n+1)); // provide space to transfer inputs during call_ext_fct
+	IssmPDouble* pdoubleEDFout=xNew<IssmPDouble>(n);           // provide space to transfer outputs during call_ext_fct
+	// call the wrapped solver through the registry entry we retrieve from parameters
+	call_ext_fct(dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p,
+	             n*(n+1), pdoubleEDFin, adoubleEDFin,
+	             n, pdoubleEDFout,X);
+	// for(int i=0; i<n;  i++) {ADOLC_DUMP_MACRO(X[i]);}
+	xDelete(adoubleEDFin);
+	xDelete(pdoubleEDFin);
+	xDelete(pdoubleEDFout);
+}
+/*}}}*/
+#endif
Index: ../trunk-jpl/src/c/toolkits/issm/issmtoolkit.h
===================================================================
--- ../trunk-jpl/src/c/toolkits/issm/issmtoolkit.h	(revision 14791)
+++ ../trunk-jpl/src/c/toolkits/issm/issmtoolkit.h	(revision 14792)
@@ -17,6 +17,7 @@
 #include "./IssmMat.h"
 #include "./IssmSeqVec.h"
 #include "./IssmVec.h"
+#include "./IssmSolver.h"
 
 #ifdef _HAVE_MPI_
 #include "./IssmMpiDenseMat.h"
Index: ../trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h
===================================================================
--- ../trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h	(revision 14791)
+++ ../trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h	(revision 14792)
@@ -271,27 +271,80 @@
 		/*}}}*/
 		/*FUNCTION Norm{{{*/
 		doubletype Norm(NormMode mode){
-			_error_("not supported yet!");
+			
+			
+			doubletype norm,local_norm;
+			doubletype absolute;
+			int i,j;
+
+			switch(mode){
+				case NORM_INF:
+					local_norm=0;
+					for(i=0;i<this->M;i++){
+						absolute=0;
+						for(j=0;j<this->N;j++){
+							absolute+=fabs(this->matrix[N*i+j]);
+						}
+						local_norm=max(local_norm,absolute);
+					}
+					MPI_Reduce(&local_norm, &norm, 1, MPI_DOUBLE, MPI_MAX, 0, IssmComm::GetComm());
+					MPI_Bcast(&norm,1,MPI_DOUBLE,0,IssmComm::GetComm());
+					return norm;
+					break; 
+				default:
+					_error_("unknown norm !");
+					break;
+			}
 		}
 		/*}}}*/
 		/*FUNCTION GetSize{{{*/
 		void GetSize(int* pM,int* pN){
-			_error_("not supported yet!");
+			*pM=M;
+			*pN=N;
 		}
 		/*}}}*/
 		/*FUNCTION GetLocalSize{{{*/
 		void GetLocalSize(int* pM,int* pN){
-			_error_("not supported yet!");
+			*pM=m;
+			*pN=N;
 		}
 		/*}}}*/
 		/*FUNCTION MatMult{{{*/
 		void MatMult(IssmAbsVec<doubletype>* Xin,IssmAbsVec<doubletype>* AXin){
-			_error_("not supported yet!");
+
+
+			int         i,j;
+			doubletype *X_serial  = NULL;
+
+
+			/*A check on the types: */
+			if(IssmVecTypeFromToolkitOptions()!=MpiEnum)_error_("MatMult operation only possible with 'mpi' vectors");
+
+			/*Now that we are sure, cast vectors: */
+			IssmMpiVec<doubletype>* X=(IssmMpiVec<doubletype>*)Xin;
+			IssmMpiVec<doubletype>* AX=(IssmMpiVec<doubletype>*)AXin;
+
+			/*Serialize input Xin: */
+			X_serial=X->ToMPISerial();
+
+			/*Every cpu has a serial version of the input vector. Use it to do the Matrix-Vector multiply 
+			 *locally and plug it into AXin: */
+			for(i=0;i<this->m;i++){
+				for(j=0;j<this->N;j++){
+					AX->vector[i]+=this->matrix[i*N+j]*X_serial[j];
+				}
+			}
+
+			/*Free ressources: */
+			xDelete<doubletype>(X_serial);
 		}
 		/*}}}*/
 		/*FUNCTION Duplicate{{{*/
 		IssmMpiDenseMat<doubletype>* Duplicate(void){
-			_error_("not supported yet!");
+
+			IssmMpiDenseMat<doubletype>* dup=new IssmMpiDenseMat<doubletype>(this->matrix,this->M,this->N,0);
+			return dup;
+
 		}
 		/*}}}*/
 		/*FUNCTION ToSerial{{{*/
@@ -313,7 +366,7 @@
 		/*}}}*/
 		/*FUNCTION Convert{{{*/
 		void Convert(MatrixType type){
-			/*do nothing: */
+			_error_("not supported yet!");
 		}
 		/*}}}*/		
 };
Index: ../trunk-jpl/src/c/toolkits/issm/IssmSolver.h
===================================================================
--- ../trunk-jpl/src/c/toolkits/issm/IssmSolver.h	(revision 0)
+++ ../trunk-jpl/src/c/toolkits/issm/IssmSolver.h	(revision 14792)
@@ -0,0 +1,38 @@
+/*!\file:  IssmSolver.h
+ * \brief main hook up from Solver toolkit object (in src/c/classes/toolkits) to the ISSM toolkit
+ */ 
+
+#ifndef _ISSM_SOLVER_H_
+#define _ISSM_SOLVER_H_
+
+/*Headers:*/
+/*{{{*/
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "../../include/types.h"
+
+/*}}}*/
+
+template <class doubletype> class IssmVec;
+template <class doubletype> class IssmMat;
+class Parameters;
+
+void IssmSolve(IssmVec<IssmDouble>** puf,IssmMat<IssmDouble>* Kff, IssmVec<IssmDouble>* pf,Parameters* parameters);
+void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n);
+void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n);
+
+#if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_)
+void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters);
+// call back functions:
+ADOLC_ext_fct EDF_for_solverx;
+ADOLC_ext_fct_fos_forward EDF_fos_forward_for_solverx;
+ADOLC_ext_fct_fos_reverse EDF_fos_reverse_for_solverx;
+ADOLC_ext_fct_fov_forward EDF_fov_forward_for_solverx;
+ADOLC_ext_fct_fov_reverse EDF_fov_reverse_for_solverx;
+#endif
+
+#endif 
Index: ../trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.h
===================================================================
--- ../trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.h	(revision 0)
+++ ../trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.h	(revision 14792)
@@ -0,0 +1,24 @@
+/*!\file:  PetscMat.h
+ */ 
+
+#ifndef _PETSC_SOLVER_H_
+#define _PETSC_SOLVER_H_
+
+/*Headers:*/
+/*{{{*/
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "../petscincludes.h"
+class Parameters;
+
+/*}}}*/
+
+void	PetscSolve(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters);
+void	SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters);
+void    DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum);
+
+#endif //#ifndef _PETSCSOLVER_H_
Index: ../trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.cpp
===================================================================
--- ../trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.cpp	(revision 0)
+++ ../trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.cpp	(revision 14792)
@@ -0,0 +1,232 @@
+/*!\file SolverxPetsc
+ * \brief Petsc implementation of solver
+ */
+
+
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "./PetscSolver.h"
+#include "../../../shared/Numerics/Verbosity.h"
+#include "../../../shared/Alloc/alloc.h"
+#include "../../../shared/Exceptions/exceptions.h"
+#include "../../../classes/IssmComm.h"
+#include "../../../EnumDefinitions/EnumDefinitions.h"
+
+void	PetscSolve(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters){ /*{{{*/
+
+	PetscVec* uf=new PetscVec();
+
+	Vec uf0_vector = NULL;
+	Vec df_vector  = NULL;
+
+	if(uf0) uf0_vector = uf0->vector;
+	if(df)  df_vector  = df->vector;
+
+	SolverxPetsc(&uf->vector, Kff->matrix, pf->vector, uf0_vector, df_vector, parameters);
+
+	*puf=uf;
+
+}
+/*}}}*/
+void	SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){ /*{{{*/
+
+	/*Output: */
+	Vec        uf               = NULL;
+
+	/*Intermediary: */
+	int        local_m,local_n,global_m,global_n;
+
+	/*Solver */
+	KSP        ksp              = NULL;
+	PC         pc               = NULL;
+	int        iteration_number;
+	int        solver_type;
+	bool       fromlocalsize    = true;
+	#if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
+	PetscTruth flag,flg;
+	#else
+	PetscBool flag,flg;
+	#endif
+
+	/*Stokes: */
+	IS         isv=NULL;
+	IS         isp=NULL;
+
+	#if _PETSC_MAJOR_ >= 3 
+	char ksp_type[50];
+	#endif
+
+	/*Display message*/
+	#if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
+	if(VerboseSolver())PetscOptionsPrint(stdout);
+	#else
+	if(VerboseSolver())PetscOptionsView(PETSC_VIEWER_STDOUT_WORLD);
+	#endif
+
+	/*First, check that f-set is not NULL, i.e. model is fully constrained:*/ 
+	_assert_(Kff);
+	MatGetSize(Kff,&global_m,&global_n); _assert_(global_m==global_n);
+	if(!global_n){
+		*puf=NewVec(0,IssmComm::GetComm()); return;
+	}
+
+	/*Initial guess */
+	/*Now, check that we are not giving an initial guess to the solver, if we are running a direct solver: */
+	#if _PETSC_MAJOR_ >= 3 
+	PetscOptionsGetString(PETSC_NULL,"-ksp_type",ksp_type,49,&flg);
+	if (strcmp(ksp_type,"preonly")==0)uf0=NULL;
+	#endif
+
+	/*If initial guess for the solution exists, use it to create uf, otherwise, 
+	 * duplicate the right hand side so that the solution vector has the same structure*/
+	if(uf0){
+		VecDuplicate(uf0,&uf); VecCopy(uf0,uf);
+	}
+	else{
+		MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,IssmComm::GetComm(),fromlocalsize);
+	}
+
+	/*Process petsc options to see if we are using special types of external solvers*/
+	PetscOptionsDetermineSolverType(&solver_type);
+
+	/*Check the solver is available*/
+	if(solver_type==MUMPSPACKAGE_LU || solver_type==MUMPSPACKAGE_CHOL){
+		#if _PETSC_MAJOR_ >=3
+			#ifndef _HAVE_MUMPS_
+			_error_("requested MUMPS solver, which was not compiled into ISSM!\n");
+			#endif
+		#endif
+	}
+
+	/*Prepare solver*/
+	KSPCreate(IssmComm::GetComm(),&ksp);
+	KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN);
+	KSPSetFromOptions(ksp);
+
+	#if _PETSC_MAJOR_==3
+	/*Specific solver?: */
+	KSPGetPC(ksp,&pc);
+	if (solver_type==MUMPSPACKAGE_LU){
+		#if _PETSC_MINOR_==1
+		PCFactorSetMatSolverPackage(pc,MAT_SOLVER_MUMPS);
+		#else
+		PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
+		#endif
+	}
+
+	/*Stokes: */
+	if (solver_type==StokesSolverEnum){
+		/*Make indices out of doftypes: */
+		if(!df)_error_("need doftypes for Stokes solver!\n");
+		DofTypesToIndexSet(&isv,&isp,df,StokesSolverEnum);
+
+		/*Set field splits: */
+		KSPGetPC(ksp,&pc);
+		#if _PETSC_MINOR_==1
+		PCFieldSplitSetIS(pc,isv);
+		PCFieldSplitSetIS(pc,isp);
+		#else
+		PCFieldSplitSetIS(pc,PETSC_NULL,isv);
+		PCFieldSplitSetIS(pc,PETSC_NULL,isp);
+		#endif
+
+	}
+	#endif
+
+	/*If there is an initial guess for the solution, use it
+	 * except if we are using the MUMPS direct solver
+	 * where any initial guess will crash Petsc*/
+	if (uf0){
+		if((solver_type!=MUMPSPACKAGE_LU) && (solver_type!=MUMPSPACKAGE_CHOL) && (solver_type!=SPOOLESPACKAGE_LU) && (solver_type!=SPOOLESPACKAGE_CHOL) && (solver_type!=SUPERLUDISTPACKAGE)){
+			KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
+		}
+	}
+
+	/*Solve: */
+	if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
+	KSPSolve(ksp,pf,uf);
+
+	/*Check convergence*/
+	KSPGetIterationNumber(ksp,&iteration_number);
+	if (iteration_number<0) _error_("Solver diverged at iteration number: " << -iteration_number);
+
+	/*Free resources:*/
+	KSPFree(&ksp);
+
+	/*Assign output pointers:*/
+	*puf=uf;
+} 
+/*}}}*/
+void DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum){ /*{{{*/
+
+	/*output: */
+	IS          isv=NULL;
+	IS          isp=NULL;
+
+	int         start,end;
+	IssmDouble*     df_local=NULL;
+	int         df_local_size;
+	int         i;
+
+	int*     pressure_indices=NULL;
+	int*     velocity_indices=NULL;
+	int      pressure_num=0;
+	int      velocity_num=0;
+	int      pressure_count=0;
+	int      velocity_count=0;
+
+	if(typeenum==StokesSolverEnum){
+
+		/*Ok, recover doftypes vector values and indices: */
+		VecGetOwnershipRange(df,&start,&end);
+		VecGetLocalSize(df,&df_local_size);
+		VecGetArray(df,&df_local);
+
+		pressure_num=0;
+		velocity_num=0;
+		for(i=0;i<df_local_size;i++){
+			if (df_local[i]==PressureEnum)pressure_num++;
+			else velocity_num++;
+		}
+
+		/*Allocate indices: */
+		if(pressure_num)pressure_indices=xNew<int>(pressure_num);
+		if(velocity_num)velocity_indices=xNew<int>(velocity_num);
+
+		pressure_count=0;
+		velocity_count=0;
+		for(i=0;i<df_local_size;i++){
+			if (df_local[i]==PressureEnum){
+				pressure_indices[pressure_count]=start+i;
+				pressure_count++;
+			}
+			if (df_local[i]==VelocityEnum){
+				velocity_indices[velocity_count]=start+i;
+				velocity_count++;
+			}
+		}
+		VecRestoreArray(df,&df_local);
+
+		/*Create indices sets: */
+		#if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
+		ISCreateGeneral(IssmComm::GetComm(),pressure_num,pressure_indices,&isp);
+		ISCreateGeneral(IssmComm::GetComm(),velocity_num,velocity_indices,&isv);
+		#else
+		ISCreateGeneral(IssmComm::GetComm(),pressure_num,pressure_indices,PETSC_COPY_VALUES,&isp);
+		ISCreateGeneral(IssmComm::GetComm(),velocity_num,velocity_indices,PETSC_COPY_VALUES,&isv);
+		#endif
+	}
+
+	/*Free ressources:*/
+	xDelete<int>(pressure_indices);
+	xDelete<int>(velocity_indices);
+
+	/*Assign output pointers:*/
+	*pisv=isv;
+	*pisp=isp;
+}
+/*}}}*/
Index: ../trunk-jpl/src/c/toolkits/petsc/objects/petscobjects.h
===================================================================
--- ../trunk-jpl/src/c/toolkits/petsc/objects/petscobjects.h	(revision 14791)
+++ ../trunk-jpl/src/c/toolkits/petsc/objects/petscobjects.h	(revision 14792)
@@ -7,5 +7,6 @@
 
 #include "./PetscMat.h"
 #include "./PetscVec.h"
+#include "./PetscSolver.h"
 
 #endif
Index: ../trunk-jpl/src/c/Makefile.am
===================================================================
--- ../trunk-jpl/src/c/Makefile.am	(revision 14791)
+++ ../trunk-jpl/src/c/Makefile.am	(revision 14792)
@@ -118,8 +118,10 @@
 					./classes/matrix/ElementMatrix.cpp\
 					./classes/matrix/ElementVector.h\
 					./classes/matrix/ElementVector.cpp\
-					./classes/matrix/Matrix.h\
-					./classes/matrix/Vector.h\
+					./classes/toolkits/toolkitobjects.h\
+					./classes/toolkits/Matrix.h\
+					./classes/toolkits/Vector.h\
+					./classes/toolkits/Solver.h\
 					./classes/objects/Params/Param.h\
 					./classes/objects/Params/GenericParam.h\
 					./classes/objects/Params/BoolParam.cpp\
@@ -228,6 +230,8 @@
 					./toolkits/issm/IssmMat.h\
 					./toolkits/issm/IssmSeqVec.h\
 					./toolkits/issm/IssmVec.h\
+					./toolkits/issm/IssmSolver.h\
+					./toolkits/issm/IssmSolver.cpp\
 					./toolkits/triangle/triangleincludes.h\
 					./toolkitsenums.h\
 					./toolkits.h\
@@ -319,7 +323,6 @@
 					./modules/ResetCoordinateSystemx/ResetCoordinateSystemx.cpp\
 					./modules/Solverx/Solverx.cpp\
 					./modules/Solverx/Solverx.h\
-					./modules/Solverx/SolverxSeq.cpp\
 					./modules/VecMergex/VecMergex.cpp\
 					./modules/VecMergex/VecMergex.h\
 					./modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp\
@@ -759,9 +762,9 @@
 					./toolkits/petsc/objects/PetscMat.cpp\
 					./toolkits/petsc/objects/PetscVec.h\
 					./toolkits/petsc/objects/PetscVec.cpp\
-					./toolkits/petsc/petscincludes.h\
-					./modules/Solverx/SolverxPetsc.cpp\
-					./modules/Solverx/DofTypesToIndexSet.cpp
+					./toolkits/petsc/objects/PetscSolver.cpp\
+					./toolkits/petsc/objects/PetscSolver.h\
+					./toolkits/petsc/petscincludes.h
 
 #}}}
 #Mpi sources  {{{
Index: ../trunk-jpl/src/c/classes/classes.h
===================================================================
--- ../trunk-jpl/src/c/classes/classes.h	(revision 14791)
+++ ../trunk-jpl/src/c/classes/classes.h	(revision 14792)
@@ -11,6 +11,9 @@
 /*matrix: */
 #include "./matrix/matrixobjects.h"
 
+/*toolkit objects: */
+#include "./toolkits/toolkitobjects.h"
+
 /*bamg: */
 #include "./bamg/bamgobjects.h"
 
Index: ../trunk-jpl/src/c/classes/matrix/Vector.h
===================================================================
--- ../trunk-jpl/src/c/classes/matrix/Vector.h	(revision 14791)
+++ ../trunk-jpl/src/c/classes/matrix/Vector.h	(revision 14792)
@@ -1,357 +0,0 @@
-/*!\file:  Vector.h
- * \brief wrapper to vector objects. The goal is to control which API (PETSc,Scalpack, Plapack?) 
- * implements our underlying vector format.
- */ 
-
-#ifndef _VECTOR_H_
-#define _VECTOR_H_
-
-/*Headers:*/
-/*{{{*/
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-#include <cstring>
-#include "../../toolkits/toolkits.h"
-#include "../../EnumDefinitions/EnumDefinitions.h"
-/*}}}*/
-
-enum vectortype { PetscVecType, IssmVecType };
-
-template <class doubletype> 
-class Vector{
-
-	public:
-
-		int  type;
-		#ifdef _HAVE_PETSC_
-		PetscVec* pvector;
-		#endif
-		IssmVec<doubletype>* ivector; 
-
-		/*Vector constructors, destructors */
-		Vector(){ /*{{{*/
-
-			InitCheckAndSetType();
-		}
-		/*}}}*/
-		Vector(int M,bool fromlocalsize=false){ /*{{{*/
-			
-			InitCheckAndSetType();
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				this->pvector=new PetscVec(M,fromlocalsize);
-				#endif
-			}
-			else this->ivector=new IssmVec<doubletype>(M,fromlocalsize);
-
-		}
-		/*}}}*/
-		Vector(int m,int M){ /*{{{*/
-
-			InitCheckAndSetType();
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-					this->pvector=new PetscVec(m,M);
-				 #endif
-			}
-			else this->ivector=new IssmVec<doubletype>(m,M);
-		}
-		/*}}}*/
-		Vector(doubletype* serial_vec,int M){ /*{{{*/
-			
-			InitCheckAndSetType();
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				this->pvector=new PetscVec(serial_vec,M);
-				#endif
-			}
-			else this->ivector=new IssmVec<doubletype>(serial_vec,M);
-		}
-		/*}}}*/
-		~Vector(){ /*{{{*/
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				delete this->pvector;
-				#endif
-			}
-			else delete this->ivector;
-		}
-		/*}}}*/
-		#ifdef _HAVE_PETSC_
-		Vector(Vec petsc_vector){ /*{{{*/
-
-			this->type=PetscVecType;
-			this->ivector=NULL;
-			this->pvector=new PetscVec(petsc_vector);
-
-		}
-		/*}}}*/
-		#endif
-		void InitCheckAndSetType(void){ /*{{{*/
-
-			char* toolkittype=NULL;
-
-			#ifdef _HAVE_PETSC_
-			pvector=NULL;
-			#endif
-			ivector=NULL;
-
-			/*retrieve toolkittype: */
-			toolkittype=ToolkitOptions::GetToolkitType();
-
-			/*set vector type: */
-			if (strcmp(toolkittype,"petsc")==0){
-				#ifdef _HAVE_PETSC_
-				type=PetscVecType; 
-				#else
-				_error_("cannot create petsc vector without PETSC compiled!");
-				#endif
-			}
-			else if(strcmp(toolkittype,"issm")==0){
-				/*let this choice stand:*/
-				type=IssmVecType;
-			}
-			else _error_("unknow toolkit type ");
-			
-			/*Free ressources: */
-			xDelete<char>(toolkittype);
-		}
-		/*}}}*/
-
-		/*Vector specific routines*/
-		/*FUNCTION Echo{{{*/
-		void Echo(void){_assert_(this);
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				this->pvector->Echo();
-				#endif
-			}
-			else this->ivector->Echo();
-
-		}
-		/*}}}*/
-		/*FUNCTION Assemble{{{*/
-		void Assemble(void){_assert_(this);
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				this->pvector->Assemble();
-				#endif
-			}
-			else this->ivector->Assemble();
-
-		}
-		/*}}}*/
-		/*FUNCTION SetValues{{{*/
-		void SetValues(int ssize, int* list, doubletype* values, InsMode mode){ _assert_(this);
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				this->pvector->SetValues(ssize,list,values,mode);
-				#endif
-			}
-			else this->ivector->SetValues(ssize,list,values,mode);
-
-		}
-		/*}}}*/
-		/*FUNCTION SetValue{{{*/
-		void SetValue(int dof, doubletype value, InsMode mode){_assert_(this);
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				this->pvector->SetValue(dof,value,mode);
-				#endif
-			}
-			else this->ivector->SetValue(dof,value,mode);
-
-		}
-		/*}}}*/
-		/*FUNCTION GetValue{{{*/
-		void GetValue(doubletype* pvalue,int dof){_assert_(this);
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				this->pvector->GetValue(pvalue,dof);
-				#endif
-			}
-			else this->ivector->GetValue(pvalue,dof);
-
-		}
-		/*}}}*/
-		/*FUNCTION GetSize{{{*/
-		void GetSize(int* pM){_assert_(this);
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				this->pvector->GetSize(pM);
-				#endif
-			}
-			else this->ivector->GetSize(pM);
-
-		}
-		/*}}}*/
-		/*FUNCTION IsEmpty{{{*/
-		bool IsEmpty(void){
-			int M;
-			
-			_assert_(this);
-			this->GetSize(&M);
-
-			if(M==0) 
-				return true;
-			else
-				return false;
-		}
-		/*}}}*/
-		/*FUNCTION GetLocalSize{{{*/
-		void GetLocalSize(int* pM){_assert_(this);
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				this->pvector->GetLocalSize(pM);
-				#endif
-			}
-			else this->ivector->GetLocalSize(pM);
-
-		}
-		/*}}}*/
-		/*FUNCTION Duplicate{{{*/
-		Vector<doubletype>* Duplicate(void){_assert_(this);
-
-			Vector<doubletype>* output=NULL;
-				
-			output=new Vector<doubletype>();
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				output->pvector=this->pvector->Duplicate();
-				#endif
-			}
-			else output->ivector=this->ivector->Duplicate();
-
-			return output;
-
-		}
-		/*}}}*/
-		/*FUNCTION Set{{{*/
-		void Set(doubletype value){_assert_(this);
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				this->pvector->Set(value);
-				#endif
-			}
-			else this->ivector->Set(value);
-
-		}
-		/*}}}*/
-		/*FUNCTION AXPY{{{*/
-		void AXPY(Vector* X, doubletype a){_assert_(this);
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				this->pvector->AXPY(X->pvector,a);
-				#endif
-			}
-			else this->ivector->AXPY(X->ivector,a);
-
-		}
-		/*}}}*/
-		/*FUNCTION AYPX{{{*/
-		void AYPX(Vector* X, doubletype a){_assert_(this);
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				this->pvector->AYPX(X->pvector,a);
-				#endif
-			}
-			else this->ivector->AYPX(X->ivector,a);
-		}
-		/*}}}*/
-		/*FUNCTION ToMPISerial{{{*/
-		doubletype* ToMPISerial(void){
-
-			doubletype* vec_serial=NULL;
-			
-			_assert_(this);
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				vec_serial=this->pvector->ToMPISerial();
-				#endif
-			}
-			else vec_serial=this->ivector->ToMPISerial();
-
-			return vec_serial;
-
-		}
-		/*}}}*/
-		/*FUNCTION Copy{{{*/
-		void Copy(Vector* to){_assert_(this);
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				this->pvector->Copy(to->pvector);
-				#endif
-			}
-			else this->ivector->Copy(to->ivector);
-		}
-		/*}}}*/
-		/*FUNCTION Norm{{{*/
-		doubletype Norm(NormMode norm_type){_assert_(this);
-
-			doubletype norm=0;
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				norm=this->pvector->Norm(norm_type);
-				#endif
-			}
-			else norm=this->ivector->Norm(norm_type);
-			return norm;
-		}
-		/*}}}*/
-		/*FUNCTION Scale{{{*/
-		void Scale(doubletype scale_factor){_assert_(this);
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				this->pvector->Scale(scale_factor);
-				#endif
-			}
-			else this->ivector->Scale(scale_factor);
-		}
-		/*}}}*/
-		/*FUNCTION Dot{{{*/
-		doubletype Dot(Vector* vector){_assert_(this);
-
-			doubletype dot;
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				dot=this->pvector->Dot(vector->pvector);
-				#endif
-			}
-			else dot=this->ivector->Dot(vector->ivector);
-			return dot;
-		}
-		/*}}}*/
-		/*FUNCTION PointwiseDivide{{{*/
-		void PointwiseDivide(Vector* x,Vector* y){_assert_(this);
-
-			if(type==PetscVecType){
-				#ifdef _HAVE_PETSC_
-				this->pvector->PointwiseDivide(x->pvector,y->pvector);
-				#endif
-			}
-			else this->ivector->PointwiseDivide(x->ivector,y->ivector);
-		}
-		/*}}}*/
-};
-#endif //#ifndef _VECTOR_H_
Index: ../trunk-jpl/src/c/classes/matrix/Matrix.h
===================================================================
--- ../trunk-jpl/src/c/classes/matrix/Matrix.h	(revision 14791)
+++ ../trunk-jpl/src/c/classes/matrix/Matrix.h	(revision 14792)
@@ -1,328 +0,0 @@
-/*!\file:  Matrix.h
- * \brief wrapper to matrix objects. The goal is to control which API (PETSc,Scalpack, Plapack?) 
- * implements our underlying matrix format.
- */ 
-
-#ifndef _MATRIX_H_
-#define _MATRIX_H_
-
-/*Headers:*/
-/*{{{*/
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-#include <cstring>
-#include "../../toolkits/toolkits.h"
-#include "../../EnumDefinitions/EnumDefinitions.h"
-/*}}}*/
-
-enum matrixtype { PetscMatType, IssmMatType };
-
-template <class doubletype> class Vector;
-
-template <class doubletype> 
-class Matrix{
-
-	public:
-
-		int       type;
-		#ifdef _HAVE_PETSC_
-		PetscMat              *pmatrix;
-		#endif
-		IssmMat<doubletype>   *imatrix;
-
-		/*Matrix constructors, destructors*/
-		/*FUNCTION Matrix(){{{*/
-		Matrix(){
-			InitCheckAndSetType();
-		}
-		/*}}}*/
-		/*FUNCTION Matrix(int M,int N){{{*/
-		Matrix(int M,int N){
-			
-			InitCheckAndSetType();
-
-			if(type==PetscMatType){
-				#ifdef _HAVE_PETSC_
-				this->pmatrix=new PetscMat(M,N);
-				#endif
-			}
-			else{
-				this->imatrix=new IssmMat<doubletype>(M,N);
-			}
-
-		}
-		/*}}}*/
-		/*FUNCTION Matrix(int m,int n,int M,int N,int* d_nnz,int* o_nnz){{{*/
-		Matrix(int m,int n,int M,int N,int* d_nnz,int* o_nnz){
-
-			InitCheckAndSetType();
-
-			if(type==PetscMatType){
-				#ifdef _HAVE_PETSC_
-				this->pmatrix=new PetscMat(m,n,M,N,d_nnz,o_nnz);
-				#endif
-			}
-			else{
-				this->imatrix=new IssmMat<doubletype>(m,n,M,N,d_nnz,o_nnz);
-			}
-
-		}
-		/*}}}*/
-		/*FUNCTION Matrix(int M,int N,IssmDouble sparsity){{{*/
-		Matrix(int M,int N,double sparsity){
-			
-			InitCheckAndSetType();
-
-			if(type==PetscMatType){
-				#ifdef _HAVE_PETSC_
-				this->pmatrix=new PetscMat(M,N,sparsity);
-				#endif
-			}
-			else{
-				this->imatrix=new IssmMat<doubletype>(M,N,sparsity);
-			}
-		}
-		/*}}}*/
-		/*FUNCTION Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity){{{*/
-		Matrix(IssmPDouble* serial_mat,int M,int N,IssmPDouble sparsity){
-			
-			InitCheckAndSetType();
-
-			if(type==PetscMatType){
-				#ifdef _HAVE_PETSC_
-				this->pmatrix=new PetscMat(serial_mat,M,N,sparsity);
-				#endif
-			}
-			else{
-				this->imatrix=new IssmMat<doubletype>(serial_mat,M,N,sparsity);
-			}
-
-		}
-		/*}}}*/
-		/*FUNCTION Matrix(int M,int N,int connectivity,int numberofdofspernode){{{*/
-		Matrix(int M,int N,int connectivity,int numberofdofspernode){
-			
-			InitCheckAndSetType();
-
-			if(type==PetscMatType){
-				#ifdef _HAVE_PETSC_
-				this->pmatrix=new PetscMat(M,N,connectivity,numberofdofspernode);
-				#endif
-			}
-			else{
-				this->imatrix=new IssmMat<doubletype>(M,N,connectivity,numberofdofspernode);
-			}
-
-		}
-		/*}}}*/
-		/*FUNCTION ~Matrix(){{{*/
-		~Matrix(){
-
-			if(type==PetscMatType){
-				#ifdef _HAVE_PETSC_
-				delete this->pmatrix;
-				#endif
-			}
-			else delete this->imatrix;
-
-		}
-		/*}}}*/
-		/*FUNCTION InitCheckAndSetType(){{{*/
-		void InitCheckAndSetType(void){
-
-			char* toolkittype=NULL;
-
-			#ifdef _HAVE_PETSC_
-			pmatrix=NULL;
-			#endif
-			imatrix=NULL;
-
-			/*retrieve toolkittype: */
-			toolkittype=ToolkitOptions::GetToolkitType();
-
-			/*set matrix type: */
-			if (strcmp(toolkittype,"petsc")==0){
-				#ifdef _HAVE_PETSC_
-				type=PetscMatType; 
-				#else
-				_error_("cannot create petsc matrix without PETSC compiled!");
-				#endif
-			}
-			else if(strcmp(toolkittype,"issm")==0){
-				/*let this choice stand:*/
-				type=IssmMatType;
-			}
-			else _error_("unknow toolkit type ");
-			
-			/*Free ressources: */
-			xDelete<char>(toolkittype);
-		}
-		/*}}}*/
-
-		/*Matrix specific routines:*/
-		/*FUNCTION Echo{{{*/
-		void Echo(void){
-			_assert_(this);
-
-			if(type==PetscMatType){
-				#ifdef _HAVE_PETSC_
-				this->pmatrix->Echo();
-				#endif
-			}
-			else{
-				this->imatrix->Echo();
-			}
-
-		}
-		/*}}}*/
-		/*FUNCTION AllocationInfo{{{*/
-		void AllocationInfo(void){
-			_assert_(this);
-			if(type==PetscMatType){
-				#ifdef _HAVE_PETSC_
-				this->pmatrix->AllocationInfo();
-				#endif
-			}
-			else{
-				this->imatrix->AllocationInfo();
-			}
-		}/*}}}*/
-		/*FUNCTION Assemble{{{*/
-		void Assemble(void){
-
-			if(type==PetscMatType){
-				#ifdef _HAVE_PETSC_
-				this->pmatrix->Assemble();
-				#endif
-			}
-			else{
-				this->imatrix->Assemble();
-			}
-		}
-		/*}}}*/
-		/*FUNCTION Norm{{{*/
-		IssmDouble Norm(NormMode norm_type){
-
-			IssmDouble norm=0;
-
-			if(type==PetscMatType){
-				#ifdef _HAVE_PETSC_
-				norm=this->pmatrix->Norm(norm_type);
-				#endif
-			}
-			else{
-				norm=this->imatrix->Norm(norm_type);
-			}
-
-			return norm;
-		}
-		/*}}}*/
-		/*FUNCTION GetSize{{{*/
-		void GetSize(int* pM,int* pN){
-
-			if(type==PetscMatType){
-				#ifdef _HAVE_PETSC_
-				this->pmatrix->GetSize(pM,pN);
-				#endif
-			}
-			else{
-				this->imatrix->GetSize(pM,pN);
-			}
-
-		}
-		/*}}}*/
-		/*FUNCTION GetLocalSize{{{*/
-		void GetLocalSize(int* pM,int* pN){
-
-			if(type==PetscMatType){
-				#ifdef _HAVE_PETSC_
-				this->pmatrix->GetLocalSize(pM,pN);
-				#endif
-			}
-			else{
-				this->imatrix->GetLocalSize(pM,pN);
-			}
-
-		}
-		/*}}}*/
-		/*FUNCTION MatMult{{{*/
-		void MatMult(Vector<doubletype>* X,Vector<doubletype>* AX){
-
-			if(type==PetscMatType){
-				#ifdef _HAVE_PETSC_
-				this->pmatrix->MatMult(X->pvector,AX->pvector);
-				#endif
-			}
-			else{
-				this->imatrix->MatMult(X->ivector,AX->ivector);
-			}
-
-		}
-		/*}}}*/
-		/*FUNCTION Duplicate{{{*/
-		Matrix<doubletype>* Duplicate(void){
-
-			Matrix<doubletype>* output=new Matrix<doubletype>();
-
-			if(type==PetscMatType){
-				#ifdef _HAVE_PETSC_
-				output->pmatrix=this->pmatrix->Duplicate();
-				#endif
-			}
-			else{
-				output->imatrix=this->imatrix->Duplicate();
-			}
-
-			return output;
-		}
-		/*}}}*/
-		/*FUNCTION ToSerial{{{*/
-		doubletype* ToSerial(void){
-
-			doubletype* output=NULL;
-
-			if(type==PetscMatType){
-				#ifdef _HAVE_PETSC_
-				output=this->pmatrix->ToSerial();
-				#endif
-			}
-			else{
-				output=this->imatrix->ToSerial();
-			}
-
-			return output;
-		}
-		/*}}}*/
-		/*FUNCTION SetValues{{{*/
-		void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){
-
-			if(type==PetscMatType){
-				#ifdef _HAVE_PETSC_
-				this->pmatrix->SetValues(m,idxm,n,idxn,values,mode);
-				#endif
-			}
-			else{
-				this->imatrix->SetValues(m,idxm,n,idxn,values,mode);
-			}
-		}
-		/*}}}*/
-		/*FUNCTION Convert{{{*/
-		void Convert(MatrixType newtype){
-
-			if(type==PetscMatType){
-				#ifdef _HAVE_PETSC_
-				this->pmatrix->Convert(newtype);
-				#endif
-			}
-			else{
-				this->imatrix->Convert(newtype);
-			}
-
-		}
-		/*}}}*/
-};
-
-#endif //#ifndef _MATRIX_H_
Index: ../trunk-jpl/src/c/classes/matrix/matrixobjects.h
===================================================================
--- ../trunk-jpl/src/c/classes/matrix/matrixobjects.h	(revision 14791)
+++ ../trunk-jpl/src/c/classes/matrix/matrixobjects.h	(revision 14792)
@@ -8,7 +8,5 @@
 /*Numerics:*/
 #include "./ElementMatrix.h"
 #include "./ElementVector.h"
-#include "./Vector.h"
-#include "./Matrix.h"
 
 #endif
Index: ../trunk-jpl/src/c/classes/toolkits/Solver.h
===================================================================
--- ../trunk-jpl/src/c/classes/toolkits/Solver.h	(revision 0)
+++ ../trunk-jpl/src/c/classes/toolkits/Solver.h	(revision 14792)
@@ -0,0 +1,99 @@
+/*!\file:  Solver.h
+ */ 
+
+#ifndef _SOLVER_CLASS_H_
+#define _SOLVER_CLASS_H_
+
+/*Headers:*/
+/*{{{*/
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+#include "../../toolkits/toolkits.h"
+/*}}}*/
+
+#include "./Matrix.h"
+#include "./Vector.h"
+class Parameters;
+
+template <class doubletype> 
+class Solver{
+
+	private:
+		Matrix<doubletype>* Kff;
+		Vector<doubletype>* pf;
+		Vector<doubletype>* uf0;
+		Vector<doubletype>* df;
+		Parameters* parameters;
+
+	public:
+
+		/*Constructors, destructors:*/
+		/*Solver(){{{*/
+		Solver(){
+		}
+		/*}}}*/
+		/*Solver(Matrix<doubletype>* Kff, Vector<doubletype>* pf, Vector<doubletype>* uf0,Vector<doubletype>* df, Parameters* parameters):{{{*/
+		Solver(Matrix<doubletype>* Kff_in, Vector<doubletype>* pf_in, Vector<doubletype>* uf0_in,Vector<doubletype>* df_in, Parameters* parameters_in){
+
+			/*In debugging mode, check that stiffness matrix and load vectors are not NULL (they can be empty)*/
+			_assert_(Kff_in);
+			_assert_(pf_in);
+
+			/*initialize fields: */
+			this->Kff=Kff_in;
+			this->pf=pf_in;
+			this->uf0=uf0_in;
+			this->df=df_in;
+			this->parameters=parameters_in;
+		}
+		/*}}}*/
+		/*~Solver(){{{*/
+		~Solver(){
+		}
+		/*}}}*/
+		
+		/*Methods: */
+		/*Vector<doubletype Solve(void): {{{*/
+		Vector<doubletype>* Solve(){
+
+			/*output: */
+			Vector<doubletype>* uf=NULL;
+
+			/*Initialize vector: */
+			uf=new Vector<doubletype>();
+
+			/*According to matrix type, use specific solvers: */
+			switch(Kff->type){
+				#ifdef _HAVE_PETSC_
+				case PetscMatType:{
+					PetscVec* uf0_vector = NULL;
+					PetscVec* df_vector  = NULL;
+					if(uf0) uf0_vector = uf0->pvector;
+					if(df)  df_vector  = df->pvector;
+					PetscSolve(&uf->pvector,Kff->pmatrix,pf->pvector,uf0_vector,df_vector,parameters);
+					break;
+								  }
+				#endif
+				case IssmMatType:{
+					IssmSolve(&uf->ivector,Kff->imatrix,pf->ivector,parameters);
+					break;
+								 }
+				default:
+					_error_("Matrix type: " << Kff->type << " not supported yet!");
+			}
+
+			/*allocate output pointer: */
+			return uf;
+
+		}
+		/*}}}*/
+
+};
+
+
+
+
+#endif //#ifndef _SOLVER_H_
Index: ../trunk-jpl/src/c/classes/toolkits/toolkitobjects.h
===================================================================
--- ../trunk-jpl/src/c/classes/toolkits/toolkitobjects.h	(revision 0)
+++ ../trunk-jpl/src/c/classes/toolkits/toolkitobjects.h	(revision 14792)
@@ -0,0 +1,13 @@
+
+/*!\file:  toolkitobjects.h
+ * \brief wrappers to toolkits
+ */ 
+
+#ifndef _TOOLKIT_OBJECTS_H_
+#define _TOOLKIT_OBJECTS_H_
+
+#include "./Vector.h"
+#include "./Matrix.h"
+#include "./Solver.h"
+
+#endif
Index: ../trunk-jpl/src/c/classes/toolkits/Vector.h
===================================================================
--- ../trunk-jpl/src/c/classes/toolkits/Vector.h	(revision 0)
+++ ../trunk-jpl/src/c/classes/toolkits/Vector.h	(revision 14792)
@@ -0,0 +1,357 @@
+/*!\file:  Vector.h
+ * \brief wrapper to vector objects. The goal is to control which API (PETSc,Scalpack, Plapack?) 
+ * implements our underlying vector format.
+ */ 
+
+#ifndef _VECTOR_H_
+#define _VECTOR_H_
+
+/*Headers:*/
+/*{{{*/
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+#include <cstring>
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+/*}}}*/
+
+enum vectortype { PetscVecType, IssmVecType };
+
+template <class doubletype> 
+class Vector{
+
+	public:
+
+		int  type;
+		#ifdef _HAVE_PETSC_
+		PetscVec* pvector;
+		#endif
+		IssmVec<doubletype>* ivector; 
+
+		/*Vector constructors, destructors */
+		Vector(){ /*{{{*/
+
+			InitCheckAndSetType();
+		}
+		/*}}}*/
+		Vector(int M,bool fromlocalsize=false){ /*{{{*/
+			
+			InitCheckAndSetType();
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				this->pvector=new PetscVec(M,fromlocalsize);
+				#endif
+			}
+			else this->ivector=new IssmVec<doubletype>(M,fromlocalsize);
+
+		}
+		/*}}}*/
+		Vector(int m,int M){ /*{{{*/
+
+			InitCheckAndSetType();
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+					this->pvector=new PetscVec(m,M);
+				 #endif
+			}
+			else this->ivector=new IssmVec<doubletype>(m,M);
+		}
+		/*}}}*/
+		Vector(doubletype* serial_vec,int M){ /*{{{*/
+			
+			InitCheckAndSetType();
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				this->pvector=new PetscVec(serial_vec,M);
+				#endif
+			}
+			else this->ivector=new IssmVec<doubletype>(serial_vec,M);
+		}
+		/*}}}*/
+		~Vector(){ /*{{{*/
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				delete this->pvector;
+				#endif
+			}
+			else delete this->ivector;
+		}
+		/*}}}*/
+		#ifdef _HAVE_PETSC_
+		Vector(Vec petsc_vector){ /*{{{*/
+
+			this->type=PetscVecType;
+			this->ivector=NULL;
+			this->pvector=new PetscVec(petsc_vector);
+
+		}
+		/*}}}*/
+		#endif
+		void InitCheckAndSetType(void){ /*{{{*/
+
+			char* toolkittype=NULL;
+
+			#ifdef _HAVE_PETSC_
+			pvector=NULL;
+			#endif
+			ivector=NULL;
+
+			/*retrieve toolkittype: */
+			toolkittype=ToolkitOptions::GetToolkitType();
+
+			/*set vector type: */
+			if (strcmp(toolkittype,"petsc")==0){
+				#ifdef _HAVE_PETSC_
+				type=PetscVecType; 
+				#else
+				_error_("cannot create petsc vector without PETSC compiled!");
+				#endif
+			}
+			else if(strcmp(toolkittype,"issm")==0){
+				/*let this choice stand:*/
+				type=IssmVecType;
+			}
+			else _error_("unknow toolkit type ");
+			
+			/*Free ressources: */
+			xDelete<char>(toolkittype);
+		}
+		/*}}}*/
+
+		/*Vector specific routines*/
+		/*FUNCTION Echo{{{*/
+		void Echo(void){_assert_(this);
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				this->pvector->Echo();
+				#endif
+			}
+			else this->ivector->Echo();
+
+		}
+		/*}}}*/
+		/*FUNCTION Assemble{{{*/
+		void Assemble(void){_assert_(this);
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				this->pvector->Assemble();
+				#endif
+			}
+			else this->ivector->Assemble();
+
+		}
+		/*}}}*/
+		/*FUNCTION SetValues{{{*/
+		void SetValues(int ssize, int* list, doubletype* values, InsMode mode){ _assert_(this);
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				this->pvector->SetValues(ssize,list,values,mode);
+				#endif
+			}
+			else this->ivector->SetValues(ssize,list,values,mode);
+
+		}
+		/*}}}*/
+		/*FUNCTION SetValue{{{*/
+		void SetValue(int dof, doubletype value, InsMode mode){_assert_(this);
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				this->pvector->SetValue(dof,value,mode);
+				#endif
+			}
+			else this->ivector->SetValue(dof,value,mode);
+
+		}
+		/*}}}*/
+		/*FUNCTION GetValue{{{*/
+		void GetValue(doubletype* pvalue,int dof){_assert_(this);
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				this->pvector->GetValue(pvalue,dof);
+				#endif
+			}
+			else this->ivector->GetValue(pvalue,dof);
+
+		}
+		/*}}}*/
+		/*FUNCTION GetSize{{{*/
+		void GetSize(int* pM){_assert_(this);
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				this->pvector->GetSize(pM);
+				#endif
+			}
+			else this->ivector->GetSize(pM);
+
+		}
+		/*}}}*/
+		/*FUNCTION IsEmpty{{{*/
+		bool IsEmpty(void){
+			int M;
+			
+			_assert_(this);
+			this->GetSize(&M);
+
+			if(M==0) 
+				return true;
+			else
+				return false;
+		}
+		/*}}}*/
+		/*FUNCTION GetLocalSize{{{*/
+		void GetLocalSize(int* pM){_assert_(this);
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				this->pvector->GetLocalSize(pM);
+				#endif
+			}
+			else this->ivector->GetLocalSize(pM);
+
+		}
+		/*}}}*/
+		/*FUNCTION Duplicate{{{*/
+		Vector<doubletype>* Duplicate(void){_assert_(this);
+
+			Vector<doubletype>* output=NULL;
+				
+			output=new Vector<doubletype>();
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				output->pvector=this->pvector->Duplicate();
+				#endif
+			}
+			else output->ivector=this->ivector->Duplicate();
+
+			return output;
+
+		}
+		/*}}}*/
+		/*FUNCTION Set{{{*/
+		void Set(doubletype value){_assert_(this);
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				this->pvector->Set(value);
+				#endif
+			}
+			else this->ivector->Set(value);
+
+		}
+		/*}}}*/
+		/*FUNCTION AXPY{{{*/
+		void AXPY(Vector* X, doubletype a){_assert_(this);
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				this->pvector->AXPY(X->pvector,a);
+				#endif
+			}
+			else this->ivector->AXPY(X->ivector,a);
+
+		}
+		/*}}}*/
+		/*FUNCTION AYPX{{{*/
+		void AYPX(Vector* X, doubletype a){_assert_(this);
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				this->pvector->AYPX(X->pvector,a);
+				#endif
+			}
+			else this->ivector->AYPX(X->ivector,a);
+		}
+		/*}}}*/
+		/*FUNCTION ToMPISerial{{{*/
+		doubletype* ToMPISerial(void){
+
+			doubletype* vec_serial=NULL;
+			
+			_assert_(this);
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				vec_serial=this->pvector->ToMPISerial();
+				#endif
+			}
+			else vec_serial=this->ivector->ToMPISerial();
+
+			return vec_serial;
+
+		}
+		/*}}}*/
+		/*FUNCTION Copy{{{*/
+		void Copy(Vector* to){_assert_(this);
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				this->pvector->Copy(to->pvector);
+				#endif
+			}
+			else this->ivector->Copy(to->ivector);
+		}
+		/*}}}*/
+		/*FUNCTION Norm{{{*/
+		doubletype Norm(NormMode norm_type){_assert_(this);
+
+			doubletype norm=0;
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				norm=this->pvector->Norm(norm_type);
+				#endif
+			}
+			else norm=this->ivector->Norm(norm_type);
+			return norm;
+		}
+		/*}}}*/
+		/*FUNCTION Scale{{{*/
+		void Scale(doubletype scale_factor){_assert_(this);
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				this->pvector->Scale(scale_factor);
+				#endif
+			}
+			else this->ivector->Scale(scale_factor);
+		}
+		/*}}}*/
+		/*FUNCTION Dot{{{*/
+		doubletype Dot(Vector* vector){_assert_(this);
+
+			doubletype dot;
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				dot=this->pvector->Dot(vector->pvector);
+				#endif
+			}
+			else dot=this->ivector->Dot(vector->ivector);
+			return dot;
+		}
+		/*}}}*/
+		/*FUNCTION PointwiseDivide{{{*/
+		void PointwiseDivide(Vector* x,Vector* y){_assert_(this);
+
+			if(type==PetscVecType){
+				#ifdef _HAVE_PETSC_
+				this->pvector->PointwiseDivide(x->pvector,y->pvector);
+				#endif
+			}
+			else this->ivector->PointwiseDivide(x->ivector,y->ivector);
+		}
+		/*}}}*/
+};
+#endif //#ifndef _VECTOR_H_
Index: ../trunk-jpl/src/c/classes/toolkits/Matrix.h
===================================================================
--- ../trunk-jpl/src/c/classes/toolkits/Matrix.h	(revision 0)
+++ ../trunk-jpl/src/c/classes/toolkits/Matrix.h	(revision 14792)
@@ -0,0 +1,329 @@
+/*!\file:  Matrix.h
+ * \brief wrapper to matrix objects. The goal is to control which API (PETSc,Scalpack, Plapack?) 
+ * implements our underlying matrix format.
+ */ 
+
+#ifndef _MATRIX_H_
+#define _MATRIX_H_
+
+/*Headers:*/
+/*{{{*/
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+#include <cstring>
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+/*}}}*/
+
+enum matrixtype { PetscMatType, IssmMatType };
+
+template <class doubletype> class Vector;
+
+template <class doubletype> 
+class Matrix{
+
+	public:
+
+		int       type;
+		#ifdef _HAVE_PETSC_
+		PetscMat              *pmatrix;
+		#endif
+		IssmMat<doubletype>   *imatrix;
+
+		/*Matrix constructors, destructors*/
+		/*FUNCTION Matrix(){{{*/
+		Matrix(){
+			InitCheckAndSetType();
+		}
+		/*}}}*/
+		/*FUNCTION Matrix(int M,int N){{{*/
+		Matrix(int M,int N){
+			
+			InitCheckAndSetType();
+
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				this->pmatrix=new PetscMat(M,N);
+				#endif
+			}
+			else{
+				this->imatrix=new IssmMat<doubletype>(M,N);
+			}
+
+		}
+		/*}}}*/
+		/*FUNCTION Matrix(int m,int n,int M,int N,int* d_nnz,int* o_nnz){{{*/
+		Matrix(int m,int n,int M,int N,int* d_nnz,int* o_nnz){
+
+			InitCheckAndSetType();
+
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				this->pmatrix=new PetscMat(m,n,M,N,d_nnz,o_nnz);
+				#endif
+			}
+			else{
+				this->imatrix=new IssmMat<doubletype>(m,n,M,N,d_nnz,o_nnz);
+			}
+
+		}
+		/*}}}*/
+		/*FUNCTION Matrix(int M,int N,IssmDouble sparsity){{{*/
+		Matrix(int M,int N,double sparsity){
+			
+			InitCheckAndSetType();
+
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				this->pmatrix=new PetscMat(M,N,sparsity);
+				#endif
+			}
+			else{
+				this->imatrix=new IssmMat<doubletype>(M,N,sparsity);
+			}
+		}
+		/*}}}*/
+		/*FUNCTION Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity){{{*/
+		Matrix(IssmPDouble* serial_mat,int M,int N,IssmPDouble sparsity){
+			
+			InitCheckAndSetType();
+
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				this->pmatrix=new PetscMat(serial_mat,M,N,sparsity);
+				#endif
+			}
+			else{
+				this->imatrix=new IssmMat<doubletype>(serial_mat,M,N,sparsity);
+			}
+
+		}
+		/*}}}*/
+		/*FUNCTION Matrix(int M,int N,int connectivity,int numberofdofspernode){{{*/
+		Matrix(int M,int N,int connectivity,int numberofdofspernode){
+			
+			InitCheckAndSetType();
+
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				this->pmatrix=new PetscMat(M,N,connectivity,numberofdofspernode);
+				#endif
+			}
+			else{
+				this->imatrix=new IssmMat<doubletype>(M,N,connectivity,numberofdofspernode);
+			}
+
+		}
+		/*}}}*/
+		/*FUNCTION ~Matrix(){{{*/
+		~Matrix(){
+
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				delete this->pmatrix;
+				#endif
+			}
+			else delete this->imatrix;
+
+		}
+		/*}}}*/
+		/*FUNCTION InitCheckAndSetType(){{{*/
+		void InitCheckAndSetType(void){
+
+			char* toolkittype=NULL;
+
+			#ifdef _HAVE_PETSC_
+			pmatrix=NULL;
+			#endif
+			imatrix=NULL;
+
+			/*retrieve toolkittype: */
+			toolkittype=ToolkitOptions::GetToolkitType();
+
+			/*set matrix type: */
+			if (strcmp(toolkittype,"petsc")==0){
+				#ifdef _HAVE_PETSC_
+				type=PetscMatType; 
+				#else
+				_error_("cannot create petsc matrix without PETSC compiled!");
+				#endif
+			}
+			else if(strcmp(toolkittype,"issm")==0){
+				/*let this choice stand:*/
+				type=IssmMatType;
+			}
+			else _error_("unknow toolkit type ");
+			
+			/*Free ressources: */
+			xDelete<char>(toolkittype);
+		}
+		/*}}}*/
+
+		/*Matrix specific routines:*/
+		/*FUNCTION Echo{{{*/
+		void Echo(void){
+			_assert_(this);
+
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				this->pmatrix->Echo();
+				#endif
+			}
+			else{
+				this->imatrix->Echo();
+			}
+
+		}
+		/*}}}*/
+		/*FUNCTION AllocationInfo{{{*/
+		void AllocationInfo(void){
+			_assert_(this);
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				this->pmatrix->AllocationInfo();
+				#endif
+			}
+			else{
+				this->imatrix->AllocationInfo();
+			}
+		}/*}}}*/
+		/*FUNCTION Assemble{{{*/
+		void Assemble(void){
+
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				this->pmatrix->Assemble();
+				#endif
+			}
+			else{
+				this->imatrix->Assemble();
+			}
+		}
+		/*}}}*/
+		/*FUNCTION Norm{{{*/
+		IssmDouble Norm(NormMode norm_type){
+
+			IssmDouble norm=0;
+
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				norm=this->pmatrix->Norm(norm_type);
+				#endif
+			}
+			else{
+				norm=this->imatrix->Norm(norm_type);
+			}
+
+			return norm;
+		}
+		/*}}}*/
+		/*FUNCTION GetSize{{{*/
+		void GetSize(int* pM,int* pN){
+
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				this->pmatrix->GetSize(pM,pN);
+				#endif
+			}
+			else{
+				this->imatrix->GetSize(pM,pN);
+			}
+
+		}
+		/*}}}*/
+		/*FUNCTION GetLocalSize{{{*/
+		void GetLocalSize(int* pM,int* pN){
+
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				this->pmatrix->GetLocalSize(pM,pN);
+				#endif
+			}
+			else{
+				this->imatrix->GetLocalSize(pM,pN);
+			}
+
+		}
+		/*}}}*/
+		/*FUNCTION MatMult{{{*/
+		void MatMult(Vector<doubletype>* X,Vector<doubletype>* AX){
+
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				this->pmatrix->MatMult(X->pvector,AX->pvector);
+				#endif
+			}
+			else{
+				this->imatrix->MatMult(X->ivector,AX->ivector);
+			}
+
+		}
+		/*}}}*/
+		/*FUNCTION Duplicate{{{*/
+		Matrix<doubletype>* Duplicate(void){
+
+			Matrix<doubletype>* output=new Matrix<doubletype>();
+
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				output->pmatrix=this->pmatrix->Duplicate();
+				#endif
+			}
+			else{
+				output->imatrix=this->imatrix->Duplicate();
+			}
+
+			return output;
+		}
+		/*}}}*/
+		/*FUNCTION ToSerial{{{*/
+		doubletype* ToSerial(void){
+
+			doubletype* output=NULL;
+
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				output=this->pmatrix->ToSerial();
+				#endif
+			}
+			else{
+				output=this->imatrix->ToSerial();
+			}
+
+			return output;
+		}
+		/*}}}*/
+		/*FUNCTION SetValues{{{*/
+		void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){
+
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				this->pmatrix->SetValues(m,idxm,n,idxn,values,mode);
+				#endif
+			}
+			else{
+				this->imatrix->SetValues(m,idxm,n,idxn,values,mode);
+			}
+		}
+		/*}}}*/
+		/*FUNCTION Convert{{{*/
+		void Convert(MatrixType newtype){
+
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				this->pmatrix->Convert(newtype);
+				#endif
+			}
+			else{
+				this->imatrix->Convert(newtype);
+			}
+
+		}
+		/*}}}*/
+
+};
+
+#endif //#ifndef _MATRIX_H_