Changeset 23251


Ignore:
Timestamp:
09/11/18 09:19:05 (7 years ago)
Author:
Mathieu Morlighem
Message:

CHG: final changes from CoDiPack

Location:
issm/trunk-jpl/src/c
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r23242 r23251  
    1717#include "../shared/Enum/Enum.h"
    1818#include "../analyses/analyses.h"
     19
     20#if _HAVE_CODIPACK_
     21#include <sstream> // for output of the CoDiPack tape
     22extern CoDi_global codi_global;
     23#endif
    1924
    2025/*module includes: {{{*/
    … …  
    18511856                        dependents=xNew<IssmPDouble>(num_dependents);
    18521857
     1858                        #if defined(_HAVE_CODIPACK_)
     1859                        auto& tape_codi = IssmDouble::getGlobalTape();
     1860                        #endif
     1861
    18531862                        /*Go through our dependent variables, and compute the response:*/
    18541863                        for(int i=0;i<dependent_objects->Size();i++){
    … …  
    18561865                                dep->Responsex(&output_value,this);
    18571866                                if (my_rank==0) {
    1858                                         output_value>>=dependents[i];
     1867                                        #if defined(_HAVE_CODIPACK_)
     1868                                                tape_codi.registerOutput(output_value);
     1869                                                dependents[i] = output_value.getValue();
     1870                                                codi_global.output_indices.push_back(output_value.getGradientData());
     1871                                        #else
     1872                                                output_value>>=dependents[i];
     1873                                        #endif
    18591874                                }
    18601875                        }

  • issm/trunk-jpl/src/c/classes/IoModel.cpp

    r23174 r23251  
    2727#include "../shared/io/io.h"
    2828#include "../shared/shared.h"
     29
     30#ifdef _HAVE_CODIPACK_
     31extern CoDi_global codi_global;
     32#include <sstream> // for output of the CoDiPack tape
     33#endif
    2934
    3035/*IoConstant class and methods*/
    … …  
    431436        if(trace || (autodiff && !iscontrol)){
    432437
    433                 #ifdef _HAVE_ADOLC_
     438                #ifdef _HAVE_AD_
     439                // FIXME codi here we should be able to execute codi version as normal
    434440                this->FetchData(&num_independent_objects,"md.autodiff.num_independent_objects");
    435441                if(num_independent_objects){
    … …  
    460466                #else
    461467                /*if we asked for AD computations, we have a problem!: */
    462                 _error_("Cannot carry out AD mode computations without support of ADOLC compiled in!");
     468                _error_("Cannot carry out AD mode computations without support of ADOLC or CoDiPack compiled in!");
    463469                #endif
    464470        }
    … …  
    17021708        if(X) Xcount=*pXcount;
    17031709
    1704         #ifdef _HAVE_ADOLC_ //cannot come here unless you are running AD mode, from DeclaredIndependents:
     1710        #ifdef _HAVE_AD_ //cannot come here unless you are running AD mode, from DeclaredIndependents:
    17051711
    17061712        /*output: */
    … …  
    17191725                /*Now, before we even broadcast this to other nodes, declare the scalar  as an independent variable!. If we
    17201726                 *have been supplied an X vector, use it instead of what we just read: */
    1721                 if(X){
    1722                         scalar<<=X[Xcount];
    1723                 }
    1724                 else{
    1725                         scalar<<=pscalar;
    1726                 }
     1727                #if defined(_HAVE_CODIPACK_)
     1728                        // FIXME codi here we just assign instead of using "operator <<="
     1729                        if(X){
     1730                                scalar=X[Xcount];
     1731                        } else {
     1732                                scalar=pscalar;
     1733                        }
     1734                        auto& tape_codi = IssmDouble::getGlobalTape();
     1735                        tape_codi.registerInput(scalar);
     1736                        codi_global.input_indices.push_back(scalar.getGradientData());
     1737                #else
     1738                        if(X){
     1739                                scalar<<=X[Xcount];
     1740                        }
     1741                        else{
     1742                                scalar<<=pscalar;
     1743                        }
     1744                #endif
    17271745        }
    17281746
    … …  
    17851803                        /*Now, before we even broadcast this to other nodes, declare the whole matrix as a independent variable!
    17861804                          If we have been supplied an X vector, use it instead of what we just read: */
    1787                         if(X){
    1788                                 for(int i=0;i<M*N;i++) matrix[i]<<=X[Xcount+i];  /*<<= ADOLC overloaded operator to declare independent*/
    1789                         }
    1790                         else{
    1791                                 for(int i=0;i<M*N;i++) matrix[i]<<=buffer[i];
    1792                         }
     1805                        #if defined(_HAVE_CODIPACK_)
     1806                                // FIXME codi here we just assign instead of using "operator <<="
     1807                                auto& tape_codi = IssmDouble::getGlobalTape();
     1808                                if(X){
     1809                                        for (int i=0;i<M*N;i++) {
     1810                                                matrix[i]=X[Xcount+i];
     1811                                                tape_codi.registerInput(matrix[i]);
     1812                                                codi_global.input_indices.push_back(matrix[i].getGradientData());
     1813                                        }
     1814                                }
     1815                                else{
     1816                                        for (int i=0;i<M*N;i++) {
     1817                                                matrix[i]=buffer[i];
     1818                                                tape_codi.registerInput(matrix[i]);
     1819                                                codi_global.input_indices.push_back(matrix[i].getGradientData());
     1820                                        }
     1821                                }
     1822                        #else /*ADOLC*/
     1823                                if(X){
     1824                                        for(int i=0;i<M*N;i++) matrix[i]<<=X[Xcount+i];  /*<<= ADOLC overloaded operator to declare independent*/
     1825                                }
     1826                                else{
     1827                                        for(int i=0;i<M*N;i++) matrix[i]<<=buffer[i];
     1828                                }
     1829                        #endif
    17931830                }
    17941831                ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    … …  
    17981835        else _error_("cannot declare the independent variable \"" << data_name <<  "\" if it's empty!");
    17991836
     1837        // FIXME codi is that at all relevant to CoDiPack or can we simply assume the same?
     1838       
    18001839        /*Add to data as independent*/
    18011840        this->AddDataIndependent(new IoData(matrix,code,layout,M,N,data_name));
    … …  
    25622601
    25632602        if(trace || (autodiff && !iscontrol)){
    2564                 #ifdef _HAVE_ADOLC_
     2603
     2604                #if defined(_HAVE_ADOLC_)
    25652605                /*Retrieve parameters: */
    25662606                this->FetchData(&keep,"md.autodiff.keep");
    … …  
    25792619                int skipFileDeletion=1;
    25802620                trace_on(my_rank,keepTaylors,reCast<size_t>(obufsize),reCast<size_t>(lbufsize),reCast<size_t>(cbufsize),reCast<size_t>(tbufsize),skipFileDeletion);
     2621
     2622                #elif defined(_HAVE_CODIPACK_)
     2623                //fprintf(stderr, "*** Codipack IoModel::StartTrace\n");
     2624                /*
     2625                 * FIXME codi
     2626                 * - ADOL-C variant uses fine grained tracing with various arguments
     2627                 * - ADOL-C variant sets a garbage collection parameter for its tape
     2628                 * -> These parameters are not read for the CoDiPack ISSM version!
     2629                 */
     2630                auto& tape_codi = IssmDouble::getGlobalTape();
     2631                tape_codi.setActive();
     2632                #if _AD_TAPE_ALLOC_
     2633                //alloc_profiler.Tag(StartInit, true);
     2634                IssmDouble x_t(1.0), y_t(1.0);
     2635                tape_codi.registerInput(y_t);
     2636                int codi_allocn = 0;
     2637                this->FetchData(&codi_allocn, AutodiffTapeAllocEnum);
     2638                for(int i = 0;i < codi_allocn;++i) {
     2639                        x_t = y_t * y_t;
     2640                }
     2641                /*
     2642                std::stringstream out_s;
     2643                IssmDouble::getGlobalTape().printStatistics(out_s);
     2644                _printf0_("StartTrace::Tape Statistics     : TapeAlloc count=[" << codi_allocn << "]\n" << out_s.str());
     2645                */
     2646                tape_codi.reset();
     2647                //alloc_profiler.Tag(FinishInit, true);
    25812648                #endif
    2582         }
    2583 
    2584 }
    2585 /*}}}*/
     2649                #endif
     2650        }
     2651
     2652}
     2653/*}}}*/

  • issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp

    r21615 r23251  
    1414#include "../../classes/Params/Parameters.h"
    1515#include "../adolc/adolcincludes.h"
     16#include "../codipack/codipackincludes.h"
    1617#include "./gslincludes.h"
    1718
    … …  
    251252/*}}}*/
    252253#endif
     254
     255#ifdef _HAVE_CODIPACK_
     256void SolverxSeq_codi_b(codi::DataStore* data) {/*{{{*/
     257  IssmDouble::TapeType& tape = IssmDouble::getGlobalTape();
     258
     259  IssmDouble::Real* valueATrans;
     260  IssmDouble::GradientData* indexATrans;
     261  IssmDouble::GradientData* indexB;
     262  IssmDouble::Real* valueX;
     263  IssmDouble::GradientData* indexX;
     264  int n;
     265
     266  data->getData(valueATrans);
     267  data->getData(indexATrans);
     268  data->getData(indexB);
     269  data->getData(valueX);
     270  data->getData(indexX);
     271  data->getData(n);
     272
     273
     274  // create the adjoint vector for x and reset the adjoint values on the tape
     275  IssmDouble::GradientValue* adjX = xNew<IssmDouble::GradientValue>(n);
     276  getVectorAdjoint(tape, indexX, adjX, n);
     277
     278  IssmDouble::GradientValue* sol  = xNew<IssmDouble::GradientValue>(n);
     279  SolverxSeq(sol, valueATrans, adjX, n);
     280
     281  updateVectorAdjoint(tape, indexB, sol, n);
     282  for(int i=0; i<n; ++i) {
     283    for (int j=0; j<n; ++j) {
     284      // we access the transposed matrix here because we stored the indices in a transposed way
     285      updateAdjoint(tape, indexATrans[i*n+j], -sol[j]*valueX[i]);
     286    }
     287  }
     288
     289  xDelete(sol);
     290  xDelete(adjX);
     291}
     292/*}}}*/
     293void SolverxSeq_codi_delete(codi::DataStore* data) {/*{{{*/
     294  IssmDouble::Real* valueATrans;
     295  IssmDouble::GradientData* indexATrans;
     296  IssmDouble::GradientData* indexB;
     297  IssmDouble::Real* valueX;
     298  IssmDouble::GradientData* indexX;
     299  int n;
     300
     301  data->getData(valueATrans);
     302  data->getData(indexATrans);
     303  data->getData(indexB);
     304  data->getData(valueX);
     305  data->getData(indexX);
     306  data->getData(n);
     307
     308  xDelete(valueATrans);
     309  xDelete(indexATrans);
     310  xDelete(indexB);
     311  xDelete(valueX);
     312  xDelete(indexX);
     313}
     314/*}}}*/
     315void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
     316  IssmDouble::TapeType& tape = IssmDouble::getGlobalTape();
     317  codi::DataStore* dataHandler = NULL;
     318
     319  if(tape.isActive()) {
     320    dataHandler = new codi::DataStore();
     321
     322    // create the index vector and the double data for A and B
     323    IssmDouble::Real* valueATrans = xNew<IssmDouble::Real>(n*n);
     324    IssmDouble::GradientData* indexATrans = xNew<IssmDouble::GradientData>(n*n);
     325
     326    // read the data for matrix in a transposed fashion
     327          for (int i=0; i<n; ++i) {
     328      for (int j=0; j<n; ++j) {
     329        getPrimalAndGradData(A[i*n+j], valueATrans[j*n+i], indexATrans[j*n+i]);
     330      }
     331    }
     332
     333    // read the data from B (primal values are not required vor B
     334    IssmDouble::GradientData* indexB = xNew<IssmDouble::GradientData>(n);
     335    getVectorGradData(B, indexB, n);
     336
     337    dataHandler->addData(valueATrans);
     338    dataHandler->addData(indexATrans);
     339    dataHandler->addData(indexB);
     340  }
     341
     342  // unpack the primal values from the matrix and the vector
     343  IssmDouble::Real* valueA = xNew<IssmDouble::Real>(n*n);
     344  IssmDouble::Real* valueB = xNew<IssmDouble::Real>(n);
     345  // read the data from A and B
     346  getVectorPrimal(A, valueA, n*n);
     347  getVectorPrimal(B, valueB, n);
     348
     349  // create the placeholder for X and solve the system
     350  IssmDouble::Real* valueX = xNew<IssmDouble::Real>(n);
     351  SolverxSeq(valueX, valueA, valueB, n);
     352
     353  // pack the values into x
     354  setVectorPrimal(X, valueX, n);
     355
     356  if(tape.isActive()) {
     357    // create the index vector X and register x as active variables
     358    IssmDouble::GradientData* indexX = xNew<IssmDouble::GradientData>(n);
     359    registerVector(X, indexX, n);
     360
     361    dataHandler->addData(valueX);
     362    dataHandler->addData(indexX);
     363
     364    // store other arguments
     365    dataHandler->addData(n);
     366
     367    tape.pushExternalFunction(&SolverxSeq_codi_b, dataHandler, &SolverxSeq_codi_delete);
     368  } else {
     369    // if the tape is active valueX is stored in the dataHandler and deleted in the reverse sweep
     370    xDelete(valueX);
     371  }
     372
     373  xDelete(valueB);
     374  xDelete(valueA);
     375}
     376/*}}}*/
     377void DenseGslSolve(/*output*/ IssmDouble** px,/*stiffness matrix:*/ IssmDouble* Kff, int Kff_M, int Kff_N, /*right hand side load vector: */ IssmDouble* pf, int pf_M, Parameters* parameters){ /*{{{*/
     378
     379        /*Intermediary: */
     380
     381        if(Kff_N!=pf_M)_error_("Right hand side vector of size " << pf_M << ", when matrix is of size " << Kff_M << "-" << Kff_N << " !");
     382        if(Kff_M!=Kff_N)_error_("Stiffness matrix should be square!");
     383
     384        IssmDouble *x  = xNew<IssmDouble>(Kff_N,"t");
     385
     386        SolverxSeq(x,Kff,pf,Kff_N,parameters);
     387
     388        /*allocate output pointers: */
     389        *px=x;
     390}
     391/*}}}*/
     392#endif

  • issm/trunk-jpl/src/c/toolkits/mpi/issmmpi.cpp

    r20204 r23251  
    6060# endif
    6161#else
    62 # ifdef _HAVE_ADOLC_
     62# ifdef _HAVE_AD_
    6363  if (sendtype==ISSM_MPI_DOUBLE) {
    6464    IssmDouble* activeSendBuf=(IssmDouble*)sendbuf;
    … …  
    9797#else
    9898  assert(sendcount==recvcounts[0]); // we handle only identical representations
    99 # ifdef _HAVE_ADOLC_
     99# ifdef _HAVE_AD_
    100100  if (sendtype==ISSM_MPI_DOUBLE) {
    101101    IssmDouble* activeSendBuf=(IssmDouble*)sendbuf;
    … …  
    129129# endif
    130130#else
    131 # ifdef _HAVE_ADOLC_
     131#ifdef _HAVE_AD_
    132132  if (datatype==ISSM_MPI_DOUBLE) {
    133133    IssmDouble* activeSendBuf=(IssmDouble*)sendbuf;
    … …  
    215215int ISSM_MPI_Finalize(void){  /*{{{*/
    216216
    217   int rc=0;
    218 #ifdef _HAVE_MPI_
    219 # ifdef _HAVE_AMPI_
    220   rc=AMPI_Finalize_NT();
    221 # else
    222   rc=MPI_Finalize();
    223 # endif
    224 #endif
     217        int rc=0;
     218        #ifdef _HAVE_MPI_
     219                #ifdef _HAVE_AMPI_
     220                        #ifdef _HAVE_ADJOINTMPI_
     221                                rc=AMPI_Finalize();
     222                        #else
     223                                rc=AMPI_Finalize_NT();
     224                        #endif
     225                #else
     226                  rc=MPI_Finalize();
     227                #endif
     228        #endif
    225229  return rc;
    226230}/*}}}*/
    … …  
    250254# endif
    251255#else
    252 # ifdef _HAVE_ADOLC_
     256# ifdef _HAVE_AD_
    253257  if (sendtype==ISSM_MPI_DOUBLE) {
    254258    IssmDouble* activeSendBuf=(IssmDouble*)sendbuf;
    … …  
    290294#else
    291295  assert(sendcnt==recvcnts[0]); // we handle only identical representations
    292 # ifdef _HAVE_ADOLC_
     296#ifdef _HAVE_AD_
    293297  if (sendtype==ISSM_MPI_DOUBLE) {
    294298    IssmDouble* activeSendBuf=(IssmDouble*)sendbuf;
    … …  
    304308int ISSM_MPI_Init(int *argc, char ***argv){  /*{{{*/
    305309
    306   int rc=0;
    307 #ifdef _HAVE_MPI_
    308 # ifdef _HAVE_AMPI_
    309   rc=AMPI_Init_NT(argc,argv);
    310 # else
    311   rc=MPI_Init(argc,argv);
    312 # endif
    313 #endif
     310        int rc=0;
     311        #ifdef _HAVE_MPI_
     312                #ifdef _HAVE_AMPI_
     313                        #ifdef _HAVE_ADJOINTMPI_
     314                                rc=AMPI_Init(argc,argv);
     315                        #else
     316                                rc=AMPI_Init_NT(argc,argv);
     317                        #endif
     318                #else
     319                        rc=MPI_Init(argc,argv);
     320                #endif
     321        #endif
    314322  return rc;
    315323}/*}}}*/
    … …  
    324332               source,
    325333               tag,
     334                         #ifndef _HAVE_ADJOINTMPI_
    326335               AMPI_FROM_SEND, // as long as there are no other variants
     336                         #endif
    327337               comm,
    328338               status);
    … …  
    365375# endif
    366376#else
    367 # ifdef _HAVE_ADOLC_
     377# ifdef _HAVE_AD_
    368378  if (datatype==ISSM_MPI_DOUBLE) {
    369379    IssmDouble* activeSendBuf=(IssmDouble*)sendbuf;
    … …  
    402412# endif
    403413#else
    404 # ifdef _HAVE_ADOLC_
     414# ifdef _HAVE_AD_
    405415  if (sendtype==ISSM_MPI_DOUBLE) {
    406416    IssmDouble* activeSendBuf=(IssmDouble*)sendbuf;
    … …  
    442452#else
    443453  assert(sendcnts[0]==recvcnt); // we handle only identical representations
    444 # ifdef _HAVE_ADOLC_
     454# ifdef _HAVE_AD_
    445455  if (sendtype==ISSM_MPI_DOUBLE) {
    446456    IssmDouble* activeSendBuf=(IssmDouble*)(sendbuf)+displs[0];
    … …  
    464474               dest,
    465475               tag,
     476                         #ifndef _HAVE_ADJOINTMPI_
    466477               AMPI_TO_RECV, // as long as there are no other variants
     478                         #endif
    467479               comm);
    468480# else
    … …  
    492504#ifdef _HAVE_ADOLC_
    493505  ensureContiguousLocations(aSize);
     506#else
     507  fprintf(stderr, "*** Codipack ISSM_MPI_ContiguousInAdolc()\n");
    494508#endif
    495509}/*}}}*/

  • issm/trunk-jpl/src/c/toolkits/mpi/issmmpi.h

    r22585 r23251  
    1919        /*Include header files: {{{*/
    2020        #ifdef _HAVE_AMPI_
    21                 #include <ampi/ampi.h>
     21                #ifdef _HAVE_ADJOINTMPI_
     22                        #include <ampi_tape.hpp>
     23                #else
     24                        #include <ampi/ampi.h>
     25                #endif
    2226        #elif  _HAVE_PETSC_MPI_ // Petsc now hides there MPI header. It can be reached through Petsc.
    2327                #include <petsc.h>
    … …  
    3539        // data types
    3640        #define ISSM_MPI_CHAR       MPI_CHAR
    37         #ifdef _HAVE_AMPI_
     41        #if defined(_HAVE_AMPI_)
    3842                #define ISSM_MPI_DOUBLE    AMPI_ADOUBLE // corresponds to IssmDouble
     43        #elif defined(_HAVE_ADJOINTMPI_)
     44                #define ISSM_MPI_DOUBLE    AMPI_DOUBLE // corresponds to IssmDouble
    3945        #else
    4046                #define ISSM_MPI_DOUBLE    MPI_DOUBLE  // corresponds to IssmDouble

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

    r23066 r23251  
    1818#include "../mpi/issmmpi.h"
    1919#include "../adolc/adolcincludes.h"
     20#include "../codipack/codipackincludes.h"
    2021#include "../issm/SparseRow.h"
    2122#include "./mumpsincludes.h"
    … …  
    9899}
    99100
    100 #ifdef _HAVE_ADOLC_
     101#ifdef _HAVE_AD_
    101102// prototype for active variant
    102103void MumpsSolve(int n,
    … …  
    354355
    355356#ifdef _HAVE_ADOLC_
    356 
    357 int mumpsSolveEDF(int iArrLength, int* iArr, int /* ignored */, IssmPDouble* dp_x, int /* ignored */, IssmPDouble* dp_y) {
     357int mumpsSolveEDF(int iArrLength, int* iArr, int /* ignored */, IssmPDouble* dp_x, int /* ignored */, IssmPDouble* dp_y){/*{{{*/
    358358  // unpack parameters
    359359  int n=iArr[0];
    … …  
    381381  xDelete(local_irn);
    382382  return 0;
    383 }
    384 
    385 void MumpsSolve(int n,int nnz,int local_nnz,int* irn_loc,int* jcn_loc,IssmDouble *a_loc,IssmDouble *rhs,Parameters* parameters){
     383}/*}}}*/
     384void MumpsSolve(int n,int nnz,int local_nnz,int* irn_loc,int* jcn_loc,IssmDouble *a_loc,IssmDouble *rhs,Parameters* parameters){/*{{{*/
    386385  int packedDimsSparseArrLength=1+1+1+local_nnz+local_nnz;
    387386  int *packedDimsSparseArr=xNew<int>(packedDimsSparseArrLength);
    … …  
    425424  xDelete(pack_A_rhs);
    426425  xDelete(packedDimsSparseArr);
    427 }
    428 
    429 int fos_reverse_mumpsSolveEDF(int iArrLength, int* iArr,
    430                               int m, IssmPDouble *dp_U,
    431                               int nPlusNz, IssmPDouble *dp_Z,
    432                               IssmPDouble *dp_x, IssmPDouble *dp_y) {
     426}/*}}}*/
     427int fos_reverse_mumpsSolveEDF(int iArrLength, int* iArr, int m, IssmPDouble *dp_U, int nPlusNz, IssmPDouble *dp_Z, IssmPDouble *dp_x, IssmPDouble *dp_y) {/*{{{*/
    433428  // unpack parameters
    434429  int n=iArr[0];
    … …  
    481476  xDelete(local_irn);
    482477  return 3;
    483 }
    484 
    485 #endif
     478}/*}}}*/
     479#endif
     480
     481#ifdef _HAVE_CODIPACK_
     482void MumpsSolve_codi_b(codi::DataStore* data) {/*{{{*/
     483  IssmDouble::TapeType& tape = IssmDouble::getGlobalTape();
     484
     485  IssmDouble::Real* valueATrans;
     486  IssmDouble::GradientData* indexATrans;
     487  int* irnATrans;
     488  int* jcnATrans;
     489  IssmDouble::GradientData* indexB;
     490  IssmDouble::Real* valueX;
     491  IssmDouble::GradientData* indexX;
     492  int n;
     493  int nnz;
     494  int local_nnz;
     495  Parameters* parameters;
     496
     497  data->getData(valueATrans);
     498  data->getData(indexATrans);
     499  data->getData(irnATrans);
     500  data->getData(jcnATrans);
     501  data->getData(indexB);
     502  data->getData(valueX);
     503  data->getData(indexX);
     504  data->getData(n);
     505  data->getData(nnz);
     506  data->getData(local_nnz);
     507  data->getData(parameters);
     508
     509  // create the adjoint vector for x and reset the adjoint values on the tape
     510  IssmDouble::GradientValue* adjX = xNew<IssmDouble::GradientValue>(n);
     511  getVectorAdjoint(tape, indexX, adjX, n);
     512
     513  MumpsSolve(n, nnz, local_nnz, irnATrans, jcnATrans, valueATrans, adjX, parameters);
     514  // adjX contains now the solution
     515
     516  updateVectorAdjoint(tape, indexB, adjX, n);
     517  for(int i=0; i<local_nnz; ++i) {
     518    // we access the transposed matrix here because we stored the indices in a transposed way
     519    // -1 is substracted because jcn and irn are stored with fortran indexing
     520    updateAdjoint(tape, indexATrans[i], -adjX[jcnATrans[i]-1]*valueX[irnATrans[i]-1]);
     521  }
     522
     523  xDelete(adjX);
     524}
     525/*}}}*/
     526void MumpsSolve_codi_delete(codi::DataStore* data) {/*{{{*/
     527  IssmDouble::Real* valueATrans;
     528  IssmDouble::GradientData* indexATrans;
     529  int* irnATrans;
     530  int* jcnATrans;
     531  IssmDouble::GradientData* indexB;
     532  IssmDouble::Real* valueX;
     533  IssmDouble::GradientData* indexX;
     534  int n;
     535  int nnz;
     536  int local_nnz;
     537  Parameters* parameters;
     538
     539  data->getData(valueATrans);
     540  data->getData(indexATrans);
     541  data->getData(irnATrans);
     542  data->getData(jcnATrans);
     543  data->getData(indexB);
     544  data->getData(valueX);
     545  data->getData(indexX);
     546  data->getData(n);
     547  data->getData(nnz);
     548  data->getData(local_nnz);
     549  data->getData(parameters);
     550
     551  xDelete(valueATrans);
     552  xDelete(indexATrans);
     553  xDelete(irnATrans);
     554  xDelete(jcnATrans);
     555  xDelete(indexB);
     556  xDelete(valueX);
     557  xDelete(indexX);
     558}
     559/*}}}*/
     560void MumpsSolve(int n,int nnz,int local_nnz,int* irn_loc,int* jcn_loc,IssmDouble *a_loc,IssmDouble *rhs,Parameters* parameters){/*{{{*/
     561  IssmDouble::TapeType& tape = IssmDouble::getGlobalTape();
     562  codi::DataStore* dataHandler = NULL;
     563
     564  if(tape.isActive()) {
     565    dataHandler = new codi::DataStore();
     566
     567    // create the index and double vector for the matrix
     568    IssmDouble::Real* valueATrans = xNew<IssmDouble::Real>(local_nnz);
     569    IssmDouble::GradientData* indexATrans = xNew<IssmDouble::GradientData>(local_nnz);
     570    int* irnATrans = xNew<int>(local_nnz);
     571    int* jcnATrans = xNew<int>(local_nnz);
     572
     573    // read the data for the matrix A in a transposed fashion
     574          for (int i=0; i<local_nnz; ++i) {
     575        getPrimalAndGradData(a_loc[i], valueATrans[i], indexATrans[i]);
     576        irnATrans[i]=jcn_loc[i];  // transposed store
     577        jcnATrans[i]=irn_loc[i];  // transposed store
     578    }
     579
     580    // create the index vector for a (primal values are not needed for a)
     581    IssmDouble::GradientData* indexB = xNew<IssmDouble::GradientData>(n);
     582    getVectorGradData(rhs, indexB, n);
     583
     584    dataHandler->addData(valueATrans);
     585    dataHandler->addData(indexATrans);
     586    dataHandler->addData(irnATrans);
     587    dataHandler->addData(jcnATrans);
     588    dataHandler->addData(indexB);
     589  }
     590
     591  // unpack the primal values from the matrix and the vector
     592  IssmDouble::Real* valueA = xNew<IssmDouble::Real>(local_nnz);
     593  IssmDouble::Real* valueB = xNew<IssmDouble::Real>(n);
     594  // read the data from A and B
     595  getVectorPrimal(a_loc, valueA, local_nnz);
     596  getVectorPrimal(rhs, valueB, n);
     597
     598  MumpsSolve(n, nnz, local_nnz, irn_loc, jcn_loc, valueA, valueB, parameters);
     599  // valueB contains now the solution
     600
     601  // pack the values into rhs
     602  setVectorPrimal(rhs, valueB, n);
     603
     604  if(tape.isActive()) {
     605    // create the index vector X and register x as active variables
     606    IssmDouble::GradientData* indexX = xNew<IssmDouble::GradientData>(n);
     607    registerVector(rhs, indexX, n);
     608
     609    dataHandler->addData(valueB); // contains the values from x
     610    dataHandler->addData(indexX);
     611
     612    // store other arguments
     613    dataHandler->addData(n);
     614    dataHandler->addData(nnz);
     615    dataHandler->addData(local_nnz);
     616    dataHandler->addData(parameters); // we assume here that parameters is still intact when the reverse run is called
     617
     618    tape.pushExternalFunction(&MumpsSolve_codi_b, dataHandler, &MumpsSolve_codi_delete);
     619  } else {
     620    // if the tape is active valueB is stored in the dataHandler and deleted in the reverse sweep
     621    xDelete(valueB);
     622  }
     623
     624  xDelete(valueA);
     625}
     626/*}}}*/
     627#endif
Note: See TracChangeset for help on using the changeset viewer.