Index: /issm/trunk-jpl/src/c/toolkits/mumps/MpiDenseMumpsSolve.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/mumps/MpiDenseMumpsSolve.cpp	(revision 16082)
+++ /issm/trunk-jpl/src/c/toolkits/mumps/MpiDenseMumpsSolve.cpp	(revision 16083)
@@ -14,5 +14,8 @@
 #include "../../shared/Exceptions/exceptions.h"
 #include "../../shared/io/Comm/IssmComm.h"
+#include "../../classes/Params/GenericParam.h"
+#include "../../classes/Params/Parameters.h"
 #include "../mpi/issmmpi.h"
+#include "../adolc/adolcincludes.h"
 #include "./mumpsincludes.h"
 
@@ -75,5 +78,5 @@
 		IssmPDouble *a_loc,
 		IssmPDouble *rhs,
-		Parameters* /*unused here*/) {
+		Parameters* parameters=0 /*unused here*/) {
 	/*Initialize mumps: {{{*/
 	DMUMPS_STRUC_C theMumpsStruc;
@@ -214,24 +217,29 @@
 #ifdef _HAVE_ADOLC_
 
-int mumpsSolveEDF(int iArrLength, int* iArr, int nPlusNz /* we can ignore it*/, double* dp_x, int m, double* dp_y) {
+int mumpsSolveEDF(int iArrLength, int* iArr, int /* ignored */, IssmPDouble* dp_x, int /* ignored */, IssmPDouble* dp_y) {
   // unpack parameters
   int n=iArr[0];
-  int nz=iArr[1];
-  int *irn=new int[nz];
-  int *jcn=new int[nz];
-  double *A=new double[nz];
-  for (int i=0;i<nz;++i) { 
-    irn[i]=iArr[2+i];
-    jcn[i]=iArr[2+nz+i];
+  int nnz=iArr[1];
+  int local_nnz=iArr[2];
+  int *local_irn=xNew<int>(local_nnz);
+  int *local_jcn=xNew<int>(local_nnz);
+  IssmPDouble *A=xNew<IssmPDouble>(local_nnz);
+  for (int i=0;i<local_nnz;++i) {
+    local_irn[i]=iArr[3+i];
+    local_jcn[i]=iArr[3+local_nnz+i];
     A[i]=dp_x[i];
   }
-  double *rhs_sol=new double[n];
+  IssmPDouble *rhs_sol=xNew<IssmPDouble>(n);
   for (int i=0;i<n;++i) { 
-    rhs_sol[i]=dp_x[nz+i];
-  }
-  mumpsSolve(n,nz,irn,jcn,A,rhs_sol);
-  for (int i=0;i<m;++i) { 
+    rhs_sol[i]=dp_x[local_nnz+i];
+  }
+  MumpsSolve(n,nnz,local_nnz,local_irn,local_jcn,A,rhs_sol);
+  for (int i=0;i<n;++i) {
     dp_y[i]=rhs_sol[i];
   }
+  xDelete(rhs_sol);
+  xDelete(A);
+  xDelete(local_jcn);
+  xDelete(local_irn);
   return 0;
 }
@@ -255,5 +263,5 @@
   }
   ensureContiguousLocations(local_nnz+n);
-  adouble *pack_A_rhs=xNew<IssmDouble>(local_nnz+n);
+  IssmDouble *pack_A_rhs=xNew<IssmDouble>(local_nnz+n);
   for (int i=0;i<local_nnz;++i) { 
     pack_A_rhs[i]=a_loc[i];
@@ -262,9 +270,9 @@
     pack_A_rhs[local_nnz+i]=rhs[i];
   }
-  double *passivePack_A_rhs=xNew<IssmPDouble>(local_nnz+n);
-  double *passiveSol=xNew<IssmPDouble>(n);
+  IssmPDouble *passivePack_A_rhs=xNew<IssmPDouble>(local_nnz+n);
+  IssmPDouble *passiveSol=xNew<IssmPDouble>(n);
   ensureContiguousLocations(n);
-  adouble *sol=xNew<IssmDouble>(n);
-  call_ext_fct(ourEDF_p,
+  IssmDouble *sol=xNew<IssmDouble>(n);
+  call_ext_fct(dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p,
 	       packedDimsSparseArrLength, packedDimsSparseArr,
 	       local_nnz+n, passivePack_A_rhs, pack_A_rhs, 
@@ -280,131 +288,55 @@
 }
 
-int fos_forward_mumpsSolveEDF(int iArrLength, int* iArr, int nPlusNz /* we can ignore it*/, 
-			      double *dp_x, double *dp_X, int m, double *dp_y, double *dp_Y) {
+int fos_reverse_mumpsSolveEDF(int iArrLength, int* iArr, 
+			      int m, IssmPDouble *dp_U,
+			      int nPlusNz, IssmPDouble *dp_Z,
+			      IssmPDouble *dp_x, IssmPDouble *dp_y) {
   // unpack parameters
   int n=iArr[0];
-  int nz=iArr[1];
-  int *irn=new int[nz];
-  int *jcn=new int[nz];
-  double *A=new double[nz];
-  for (int i=0;i<nz;++i) { 
-    irn[i]=iArr[2+i];
-    jcn[i]=iArr[2+nz+i];
-    A[i]=dp_x[i];
-  }
-  double *rhs_sol=new double[n];
-  for (int i=0;i<n;++i) { 
-    rhs_sol[i]=dp_x[nz+i];
-  }
-  DMUMPS_STRUC_C id;
-  id.par = 1; // one processor=sequential code
-  id.sym = 0; // asymmetric
-  id.job = JOB_INIT;
-  dmumps_c(&id);
-
-  id.icntl[1-1] = 6; //error verbose
-  id.icntl[2-1] = 0; //std verbose
-  id.icntl[3-1] = 0; // 
-  id.icntl[4-1] = 0; // 
-  id.icntl[5-1] = 0; // matrix is assembled
-  id.icntl[18-1] = 0; // centralized
-  id.icntl[20-1] = 0; // rhs is dense and centralized
-  id.icntl[21-1] = 0; // solution is centralized
-  id.n=n;
-  id.nz=nz;
-  id.irn=irn;
-  id.jcn=jcn;
-  id.a=A;
-  id.job = JOB_ANALYSIS;
-  dmumps_c(&id);
-  id.job = JOB_FACTORIZATION; 
-  dmumps_c (&id);
-  // solve the orifginal system
-  id.rhs=rhs_sol;
-  id.nrhs=1;
-  id.lrhs=1;
-  id.job = JOB_BACKSUBST; 
-  dmumps_c (&id);
-  for (int i=0;i<m;++i) { 
-    dp_y[i]=rhs_sol[i];
-  }
-  // solve for the derivative
-  for (int i=0;i<n;++i) { 
-    rhs_sol[i]=dp_X[nz+i]; 
-  }
-  for (int i=0;i<nz;++i) { 
-    rhs_sol[irn[i]-1]-=dp_X[i]*dp_y[jcn[i]-1];
-  }
-  dmumps_c (&id);
-  for (int i=0;i<m;++i) { 
-    dp_Y[i]=rhs_sol[i];
-  }
-  id.job = JOB_END; 
-  dmumps_c (&id);
-  return 3;
-}
-
-int fos_reverse_mumpsSolveEDF(int iArrLength, int* iArr, 
-			      int m, double *dp_U, 
-			      int nPlusNz, double *dp_Z, 
-			      double *dp_x, double *dp_y) {
-  // unpack parameters
-  int n=iArr[0];
-  int nz=iArr[1];
-  int *irn=new int[nz];
-  int *jcn=new int[nz];
-  double *A=new double[nz];
-  for (int i=0;i<nz;++i) { 
-    irn[i]=iArr[2+i];
-    jcn[i]=iArr[2+nz+i];
-    A[i]=dp_x[i];
-  }
-  DMUMPS_STRUC_C id;
-  id.par = 1; // one processor=sequential code
-  id.sym = 0; // asymmetric
-  id.job = JOB_INIT;
-  dmumps_c(&id);
-
-  id.icntl[1-1] = 6; //error verbose
-  id.icntl[2-1] = 0; //std verbose
-  id.icntl[3-1] = 0; // 
-  id.icntl[4-1] = 0; // 
-  id.icntl[5-1] = 0; // matrix is assembled
-  id.icntl[9-1] = 0; //solve for the transpose
-  id.icntl[18-1] = 0; // centralized
-  id.icntl[20-1] = 0; // rhs is dense and centralized
-  id.icntl[21-1] = 0; // solution is centralized
-  id.n=n;
-  id.nz=nz;
-  id.irn=irn;
-  id.jcn=jcn;
-  id.a=A;
-  id.job = JOB_ANALYSIS;
-  dmumps_c(&id);
-  id.job = JOB_FACTORIZATION; 
-  dmumps_c (&id);
-  double *rhs_sol=new double[n];
+  int nnz=iArr[1];
+  int local_nnz=iArr[2];
+  int *local_irn=xNew<int>(local_nnz);
+  int *local_jcn=xNew<int>(local_nnz);
+  IssmPDouble *a_loc=xNew<IssmPDouble>(local_nnz);
+  for (int i=0;i<local_nnz;++i) {
+	  local_irn[i]=iArr[3+i];
+	  local_jcn[i]=iArr[3+local_nnz+i];
+    a_loc[i]=dp_x[i];
+  }
+  IssmPDouble *rhs_sol=xNew<IssmPDouble>(n);
   for (int i=0;i<n;++i) { 
     rhs_sol[i]=dp_U[i];
   }
-  id.rhs=rhs_sol;
-  id.nrhs=1;
-  id.lrhs=1;
-  id.job = JOB_BACKSUBST; 
-  dmumps_c (&id);
-  // update the adhoint of the rhs: 
+  DMUMPS_STRUC_C theMumpsStruc;
+  MumpsInit(theMumpsStruc);
+  MumpsSettings(theMumpsStruc);
+  theMumpsStruc.n=n;
+  theMumpsStruc.nz=nnz;
+  theMumpsStruc.nz_loc=local_nnz;
+  theMumpsStruc.irn_loc=local_irn;
+  theMumpsStruc.jcn_loc=local_jcn;
+  theMumpsStruc.a_loc=a_loc;
+  theMumpsStruc.rhs=rhs_sol;
+  theMumpsStruc.nrhs=1;
+  theMumpsStruc.lrhs=1;
+  theMumpsStruc.icntl[9-1] = 0; //solve for the transpose
+  MumpsAnalyze(theMumpsStruc);
+  MumpsFactorize(theMumpsStruc);
+  MumpsBacksubstitute(theMumpsStruc);
+  MumpsFinalize(theMumpsStruc);
+  // update the adjoint of the rhs:
   for (int i=0;i<m;++i) { 
-    dp_Z[nz+i]+=rhs_sol[i];
+    dp_Z[local_nnz+i]+=rhs_sol[i];
   }
   // update the adjoint of the matrix: 
-  for (int i=0;i<nz;++i) { 
-    dp_Z[i]+=-dp_U[irn[i]-1]*dp_y[jcn[i]-1];
-  }
+  for (int i=0;i<local_nnz;++i) {
+    dp_Z[i]+=-dp_U[local_irn[i]-1]*dp_y[local_jcn[i]-1];
+  }
+  xDelete(rhs_sol);
+  xDelete(a_loc);
+  xDelete(local_jcn);
+  xDelete(local_irn);
   return 3;
 }
 
-void MpiDenseMumpsSolve( /*output: */ IssmDouble* uf, int uf_M, int uf_m, /*matrix input: */ IssmDouble* Kff, int Kff_M, int Kff_N, int Kff_m, /*right hand side vector: */ IssmDouble* pf, int pf_M, int pf_m){ /*{{{*/
-	_error_("not supported yet!");
-} /*}}}*/
-
 #endif