Changeset 16132

Timestamp:

09/13/13 14:02:47 (11 years ago)

Author:

utke

Message:

CHG ensure contigous explicit calls are no longer needed here because
they are embedded in the xNew specialization; m should be ==n but I
changed it anyway; in order to update the adjoint of the matrix I need
first broadcast the original solution and the adjoint of the right hand
side because the matrix is assembled and therefore the individual matrix
partitions must be updated on the respective ranks; also to get the
matrix elements right we use the original values for the indices
in the iArr because the local_irn/local_jcn may be changed during the
factorization.

File:

• issm/trunk-jpl/src/c/toolkits/mumps/MpiDenseMumpsSolve.cpp (modified) (4 diffs)

issm/trunk-jpl/src/c/toolkits/mumps/MpiDenseMumpsSolve.cpp

@@ -262 +262 @@
     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 +271 @@
   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 +300 @@
           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 +324 @@
   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);