Index: /issm/trunk-jpl/src/c/toolkits/mumps/MpiDenseMumpsSolve.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/mumps/MpiDenseMumpsSolve.cpp	(revision 16131)
+++ /issm/trunk-jpl/src/c/toolkits/mumps/MpiDenseMumpsSolve.cpp	(revision 16132)
@@ -262,5 +262,4 @@
     packedDimsSparseArr[3+local_nnz+i]=jcn_loc[i];
   }
-  ensureContiguousLocations(local_nnz+n);
   IssmDouble *pack_A_rhs=xNew<IssmDouble>(local_nnz+n);
   for (int i=0;i<local_nnz;++i) { 
@@ -272,5 +271,4 @@
   IssmPDouble *passivePack_A_rhs=xNew<IssmPDouble>(local_nnz+n);
   IssmPDouble *passiveSol=xNew<IssmPDouble>(n);
-  ensureContiguousLocations(n);
   IssmDouble *sol=xNew<IssmDouble>(n);
   call_ext_fct(dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p,
@@ -302,5 +300,5 @@
 	  local_irn[i]=iArr[3+i];
 	  local_jcn[i]=iArr[3+local_nnz+i];
-    a_loc[i]=dp_x[i];
+	  a_loc[i]=dp_x[i];
   }
   IssmPDouble *rhs_sol=xNew<IssmPDouble>(n);
@@ -326,10 +324,14 @@
   MumpsFinalize(theMumpsStruc);
   // update the adjoint of the rhs:
-  for (int i=0;i<m;++i) { 
+  for (int i=0;i<n;++i) {
     dp_Z[local_nnz+i]+=rhs_sol[i];
   }
-  // update the adjoint of the matrix: 
+  // bcast dp_y (the solution of the forward system)
+  ISSM_MPI_Bcast(dp_y,n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+  // bcast the adjoint of the right-hand-side, i.e. this solution
+  ISSM_MPI_Bcast(rhs_sol,n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+  // update the adjoint of the matrix with the outer product of right-hand-side adjoint and the original solution
   for (int i=0;i<local_nnz;++i) {
-    dp_Z[i]+=-dp_U[local_irn[i]-1]*dp_y[local_jcn[i]-1];
+    dp_Z[i]-=rhs_sol[iArr[3+i]-1]*dp_y[iArr[3+local_nnz+i]-1];
   }
   xDelete(rhs_sol);