Ice Sheet System Model  4.18
Code documentation
MumpsSolve.cpp
Go to the documentation of this file.
1 
5 /*Header files: */
6 #ifdef HAVE_CONFIG_H
7  #include <config.h>
8 #else
9 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
10 #endif
11 
12 #include "../../shared/Numerics/types.h"
13 #include "../../shared/MemOps/MemOps.h"
14 #include "../../shared/Exceptions/exceptions.h"
15 #include "../../shared/io/Comm/IssmComm.h"
16 #include "../../classes/Params/GenericParam.h"
17 #include "../../classes/Params/Parameters.h"
18 #include "../mpi/issmmpi.h"
19 #include "../adolc/adolcincludes.h"
20 #include "../codipack/codipackincludes.h"
21 #include "../issm/SparseRow.h"
22 #include "./mumpsincludes.h"
23 
24 /*Mumps header files: */
25 #include <dmumps_c.h>
26 
27 void MumpsInit(DMUMPS_STRUC_C &theMumpsStruc){
28  theMumpsStruc.n = 0;
29  theMumpsStruc.nz = 0;
30  theMumpsStruc.a = NULL;
31  theMumpsStruc.jcn = NULL;
32  theMumpsStruc.irn = NULL;
33  theMumpsStruc.par = 1;
34  theMumpsStruc.sym = 0;
35  #ifdef _HAVE_MPI_
36  theMumpsStruc.comm_fortran = MPI_Comm_c2f(IssmComm::GetComm());
37  #endif
38  theMumpsStruc.job = -1;
39  dmumps_c(&theMumpsStruc);
40 }
41 
42 // must be preceded by a call to MumpsInit
43 void MumpsSettings(DMUMPS_STRUC_C &theMumpsStruc) {
44  /*Control statements:{{{ */
45  theMumpsStruc.icntl[1-1] = 0; //error verbose: default 6, 0 or negative -> suppressed
46  theMumpsStruc.icntl[2-1] = 0; //std verbose: default 1, 0 or negative -> suppressed
47  theMumpsStruc.icntl[3-1] = 0; //global information verbose: default 6, 0 or negative -> suppressed
48  theMumpsStruc.icntl[4-1] = 0; //verbose everything: default is 4
49  theMumpsStruc.icntl[5-1] = 0;
50  theMumpsStruc.icntl[18-1] = 3;
51 
52  theMumpsStruc.icntl[20-1] = 0;
53  theMumpsStruc.icntl[21-1] = 0;
54  theMumpsStruc.icntl[30-1] = 0;
55  /*}}}*/
56 }
57 
58 // must be preceded by a call to MumpsInit
59 void MumpsAnalyze(DMUMPS_STRUC_C &theMumpsStruc) {
60  theMumpsStruc.job = 1;
61  dmumps_c(&theMumpsStruc);
62 }
63 
64 // must be preceded by a call to MumpsAnalyze
65 void MumpsFactorize(DMUMPS_STRUC_C &theMumpsStruc) {
66  theMumpsStruc.job = 2;
67  dmumps_c(&theMumpsStruc);
68 }
69 
70 // must be preceded by a call to MumpsFactorize
71 void MumpsBacksubstitute(DMUMPS_STRUC_C &theMumpsStruc) {
72  theMumpsStruc.job = 3;
73  dmumps_c(&theMumpsStruc);
74 }
75 
76 // must be preceded at least by a call to MumpsInit
77 void MumpsFinalize(DMUMPS_STRUC_C &theMumpsStruc) {
78  theMumpsStruc.job = -2;
79  dmumps_c(&theMumpsStruc);
80 }
81 
82 void MumpsSolve(int n,int nnz,int local_nnz,int* irn_loc,int* jcn_loc, IssmPDouble *a_loc, IssmPDouble *rhs, Parameters* parameters=0 /*unused here*/){
83  /*Initialize mumps*/
84  DMUMPS_STRUC_C theMumpsStruc;
85  MumpsInit(theMumpsStruc);
86  MumpsSettings(theMumpsStruc);
87 
88  /*now setup the rest of theMumpsStruc */
89  theMumpsStruc.n=n;
90  theMumpsStruc.nz=nnz;
91  theMumpsStruc.nz_loc=local_nnz;
92  theMumpsStruc.irn_loc=irn_loc;
93  theMumpsStruc.jcn_loc=jcn_loc;
94  theMumpsStruc.a_loc=a_loc;
95  theMumpsStruc.rhs=rhs;
96  theMumpsStruc.nrhs=1;
97  theMumpsStruc.lrhs=1;
98 
99  /*Solve system*/
100  MumpsAnalyze(theMumpsStruc);
101  MumpsFactorize(theMumpsStruc);
102  MumpsBacksubstitute(theMumpsStruc);
103  MumpsFinalize(theMumpsStruc);
104 }
105 
106 #ifdef _HAVE_AD_
107 // prototype for active variant
108 void MumpsSolve(int n,
109  int nnz,
110  int local_nnz,
111  int* irn_loc,
112  int* jcn_loc,
113  IssmDouble *a_loc,
114  IssmDouble *rhs,
115  Parameters* parameters);
116 #endif
117 
118 #ifdef _HAVE_MPI_
119 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, Parameters* parameters){ /*{{{*/
120 
121  /*Variables*/
122  ISSM_MPI_Comm comm;
123  int num_procs;
124  int i,j;
125  int nnz,local_nnz;
126  int *irn_loc = NULL;
127  int *jcn_loc = NULL;
128  IssmDouble *a_loc = NULL;
129  int count;
130  int lower_row;
131  int upper_row;
132  IssmDouble *rhs = NULL;
133  int *recvcounts = NULL;
134  int *displs = NULL;
135 
136  /*Communicator info */
137  num_procs = IssmComm::GetSize();
138  comm = IssmComm::GetComm();
139 
140  /*First, some checks*/
141  if (Kff_M!=Kff_N)_error_("stiffness matrix Kff should be square");
142  if (uf_M!=Kff_M || uf_M!=pf_M)_error_("solution vector should be the same size as stiffness matrix Kff and load vector pf");
143  if (uf_m!=Kff_m || uf_m!=pf_m)_error_("solution vector should be locally the same size as stiffness matrix Kff and load vector pf");
144 
145  /*Initialize matrix */
146  /*figure out number of non-zero entries: */
147  local_nnz=0;
148  for(i=0;i<Kff_m;i++){
149  for(j=0;j<Kff_N;j++){
150  if (Kff[i*Kff_N+j]!=0)local_nnz++;
151  }
152  }
153 
154  ISSM_MPI_Reduce(&local_nnz,&nnz,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,comm);
155  ISSM_MPI_Bcast(&nnz,1,ISSM_MPI_INT,0,comm);
156 
157  /*Allocate: */
158  if(local_nnz){
159  irn_loc=xNew<int>(local_nnz);
160  jcn_loc=xNew<int>(local_nnz);
161  a_loc=xNew<IssmDouble>(local_nnz);
162  }
163 
164  /*Populate the triplets: */
165  GetOwnershipBoundariesFromRange(&lower_row,&upper_row,Kff_m,comm);
166  count=0;
167  for(i=0;i<Kff_m;i++){
168  for(j=0;j<Kff_N;j++){
169  if (Kff[i*Kff_N+j]!=0){
170  irn_loc[count]=lower_row+i+1; //fortran indexing
171  jcn_loc[count]=j+1; //fortran indexing
172  a_loc[count]=Kff[i*Kff_N+j];
173  count++;
174  }
175  }
176  }
177  /*Deal with right hand side. We need to ISSM_MPI_Gather it onto cpu 0: */
178 // AD performance is sensitive to calls to ensurecontiguous.
179 // Providing "t" will cause ensurecontiguous to be called.
180 #ifdef _HAVE_AD_
181  rhs=xNew<IssmDouble>(pf_M,"t");
182 #else
183  rhs=xNew<IssmDouble>(pf_M);
184 #endif
185 
186  recvcounts=xNew<int>(num_procs);
187  displs=xNew<int>(num_procs);
188 
189  /*recvcounts:*/
190  ISSM_MPI_Allgather(&pf_m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,comm);
191 
192  /*displs: */
193  ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,comm);
194 
195  /*Gather:*/
196  ISSM_MPI_Gatherv(pf, pf_m, ISSM_MPI_DOUBLE, rhs, recvcounts, displs, ISSM_MPI_DOUBLE,0,comm);
197 
198  MumpsSolve(Kff_M,nnz,local_nnz,irn_loc,jcn_loc,a_loc,rhs,parameters);
199 
200  /*Now scatter from cpu 0 to all other cpus*/
201  ISSM_MPI_Scatterv( rhs, recvcounts, displs, ISSM_MPI_DOUBLE, uf, uf_m, ISSM_MPI_DOUBLE, 0, comm);
202 
203  /*Cleanup*/
204  xDelete<int>(irn_loc);
205  xDelete<int>(jcn_loc);
206  xDelete<IssmDouble>(a_loc);
207  xDelete<IssmDouble>(rhs);
208  xDelete<int>(recvcounts);
209  xDelete<int>(displs);
210 } /*}}}*/
211 void MpiSparseMumpsSolve( /*output: */ IssmDouble* uf, int uf_M, int uf_m, /*matrix input: */ SparseRow<IssmDouble>** Kff, int Kff_M, int Kff_N, int Kff_m, /*right hand side vector: */ IssmDouble* pf, int pf_M, int pf_m, Parameters* parameters){ /*{{{*/
212 
213  /*Variables*/
214  ISSM_MPI_Comm comm;
215  int num_procs;
216  int i;
217  int nnz,local_nnz;
218  int *irn_loc = NULL;
219  int *jcn_loc = NULL;
220  IssmDouble *a_loc = NULL;
221  int count;
222  int lower_row;
223  int upper_row;
224  IssmDouble *rhs = NULL;
225  int *recvcounts = NULL;
226  int *displs = NULL;
227 
228  /*Communicator info */
229  num_procs = IssmComm::GetSize();
230  comm = IssmComm::GetComm();
231 
232  /*First, some checks*/
233  if(Kff_M!=Kff_N)_error_("stiffness matrix Kff should be square");
234  if(uf_M!=Kff_M || uf_M!=pf_M)_error_("solution vector should be the same size as stiffness matrix Kff and load vector pf");
235  if(uf_m!=Kff_m || uf_m!=pf_m)_error_("solution vector should be locally the same size as stiffness matrix Kff and load vector pf");
236 
237  /*Initialize matrix */
238  /*figure out number of non-zero entries: */
239  local_nnz=0;
240  for(i=0;i<Kff_m;i++){
241  local_nnz+=Kff[i]->Nnz();
242  }
243 
244  ISSM_MPI_Reduce(&local_nnz,&nnz,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,comm);
245  ISSM_MPI_Bcast(&nnz,1,ISSM_MPI_INT,0,comm);
246 
247  /*Allocate: */
248  if(local_nnz){
249  irn_loc=xNew<int>(local_nnz);
250  jcn_loc=xNew<int>(local_nnz);
251  a_loc=xNew<IssmDouble>(local_nnz);
252  }
253 
254  /*Populate the triplets: */
255  GetOwnershipBoundariesFromRange(&lower_row,&upper_row,Kff_m,comm);
256  count=0;
257  for(i=0;i<Kff_m;i++){
258  Kff[i]->SetIrnJcnA(irn_loc,jcn_loc,a_loc,lower_row+i+1,count);
259  count+=Kff[i]->Nnz();
260  }
261  /*Deal with right hand side. We need to ISSM_MPI_Gather it onto cpu 0: */
262 // AD performance is sensitive to calls to ensurecontiguous.
263 // Providing "t" will cause ensurecontiguous to be called.
264 #ifdef _HAVE_AD_
265  rhs=xNew<IssmDouble>(pf_M,"t");
266 #else
267  rhs=xNew<IssmDouble>(pf_M);
268 #endif
269 
270  recvcounts=xNew<int>(num_procs);
271  displs=xNew<int>(num_procs);
272 
273  /*recvcounts:*/
274  ISSM_MPI_Allgather(&pf_m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,comm);
275 
276  /*displs: */
277  ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,comm);
278 
279  /*Gather:*/
280  ISSM_MPI_Gatherv(pf, pf_m, ISSM_MPI_DOUBLE, rhs, recvcounts, displs, ISSM_MPI_DOUBLE,0,comm);
281 
282  MumpsSolve(Kff_M,nnz,local_nnz,irn_loc,jcn_loc,a_loc,rhs,parameters);
283 
284  /*Now scatter from cpu 0 to all other cpus*/
285  ISSM_MPI_Scatterv( rhs, recvcounts, displs, ISSM_MPI_DOUBLE, uf, uf_m, ISSM_MPI_DOUBLE, 0, comm);
286 
287  /*Cleanup*/
288  xDelete<int>(irn_loc);
289  xDelete<int>(jcn_loc);
290  xDelete<IssmDouble>(a_loc);
291  xDelete<IssmDouble>(rhs);
292  xDelete<int>(recvcounts);
293  xDelete<int>(displs);
294 } /*}}}*/
295 #endif
296 void SeqDenseMumpsSolve(IssmDouble* uf,int uf_M,int uf_m, IssmDouble* Kff,int Kff_M, int Kff_N, int Kff_m, IssmDouble* pf, int pf_M, int pf_m, Parameters* parameters){/*{{{*/
297 
298  /*Variables*/
299  int *irn_loc = NULL;
300  int *jcn_loc = NULL;
301  IssmDouble *a_loc = NULL;
302  IssmDouble *rhs = NULL;
303 
304  /*First, some checks*/
305  if (Kff_M!=Kff_N)_error_("stiffness matrix Kff should be square");
306  if (Kff_M!=Kff_m)_error_("stiffness matrix Kff is not serial");
307  if (uf_M!=Kff_M || uf_M!=pf_M)_error_("solution vector should be the same size as stiffness matrix Kff and load vector pf");
308  if (uf_m!=Kff_m || uf_m!=pf_m)_error_("solution vector should be locally the same size as stiffness matrix Kff and load vector pf");
309 
310  /*Initialize matrix */
311  /*figure out number of non-zero entries: */
312  int nnz = 0;
313  for(int i=0;i<Kff_M;i++){
314  for(int j=0;j<Kff_N;j++){
315  if(Kff[i*Kff_N+j]!=0)nnz++;
316  }
317  }
318 
319  /*Allocate: */
320  if(nnz){
321  irn_loc = xNew<int>(nnz);
322  jcn_loc = xNew<int>(nnz);
323  a_loc = xNew<IssmDouble>(nnz);
324  }
325 
326  /*Populate the triplets: */
327  int count=0;
328  for(int i=0;i<Kff_M;i++){
329  for(int j=0;j<Kff_N;j++){
330  if(Kff[i*Kff_N+j]!=0){
331  irn_loc[count] = i+1; //fortran indexing
332  jcn_loc[count] = j+1; //fortran indexing
333  a_loc[count] = Kff[i*Kff_N+j];
334  count++;
335  }
336  }
337  }
338 
339  /*Deal with right hand side. We need to ISSM_MPI_Gather it onto cpu 0: */
340 // AD performance is sensitive to calls to ensurecontiguous.
341 // Providing "t" will cause ensurecontiguous to be called.
342 #ifdef _HAVE_AD_
343  rhs=xNew<IssmDouble>(pf_M,"t");
344 #else
345  rhs=xNew<IssmDouble>(pf_M);
346 #endif
347  xMemCpy<IssmDouble>(rhs,pf,Kff_M);
348 
349  MumpsSolve(Kff_M,nnz,nnz,irn_loc,jcn_loc,a_loc,rhs,parameters);
350 
351  /*Now scatter from cpu 0 to all other cpus*/
352  xMemCpy<IssmDouble>(uf,rhs,Kff_M);
353 
354  /*Cleanup*/
355  xDelete<int>(irn_loc);
356  xDelete<int>(jcn_loc);
357  xDelete<IssmDouble>(a_loc);
358  xDelete<IssmDouble>(rhs);
359 }/*}}}*/
360 
361 #ifdef _HAVE_ADOLC_
362 int mumpsSolveEDF(int iArrLength, int* iArr, int /* ignored */, IssmPDouble* dp_x, int /* ignored */, IssmPDouble* dp_y){/*{{{*/
363  // unpack parameters
364  int n=iArr[0];
365  int nnz=iArr[1];
366  int local_nnz=iArr[2];
367  int *local_irn=xNew<int>(local_nnz);
368  int *local_jcn=xNew<int>(local_nnz);
369  IssmPDouble *A=xNew<IssmPDouble>(local_nnz);
370  for (int i=0;i<local_nnz;++i) {
371  local_irn[i]=iArr[3+i];
372  local_jcn[i]=iArr[3+local_nnz+i];
373  A[i]=dp_x[i];
374  }
375  IssmPDouble *rhs_sol=xNew<IssmPDouble>(n);
376  for (int i=0;i<n;++i) {
377  rhs_sol[i]=dp_x[local_nnz+i];
378  }
379  MumpsSolve(n,nnz,local_nnz,local_irn,local_jcn,A,rhs_sol);
380  for (int i=0;i<n;++i) {
381  dp_y[i]=rhs_sol[i];
382  }
383  xDelete(rhs_sol);
384  xDelete(A);
385  xDelete(local_jcn);
386  xDelete(local_irn);
387  return 0;
388 }/*}}}*/
389 void MumpsSolve(int n,int nnz,int local_nnz,int* irn_loc,int* jcn_loc,IssmDouble *a_loc,IssmDouble *rhs,Parameters* parameters){/*{{{*/
390  int packedDimsSparseArrLength=1+1+1+local_nnz+local_nnz;
391  int *packedDimsSparseArr=xNew<int>(packedDimsSparseArrLength);
392  packedDimsSparseArr[0]=n;
393  packedDimsSparseArr[1]=nnz;
394  packedDimsSparseArr[2]=local_nnz;
395  for (int i=0;i<local_nnz;++i) {
396  packedDimsSparseArr[3+i]=irn_loc[i];
397  packedDimsSparseArr[3+local_nnz+i]=jcn_loc[i];
398  }
399 // AD performance is sensitive to calls to ensurecontiguous.
400 // Providing "t" will cause ensurecontiguous to be called.
401 #ifdef _HAVE_AD_
402  IssmDouble *pack_A_rhs=xNew<IssmDouble>(local_nnz+n,"t");
403 #else
404  IssmDouble *pack_A_rhs=xNew<IssmDouble>(local_nnz+n);
405 #endif
406 
407  for (int i=0;i<local_nnz;++i) {
408  pack_A_rhs[i]=a_loc[i];
409  }
410  for (int i=0;i<n;++i) {
411  pack_A_rhs[local_nnz+i]=rhs[i];
412  }
413 // AD performance is sensitive to calls to ensurecontiguous.
414 // Providing "t" will cause ensurecontiguous to be called.
415 #ifdef _HAVE_AD_
416  IssmDouble *sol=xNew<IssmDouble>(n,"t");
417 #else
418  IssmDouble *sol=xNew<IssmDouble>(n);
419 #endif
420 
421  call_ext_fct(xDynamicCast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p,
422  packedDimsSparseArrLength, packedDimsSparseArr,
423  local_nnz+n, pack_A_rhs,
424  n,sol);
425  for (int i=0;i<n;++i) {
426  rhs[i]=sol[i];
427  }
428  xDelete(sol);
429  xDelete(pack_A_rhs);
430  xDelete(packedDimsSparseArr);
431 }/*}}}*/
432 int fos_reverse_mumpsSolveEDF(int iArrLength, int* iArr, int m, IssmPDouble *dp_U, int nPlusNz, IssmPDouble *dp_Z, IssmPDouble *dp_x, IssmPDouble *dp_y) {/*{{{*/
433  // unpack parameters
434  int n=iArr[0];
435  int nnz=iArr[1];
436  int local_nnz=iArr[2];
437  int *local_irn=xNew<int>(local_nnz);
438  int *local_jcn=xNew<int>(local_nnz);
439  IssmPDouble *a_loc=xNew<IssmPDouble>(local_nnz);
440  for (int i=0;i<local_nnz;++i) {
441  local_irn[i]=iArr[3+i];
442  local_jcn[i]=iArr[3+local_nnz+i];
443  a_loc[i]=dp_x[i];
444  }
445  IssmPDouble *rhs_sol=xNew<IssmPDouble>(n);
446  for (int i=0;i<n;++i) {
447  rhs_sol[i]=dp_U[i];
448  }
449  DMUMPS_STRUC_C theMumpsStruc;
450  MumpsInit(theMumpsStruc);
451  MumpsSettings(theMumpsStruc);
452  theMumpsStruc.n=n;
453  theMumpsStruc.nz=nnz;
454  theMumpsStruc.nz_loc=local_nnz;
455  theMumpsStruc.irn_loc=local_irn;
456  theMumpsStruc.jcn_loc=local_jcn;
457  theMumpsStruc.a_loc=a_loc;
458  theMumpsStruc.rhs=rhs_sol;
459  theMumpsStruc.nrhs=1;
460  theMumpsStruc.lrhs=1;
461  theMumpsStruc.icntl[9-1] = 0; //solve for the transpose
462  MumpsAnalyze(theMumpsStruc);
463  MumpsFactorize(theMumpsStruc);
464  MumpsBacksubstitute(theMumpsStruc);
465  MumpsFinalize(theMumpsStruc);
466  // update the adjoint of the rhs:
467  for (int i=0;i<n;++i) {
468  dp_Z[local_nnz+i]+=rhs_sol[i];
469  }
470  // bcast dp_y (the solution of the forward system)
471  ISSM_MPI_Bcast(dp_y,n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
472  // bcast the adjoint of the right-hand-side, i.e. this solution
473  ISSM_MPI_Bcast(rhs_sol,n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
474  // update the adjoint of the matrix with the outer product of right-hand-side adjoint and the original solution
475  for (int i=0;i<local_nnz;++i) {
476  dp_Z[i]-=rhs_sol[iArr[3+i]-1]*dp_y[iArr[3+local_nnz+i]-1];
477  }
478  xDelete(rhs_sol);
479  xDelete(a_loc);
480  xDelete(local_jcn);
481  xDelete(local_irn);
482  return 3;
483 }/*}}}*/
484 #endif
485 
486 #ifdef _HAVE_CODIPACK_
487 void MumpsSolve_codi_b(void* tape_in,void* data_in,void* ra) {/*{{{*/
488 
489  /*recast data_in and tape*/
490  codi::DataStore* data = (codi::DataStore*)data_in;
491  //IssmDouble::TapeType& tape = (IssmDouble::TapeType&)tape_in;
492  IssmDouble::TapeType& tape = IssmDouble::getGlobalTape();
493 
494 
495  IssmDouble::Real* valueATrans;
496  IssmDouble::GradientData* indexATrans;
497  int* irnATrans;
498  int* jcnATrans;
499  IssmDouble::GradientData* indexB;
500  IssmDouble::Real* valueX;
501  IssmDouble::GradientData* indexX;
502  int n;
503  int nnz;
504  int local_nnz;
505  Parameters* parameters;
506 
507  data->getData(valueATrans);
508  data->getData(indexATrans);
509  data->getData(irnATrans);
510  data->getData(jcnATrans);
511  data->getData(indexB);
512  data->getData(valueX);
513  data->getData(indexX);
514  data->getData(n);
515  data->getData(nnz);
516  data->getData(local_nnz);
517  data->getData(parameters);
518 
519  // create the adjoint vector for x and reset the adjoint values on the tape
520  IssmDouble::GradientValue* adjX = xNew<IssmDouble::GradientValue>(n);
521  getVectorAdjoint(tape, indexX, adjX, n);
522 
523  MumpsSolve(n, nnz, local_nnz, irnATrans, jcnATrans, valueATrans, adjX, parameters);
524  // adjX contains now the solution
525 
526  updateVectorAdjoint(tape, indexB, adjX, n);
527 
528  // bcast dp_y (the solution of the forward system)
529  ISSM_MPI_Bcast(valueX,n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
530  // bcast the adjoint of the right-hand-side, i.e. this solution
531  ISSM_MPI_Bcast(adjX,n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
532 
533  for(int i=0; i<local_nnz; ++i) {
534  // we access the transposed matrix here because we stored the indices in a transposed way
535  // -1 is substracted because jcn and irn are stored with fortran indexing
536  updateAdjoint(tape, indexATrans[i], -adjX[jcnATrans[i]-1]*valueX[irnATrans[i]-1]);
537  }
538 
539  xDelete(adjX);
540 }
541 /*}}}*/
542 void MumpsSolve_codi_delete(void* tape_in,void* data_in) {/*{{{*/
543 
544  /*recast data_in*/
545  codi::DataStore* data = (codi::DataStore*)data_in;
546 
547  IssmDouble::Real* valueATrans;
548  IssmDouble::GradientData* indexATrans;
549  int* irnATrans;
550  int* jcnATrans;
551  IssmDouble::GradientData* indexB;
552  IssmDouble::Real* valueX;
553  IssmDouble::GradientData* indexX;
554  int n;
555  int nnz;
556  int local_nnz;
557  Parameters* parameters;
558 
559  data->getData(valueATrans);
560  data->getData(indexATrans);
561  data->getData(irnATrans);
562  data->getData(jcnATrans);
563  data->getData(indexB);
564  data->getData(valueX);
565  data->getData(indexX);
566  data->getData(n);
567  data->getData(nnz);
568  data->getData(local_nnz);
569  data->getData(parameters);
570 
571  xDelete(valueATrans);
572  xDelete(indexATrans);
573  xDelete(irnATrans);
574  xDelete(jcnATrans);
575  xDelete(indexB);
576  xDelete(valueX);
577  xDelete(indexX);
578 }
579 /*}}}*/
580 void MumpsSolve(int n,int nnz,int local_nnz,int* irn_loc,int* jcn_loc,IssmDouble *a_loc,IssmDouble *rhs,Parameters* parameters){/*{{{*/
581  IssmDouble::TapeType& tape = IssmDouble::getGlobalTape();
582  codi::DataStore* dataHandler = NULL;
583 
584  if(tape.isActive()) {
585  dataHandler = new codi::DataStore();
586 
587  // create the index and double vector for the matrix
588  IssmDouble::Real* valueATrans = xNew<IssmDouble::Real>(local_nnz);
589  IssmDouble::GradientData* indexATrans = xNew<IssmDouble::GradientData>(local_nnz);
590  int* irnATrans = xNew<int>(local_nnz);
591  int* jcnATrans = xNew<int>(local_nnz);
592 
593  // read the data for the matrix A in a transposed fashion
594  for (int i=0; i<local_nnz; ++i) {
595  getPrimalAndGradData(a_loc[i], valueATrans[i], indexATrans[i]);
596  irnATrans[i]=jcn_loc[i]; // transposed store
597  jcnATrans[i]=irn_loc[i]; // transposed store
598  }
599 
600  // create the index vector for a (primal values are not needed for a)
601  IssmDouble::GradientData* indexB = xNew<IssmDouble::GradientData>(n);
602  getVectorGradData(rhs, indexB, n);
603 
604  dataHandler->addData(valueATrans);
605  dataHandler->addData(indexATrans);
606  dataHandler->addData(irnATrans);
607  dataHandler->addData(jcnATrans);
608  dataHandler->addData(indexB);
609  }
610 
611  // unpack the primal values from the matrix and the vector
612  IssmDouble::Real* valueA = xNew<IssmDouble::Real>(local_nnz);
613  IssmDouble::Real* valueB = xNew<IssmDouble::Real>(n);
614  // read the data from A and B
615  getVectorPrimal(a_loc, valueA, local_nnz);
616  getVectorPrimal(rhs, valueB, n);
617 
618  MumpsSolve(n, nnz, local_nnz, irn_loc, jcn_loc, valueA, valueB, parameters);
619  // valueB contains now the solution
620 
621  // pack the values into rhs
622  setVectorPrimal(rhs, valueB, n);
623 
624  if(tape.isActive()) {
625  // create the index vector X and register x as active variables
626  IssmDouble::GradientData* indexX = xNew<IssmDouble::GradientData>(n);
627  registerVector(rhs, indexX, n);
628 
629  dataHandler->addData(valueB); // contains the values from x
630  dataHandler->addData(indexX);
631 
632  // store other arguments
633  dataHandler->addData(n);
634  dataHandler->addData(nnz);
635  dataHandler->addData(local_nnz);
636  dataHandler->addData(parameters); // we assume here that parameters is still intact when the reverse run is called
637 
638  //tape.pushExternalFunction(&MumpsSolve_codi_b, dataHandler, &MumpsSolve_codi_delete);
639  tape.pushExternalFunctionHandle(&MumpsSolve_codi_b,(void*)dataHandler, &MumpsSolve_codi_delete);
640  } else {
641  // if the tape is active valueB is stored in the dataHandler and deleted in the reverse sweep
642  xDelete(valueB);
643  }
644 
645  xDelete(valueA);
646 }
647 /*}}}*/
648 #endif
SeqDenseMumpsSolve
void SeqDenseMumpsSolve(IssmDouble *uf, int uf_M, int uf_m, IssmDouble *Kff, int Kff_M, int Kff_N, int Kff_m, IssmDouble *pf, int pf_M, int pf_m, Parameters *parameters)
Definition: MumpsSolve.cpp:296
xDynamicCast
To xDynamicCast(const From &from)
Definition: recast.h:51
IssmDouble
double IssmDouble
Definition: types.h:37
MumpsSolve
void MumpsSolve(int n, int nnz, int local_nnz, int *irn_loc, int *jcn_loc, IssmPDouble *a_loc, IssmPDouble *rhs, Parameters *parameters=0)
Definition: MumpsSolve.cpp:82
mumpsincludes.h
ISSM_MPI_Gatherv
int ISSM_MPI_Gatherv(void *sendbuf, int sendcnt, ISSM_MPI_Datatype sendtype, void *recvbuf, int *recvcnts, int *displs, ISSM_MPI_Datatype recvtype, int root, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:279
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
ISSM_MPI_SUM
#define ISSM_MPI_SUM
Definition: issmmpi.h:134
IssmComm::GetComm
static ISSM_MPI_Comm GetComm(void)
Definition: IssmComm.cpp:30
MumpsBacksubstitute
void MumpsBacksubstitute(DMUMPS_STRUC_C &theMumpsStruc)
Definition: MumpsSolve.cpp:71
GetOwnershipBoundariesFromRange
void GetOwnershipBoundariesFromRange(int *plower_row, int *pupper_row, int range, ISSM_MPI_Comm comm)
Definition: GetOwnershipBoundariesFromRange.cpp:16
MumpsInit
void MumpsInit(DMUMPS_STRUC_C &theMumpsStruc)
Definition: MumpsSolve.cpp:27
ISSM_MPI_DOUBLE
#define ISSM_MPI_DOUBLE
Definition: issmmpi.h:125
ISSM_MPI_INT
#define ISSM_MPI_INT
Definition: issmmpi.h:127
GenericParam::GetParameterValue
P & GetParameterValue()
Definition: GenericParam.h:62
xDelete
void xDelete(T **&aT_pp)
Definition: MemOps.h:97
SparseRow
Definition: SparseRow.h:13
IssmComm::GetSize
static int GetSize(void)
Definition: IssmComm.cpp:46
SparseRow::Nnz
int Nnz(void)
Definition: SparseRow.h:138
MumpsAnalyze
void MumpsAnalyze(DMUMPS_STRUC_C &theMumpsStruc)
Definition: MumpsSolve.cpp:59
ISSM_MPI_PDOUBLE
#define ISSM_MPI_PDOUBLE
Definition: issmmpi.h:126
ISSM_MPI_Bcast
int ISSM_MPI_Bcast(void *buffer, int count, ISSM_MPI_Datatype datatype, int root, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:162
ISSM_MPI_Allgather
int ISSM_MPI_Allgather(void *sendbuf, int sendcount, ISSM_MPI_Datatype sendtype, void *recvbuf, int recvcount, ISSM_MPI_Datatype recvtype, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:45
SparseRow::SetIrnJcnA
void SetIrnJcnA(int *irn_loc, int *jcn_loc, doubletype *a_loc, int i_index, int count)
Definition: SparseRow.h:143
MpiDenseMumpsSolve
void MpiDenseMumpsSolve(IssmDouble *uf, int uf_M, int uf_m, IssmDouble *Kff, int Kff_M, int Kff_N, int Kff_m, IssmDouble *pf, int pf_M, int pf_m, Parameters *parameters)
Definition: MpiDenseMumpsSolve.cpp:109
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
MumpsSettings
void MumpsSettings(DMUMPS_STRUC_C &theMumpsStruc)
Definition: MumpsSolve.cpp:43
MumpsFinalize
void MumpsFinalize(DMUMPS_STRUC_C &theMumpsStruc)
Definition: MumpsSolve.cpp:77
ISSM_MPI_Comm
int ISSM_MPI_Comm
Definition: issmmpi.h:118
AdolcParamEnum
@ AdolcParamEnum
Definition: EnumDefinitions.h:12
ISSM_MPI_Reduce
int ISSM_MPI_Reduce(void *sendbuf, void *recvbuf, int count, ISSM_MPI_Datatype datatype, ISSM_MPI_Op op, int root, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:373
IssmPDouble
IssmDouble IssmPDouble
Definition: types.h:38
Parameters::FindParamObject
Param * FindParamObject(int enum_type)
Definition: Parameters.cpp:588
MumpsFactorize
void MumpsFactorize(DMUMPS_STRUC_C &theMumpsStruc)
Definition: MumpsSolve.cpp:65
ISSM_MPI_Scatterv
int ISSM_MPI_Scatterv(void *sendbuf, int *sendcnts, int *displs, ISSM_MPI_Datatype sendtype, void *recvbuf, int recvcnt, ISSM_MPI_Datatype recvtype, int root, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:444
GenericParam
Definition: GenericParam.h:26
Generated on Thu Jul 2 2020 08:09:22 for Ice Sheet System Model by   doxygen 1.8.19