Changeset 16066


Ignore:
Timestamp:
09/04/13 12:13:15 (12 years ago)
Author:
utke
Message:

CHG refactored mumps logic and setup skeleton for mumps solver adjoint

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r15839 r16066  
    1818/*Mumps header files: */
    1919#include "dmumps_c.h"
    20 
    21 #define JOB_INIT -1
    22 #define JOB_SOLVE 6
    23 #define JOB_END  -2
    2420/*}}}*/
    2521
    26 void MpiDenseMumpsSolve( /*output: */ IssmPDouble* uf, int uf_M, int uf_m, /*matrix input: */ IssmPDouble* Kff, int Kff_M, int Kff_N, int Kff_m, /*right hand side vector: */ IssmPDouble* pf, int pf_M, int pf_m){ /*{{{*/
     22void MumpsSettings(DMUMPS_STRUC_C &theMumpsStruc,ISSM_MPI_Comm &comm) {
     23        theMumpsStruc.par          = 1; 
     24        theMumpsStruc.sym          = 0;
     25        theMumpsStruc.comm_fortran = MPI_Comm_c2f(comm);
     26        /*Control statements:{{{ */
     27        theMumpsStruc.icntl[1-1] = 6; //error verbose
     28        theMumpsStruc.icntl[2-1] = 1; //std verbose
     29        theMumpsStruc.icntl[4-1] = 4; //verbose everything
     30        theMumpsStruc.icntl[5-1] = 0;
     31        theMumpsStruc.icntl[18-1] = 3;
     32
     33        theMumpsStruc.icntl[20-1] = 0;
     34        theMumpsStruc.icntl[21-1] = 0;
     35        theMumpsStruc.icntl[30-1] = 0;
     36        /*}}}*/
     37}
     38
     39// must be preceded by a call to MumpsSettings
     40void MumpsInit(DMUMPS_STRUC_C &theMumpsStruc) {
     41        theMumpsStruc.job          = -1;
     42        dmumps_c(&theMumpsStruc);
     43}
     44
     45// must be preceded by a call to MumpsInit
     46void MumpsAnalyze(DMUMPS_STRUC_C &theMumpsStruc) {
     47        theMumpsStruc.job          = 1;
     48        dmumps_c(&theMumpsStruc);
     49}
     50
     51// must be preceded by a call to MumpsAnalyze
     52void MumpsFactorize(DMUMPS_STRUC_C &theMumpsStruc) {
     53        theMumpsStruc.job          = 2;
     54        dmumps_c(&theMumpsStruc);
     55}
     56
     57// must be preceded by a call to MumpsFactorize
     58void MumpsBacksubstitute(DMUMPS_STRUC_C &theMumpsStruc) {
     59        theMumpsStruc.job          = 3;
     60        dmumps_c(&theMumpsStruc);
     61}
     62
     63// must be preceded at least  by a call to MumpsInit
     64void MumpsFinalize(DMUMPS_STRUC_C &theMumpsStruc) {
     65        theMumpsStruc.job          = -2;
     66        dmumps_c(&theMumpsStruc);
     67}
     68
     69void MumpsSolve(int n,
     70                int nnz,
     71                int local_nnz,
     72                int* irn_loc,
     73                int* jcn_loc,
     74                IssmPDouble *a_loc,
     75                IssmPDouble *rhs) {
     76        /*Initialize mumps: {{{*/
     77        ISSM_MPI_Comm   comm=IssmComm::GetComm();
     78        DMUMPS_STRUC_C theMumpsStruc;
     79        MumpsSettings(theMumpsStruc,comm);
     80        MumpsInit(theMumpsStruc);
     81        /*}}}*/
     82        // now setup the rest of theMumpsStruc
     83        theMumpsStruc.n=n;
     84        theMumpsStruc.nz=nnz;
     85        theMumpsStruc.nz_loc=local_nnz;
     86        theMumpsStruc.irn_loc=irn_loc;
     87        theMumpsStruc.jcn_loc=jcn_loc;
     88        theMumpsStruc.a_loc=a_loc;
     89        theMumpsStruc.rhs=rhs;
     90        theMumpsStruc.nrhs=1;
     91        theMumpsStruc.lrhs=1;
     92        /*Solve system: {{{*/
     93        MumpsAnalyze(theMumpsStruc);
     94        MumpsFactorize(theMumpsStruc);
     95        MumpsBacksubstitute(theMumpsStruc);
     96        /*}}}*/
     97        MumpsFinalize(theMumpsStruc);
     98}
     99
     100#ifdef _HAVE_ADOLC_
     101// prototype for active variant
     102void MumpsSolve(int n,
     103                int nnz,
     104                int local_nnz,
     105                int* irn_loc,
     106                int* jcn_loc,
     107                IssmDouble *a_loc,
     108                IssmDouble *rhs);
     109#endif
     110
     111void 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){ /*{{{*/
    27112
    28113        /*Variables: {{{*/
    … …  
    36121        int        *irn_loc = NULL;
    37122        int        *jcn_loc = NULL;
    38         IssmPDouble *a_loc   = NULL;
     123        IssmDouble *a_loc   = NULL;
    39124        int         count;
    40125        int         lower_row;
    41126        int         upper_row;
    42         IssmPDouble* rhs=NULL;
     127        IssmDouble* rhs=NULL;
    43128        int*        recvcounts=NULL;
    44129        int*        displs=NULL;
    … …  
    54139        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");
    55140        /*}}}*/
    56         /*Initialize mumps: {{{*/
    57         DMUMPS_STRUC_C id;
    58         id.job          = JOB_INIT;
    59         id.par          = 1; 
    60         id.sym          = 0;
    61         id.comm_fortran = MPI_Comm_c2f(comm);
    62         dmumps_c(&id);
    63         /*}}}*/
    64         /*Control statements:{{{ */
    65         id.icntl[1-1] = 6; //error verbose
    66         id.icntl[2-1] = 1; //std verbose
    67         id.icntl[4-1] = 4; //verbose everything
    68         id.icntl[5-1] = 0;
    69         id.icntl[18-1] = 3;
    70 
    71         id.icntl[20-1] = 0;
    72         id.icntl[21-1] = 0;
    73         id.icntl[30-1] = 0;
    74         /*}}}*/
    75141        /*Initialize matrix:{{{ */
    76         id.n=Kff_M;
    77142
    78143        /*figure out number of non-zero entries: */
    … …  
    86151        ISSM_MPI_Reduce(&local_nnz,&nnz,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,comm);
    87152        ISSM_MPI_Bcast(&nnz,1,ISSM_MPI_INT,0,comm);
    88         id.nz=nnz;
    89         id.nz_loc=local_nnz;
    90153
    91154        /*Allocate: */
    … …  
    93156                irn_loc=xNew<int>(local_nnz);
    94157                jcn_loc=xNew<int>(local_nnz);
    95                 a_loc=xNew<IssmPDouble>(local_nnz);
     158                a_loc=xNew<IssmDouble>(local_nnz);
    96159        }
    97160
    … …  
    109172                }
    110173        }
    111         id.irn_loc=irn_loc;
    112         id.jcn_loc=jcn_loc;
    113         id.a_loc=a_loc;
    114 
    115174        /*Deal with right hand side. We need to ISSM_MPI_Gather it onto cpu 0: */
    116         rhs=xNew<IssmPDouble>(pf_M);
     175        rhs=xNew<IssmDouble>(pf_M);
    117176
    118177        recvcounts=xNew<int>(num_procs);
    … …  
    127186        /*Gather:*/
    128187        ISSM_MPI_Gatherv(pf, pf_m, ISSM_MPI_DOUBLE, rhs, recvcounts, displs, ISSM_MPI_DOUBLE,0,comm);
    129         id.rhs=rhs;
    130         id.nrhs=1;
    131         id.lrhs=1;
    132 
    133         /*}}}*/
    134         /*Solve system: {{{*/
    135         id.job = JOB_SOLVE;
    136         dmumps_c (&id);
    137         /*}}}*/
     188
     189        MumpsSolve(Kff_M,
     190                   nnz,
     191                   local_nnz,
     192                   irn_loc,
     193                   jcn_loc,
     194                   a_loc,
     195                   rhs);
     196
    138197        /*Now scatter from cpu 0 to all other cpus: {{{*/
    139198        ISSM_MPI_Scatterv( rhs, recvcounts, displs, ISSM_MPI_DOUBLE, uf, uf_m, ISSM_MPI_DOUBLE, 0, comm);
    … …  
    141200        /*}}}*/
    142201        /*Cleanup: {{{*/
    143         id.job = JOB_END;
    144         dmumps_c (&id);
    145 
    146202        xDelete<int>(irn_loc);
    147203        xDelete<int>(jcn_loc);
    148         xDelete<IssmPDouble>(a_loc);
    149         xDelete<IssmPDouble>(rhs);
     204        xDelete<IssmDouble>(a_loc);
     205        xDelete<IssmDouble>(rhs);
    150206        xDelete<int>(recvcounts);
    151207        xDelete<int>(displs);
    152 
    153208        /*}}}*/
    154209} /*}}}*/
    155210
    156211#ifdef _HAVE_ADOLC_
     212
     213int mumpsSolveEDF(int iArrLength, int* iArr, int nPlusNz /* we can ignore it*/, double* dp_x, int m, double* dp_y) {
     214  // unpack parameters
     215  int n=iArr[0];
     216  int nz=iArr[1];
     217  int *irn=new int[nz];
     218  int *jcn=new int[nz];
     219  double *A=new double[nz];
     220  for (int i=0;i<nz;++i) {
     221    irn[i]=iArr[2+i];
     222    jcn[i]=iArr[2+nz+i];
     223    A[i]=dp_x[i];
     224  }
     225  double *rhs_sol=new double[n];
     226  for (int i=0;i<n;++i) {
     227    rhs_sol[i]=dp_x[nz+i];
     228  }
     229  mumpsSolve(n,nz,irn,jcn,A,rhs_sol);
     230  for (int i=0;i<m;++i) {
     231    dp_y[i]=rhs_sol[i];
     232  }
     233  return 0;
     234}
     235
     236void MumpsSolve(int n,
     237                int nnz,
     238                int local_nnz,
     239                int* irn_loc,
     240                int* jcn_loc,
     241                IssmDouble *a_loc,
     242                IssmDouble *rhs) {
     243  int packedDimsSparseArrLength=1+1+1+local_nnz+local_nnz;
     244  int *packedDimsSparseArr=xNew<int>(packedDimsSparseArrLength);
     245  packedDimsSparseArr[0]=n;
     246  packedDimsSparseArr[1]=nnz;
     247  packedDimsSparseArr[2]=local_nnz;
     248  for (int i=0;i<local_nnz;++i) {
     249    packedDimsSparseArr[3+i]=irn_loc[i];
     250    packedDimsSparseArr[3+local_nnz+i]=jcn_loc[i];
     251  }
     252  ensureContiguousLocations(local_nnz+n);
     253  adouble *pack_A_rhs=xNew<IssmDouble>(local_nnz+n);
     254  for (int i=0;i<local_nnz;++i) {
     255    pack_A_rhs[i]=a_loc[i];
     256  }
     257  for (int i=0;i<n;++i) {
     258    pack_A_rhs[local_nnz+i]=rhs[i];
     259  }
     260  double *passivePack_A_rhs=xNew<IssmPDouble>(local_nnz+n);
     261  double *passiveSol=xNew<IssmPDouble>(n);
     262  ensureContiguousLocations(n);
     263  adouble *sol=xNew<IssmDouble>(n);
     264  call_ext_fct(ourEDF_p,
     265               packedDimsSparseArrLength, packedDimsSparseArr,
     266               local_nnz+n, passivePack_A_rhs, pack_A_rhs,
     267               n, passiveSol,sol);
     268  for (int i=0;i<n;++i) {
     269    rhs[i]=sol[i];
     270  }
     271  xDelete(sol);
     272  xDelete(passiveSol);
     273  xDelete(passivePack_A_rhs);
     274  xDelete(pack_A_rhs);
     275  xDelete(packedDimsSparseArr);
     276}
     277
     278int fos_forward_mumpsSolveEDF(int iArrLength, int* iArr, int nPlusNz /* we can ignore it*/,
     279                              double *dp_x, double *dp_X, int m, double *dp_y, double *dp_Y) {
     280  // unpack parameters
     281  int n=iArr[0];
     282  int nz=iArr[1];
     283  int *irn=new int[nz];
     284  int *jcn=new int[nz];
     285  double *A=new double[nz];
     286  for (int i=0;i<nz;++i) {
     287    irn[i]=iArr[2+i];
     288    jcn[i]=iArr[2+nz+i];
     289    A[i]=dp_x[i];
     290  }
     291  double *rhs_sol=new double[n];
     292  for (int i=0;i<n;++i) {
     293    rhs_sol[i]=dp_x[nz+i];
     294  }
     295  DMUMPS_STRUC_C id;
     296  id.par = 1; // one processor=sequential code
     297  id.sym = 0; // asymmetric
     298  id.job = JOB_INIT;
     299  dmumps_c(&id);
     300
     301  id.icntl[1-1] = 6; //error verbose
     302  id.icntl[2-1] = 0; //std verbose
     303  id.icntl[3-1] = 0; //
     304  id.icntl[4-1] = 0; //
     305  id.icntl[5-1] = 0; // matrix is assembled
     306  id.icntl[18-1] = 0; // centralized
     307  id.icntl[20-1] = 0; // rhs is dense and centralized
     308  id.icntl[21-1] = 0; // solution is centralized
     309  id.n=n;
     310  id.nz=nz;
     311  id.irn=irn;
     312  id.jcn=jcn;
     313  id.a=A;
     314  id.job = JOB_ANALYSIS;
     315  dmumps_c(&id);
     316  id.job = JOB_FACTORIZATION;
     317  dmumps_c (&id);
     318  // solve the orifginal system
     319  id.rhs=rhs_sol;
     320  id.nrhs=1;
     321  id.lrhs=1;
     322  id.job = JOB_BACKSUBST;
     323  dmumps_c (&id);
     324  for (int i=0;i<m;++i) {
     325    dp_y[i]=rhs_sol[i];
     326  }
     327  // solve for the derivative
     328  for (int i=0;i<n;++i) {
     329    rhs_sol[i]=dp_X[nz+i];
     330  }
     331  for (int i=0;i<nz;++i) {
     332    rhs_sol[irn[i]-1]-=dp_X[i]*dp_y[jcn[i]-1];
     333  }
     334  dmumps_c (&id);
     335  for (int i=0;i<m;++i) {
     336    dp_Y[i]=rhs_sol[i];
     337  }
     338  id.job = JOB_END;
     339  dmumps_c (&id);
     340  return 3;
     341}
     342
     343int fos_reverse_mumpsSolveEDF(int iArrLength, int* iArr,
     344                              int m, double *dp_U,
     345                              int nPlusNz, double *dp_Z,
     346                              double *dp_x, double *dp_y) {
     347  // unpack parameters
     348  int n=iArr[0];
     349  int nz=iArr[1];
     350  int *irn=new int[nz];
     351  int *jcn=new int[nz];
     352  double *A=new double[nz];
     353  for (int i=0;i<nz;++i) {
     354    irn[i]=iArr[2+i];
     355    jcn[i]=iArr[2+nz+i];
     356    A[i]=dp_x[i];
     357  }
     358  DMUMPS_STRUC_C id;
     359  id.par = 1; // one processor=sequential code
     360  id.sym = 0; // asymmetric
     361  id.job = JOB_INIT;
     362  dmumps_c(&id);
     363
     364  id.icntl[1-1] = 6; //error verbose
     365  id.icntl[2-1] = 0; //std verbose
     366  id.icntl[3-1] = 0; //
     367  id.icntl[4-1] = 0; //
     368  id.icntl[5-1] = 0; // matrix is assembled
     369  id.icntl[9-1] = 0; //solve for the transpose
     370  id.icntl[18-1] = 0; // centralized
     371  id.icntl[20-1] = 0; // rhs is dense and centralized
     372  id.icntl[21-1] = 0; // solution is centralized
     373  id.n=n;
     374  id.nz=nz;
     375  id.irn=irn;
     376  id.jcn=jcn;
     377  id.a=A;
     378  id.job = JOB_ANALYSIS;
     379  dmumps_c(&id);
     380  id.job = JOB_FACTORIZATION;
     381  dmumps_c (&id);
     382  double *rhs_sol=new double[n];
     383  for (int i=0;i<n;++i) {
     384    rhs_sol[i]=dp_U[i];
     385  }
     386  id.rhs=rhs_sol;
     387  id.nrhs=1;
     388  id.lrhs=1;
     389  id.job = JOB_BACKSUBST;
     390  dmumps_c (&id);
     391  // update the adhoint of the rhs:
     392  for (int i=0;i<m;++i) {
     393    dp_Z[nz+i]+=rhs_sol[i];
     394  }
     395  // update the adjoint of the matrix:
     396  for (int i=0;i<nz;++i) {
     397    dp_Z[i]+=-dp_U[irn[i]-1]*dp_y[jcn[i]-1];
     398  }
     399  return 3;
     400}
     401
    157402void 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){ /*{{{*/
    158403        _error_("not supported yet!");
    159404} /*}}}*/
     405
    160406#endif
Note: See TracChangeset for help on using the changeset viewer.