Changeset 14877


Ignore:
Timestamp:
05/03/13 17:00:08 (12 years ago)
Author:
Eric.Larour
Message:

CHG: rearranged the adolc solvers to actually be under the gsl toolkit.

Location:
issm/trunk-jpl/src/c
Files:
3 added
5 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r14864 r14877  
    773773                                        ./toolkits/mumps/MpiDenseMumpsSolve.cpp
    774774#}}}
     775#Gsl sources  {{{
     776gsl_sources=      ./toolkits/gsl\
     777                                        ./toolkits/gsl/gslincludes.h\
     778                                        ./toolkits/gsl/DenseGslSolve.cpp
     779#}}}
    775780#Mpi sources  {{{
    776781mpi_sources= ./toolkits/mpi/mpiincludes.h\
    … …  
    875880endif
    876881
     882if GSL
     883issm_sources  +=  $(gsl_sources)
     884endif
     885
    877886if TRANSIENT
    878887issm_sources  +=  $(transient_sources)

  • issm/trunk-jpl/src/c/toolkits/issm/IssmDenseMat.h

    r14876 r14877  
    2323#include "../../include/macros.h"
    2424
    25 #ifdef _HAVE_ADOLC_
    26 #include "adolc.h"
    27 #endif
    28 
    2925#include <math.h>
     26
    3027
    3128/*}}}*/
    … …  
    273270                        IssmSeqVec<IssmDouble>* pf = NULL;
    274271                        IssmSeqVec<IssmDouble> *uf = NULL;
     272                        IssmDouble* x=NULL;
     273
     274                        /*Assume we are getting an IssmMpiVec in input, downcast: */
     275                        pf=(IssmSeqVec<IssmDouble>*)pfin;
    275276
    276277                        #ifdef _HAVE_GSL_
    277                                 /*Assume we are getting an IssmSeqVecVec in input, downcast: */
    278                                 pf=(IssmSeqVec<IssmDouble>*)pfin;
    279 
    280                                 /*Intermediary: */
    281                                 int N2;
    282 
    283                                 pf->GetSize(&N2);
    284 
    285                                 if(N!=N2)_error_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !");
    286                                 if(M!=N)_error_("Stiffness matrix should be square!");
    287                                
    288                                 #ifdef _HAVE_ADOLC_
    289                                         ensureContiguousLocations(N);
    290                                 #endif
    291                                 IssmDouble *x  = xNew<IssmDouble>(N);
    292 
    293                                 #ifdef _HAVE_ADOLC_
    294                                         SolverxSeq(x,this->matrix,pf->vector,N,parameters);
    295                                 #else
    296                                         SolverxSeq(x,this->matrix,pf->vector,N);
    297                                 #endif
    298 
    299                                 uf=new IssmSeqVec<IssmDouble>(x,N);
    300 
    301                                 xDelete(x);
    302 
    303                                 /*return: */
    304                                 return uf;
     278                        DenseGslSolve(/*output*/ &x,/*stiffness matrix:*/ this->matrix,this->M,this->N, /*right hand side load vector: */ pf->vector,pf->M,parameters);
     279
     280                        uf=new IssmSeqVec<IssmDouble>(x,this->N); xDelete(x);
     281                        return uf;
    305282                        #else
    306283                                _error_("GSL support not compiled in!");

  • issm/trunk-jpl/src/c/toolkits/issm/IssmSolver.cpp

    r14834 r14877  
    1010#include <cstring>
    1111
    12 #include "./IssmSolver.h"
    1312#include "../../shared/shared.h"
    1413#include "../../classes/classes.h"
    1514#include "../../include/include.h"
    1615#include "../../io/io.h"
    17 
    18 #ifdef _HAVE_GSL_
    19 #include <gsl/gsl_linalg.h>
    20 #endif
    2116
    2217void IssmSolve(IssmVec<IssmDouble>** pout,IssmMat<IssmDouble>* Kff, IssmVec<IssmDouble>* pf, Parameters* parameters){/*{{{*/
    … …  
    3025}
    3126/*}}}*/
    32 void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
    33 
    34         /*Allocate output*/
    35         double* X = xNew<double>(n);
    36 
    37         /*Solve*/
    38         SolverxSeq(X,A,B,n);
    39 
    40         /*Assign output pointer*/
    41         *pX=X;
    42 }
    43 /*}}}*/
    44 void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
    45 #ifdef _HAVE_GSL_
    46         /*GSL Matrices and vectors: */
    47         int              s;
    48         gsl_matrix_view  a;
    49         gsl_vector_view  b,x;
    50         gsl_permutation *p = NULL;
    51 //      for (int i=0; i<n*n; ++i) std::cout << "SolverxSeq A["<< i << "]=" << A[i] << std::endl;
    52 //      for (int i=0; i<n; ++i) std::cout << "SolverxSeq b["<< i << "]=" << B[i] << std::endl;
    53         /*A will be modified by LU decomposition. Use copy*/
    54         double* Acopy = xNew<double>(n*n);
    55         xMemCpy(Acopy,A,n*n);
    56 
    57         /*Initialize gsl matrices and vectors: */
    58         a = gsl_matrix_view_array (Acopy,n,n);
    59         b = gsl_vector_view_array (B,n);
    60         x = gsl_vector_view_array (X,n);
    61 
    62         /*Run LU and solve: */
    63         p = gsl_permutation_alloc (n);
    64         gsl_linalg_LU_decomp (&a.matrix, p, &s);
    65         gsl_linalg_LU_solve (&a.matrix, p, &b.vector, &x.vector);
    66 
    67         /*Clean up and assign output pointer*/
    68         xDelete(Acopy);
    69         gsl_permutation_free(p);
    70 #endif
    71 }
    72 /*}}}*/
    73 
    74 #ifdef _HAVE_ADOLC_
    75 int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y){ /*{{{*/
    76         SolverxSeq(y,x, x+m*m, m); // x is where the matrix starts, x+m*m is where the right-hand side starts
    77         return 0;
    78 } /*}}}*/
    79 int EDF_fos_forward_for_solverx(int n, IssmPDouble *inVal, IssmPDouble *inDeriv, int m, IssmPDouble *outVal, IssmPDouble *outDeriv) { /*{{{*/
    80 #ifdef _HAVE_GSL_
    81         //  for (int i=0; i<m*m; ++i) std::cout << "EDF_fos_forward_for_solverx A["<< i << "]=" << inVal[i] << std::endl;
    82         //  for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx b["<< i << "]=" << inVal[i+m*m] << std::endl;
    83         // the matrix will be modified by LU decomposition. Use gsl_A copy
    84         double* Acopy = xNew<double>(m*m);
    85         xMemCpy(Acopy,inVal,m*m);
    86         /*Initialize gsl matrices and vectors: */
    87         gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
    88         gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
    89         gsl_permutation *perm_p = gsl_permutation_alloc (m);
    90         int  signPerm;
    91         // factorize
    92         gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
    93         gsl_vector *gsl_x_p = gsl_vector_alloc (m);
    94         // solve for the value
    95         gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
    96         /*Copy result*/
    97         xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
    98         gsl_vector_free(gsl_x_p);
    99         //  for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx x["<< i << "]=" << outVal[i] << std::endl;
    100         // solve for the derivatives acc. to A * dx = r  with r=db - dA * x
    101         // compute the RHS
    102         double* r=xNew<double>(m);
    103         for (int i=0; i<m; i++) {
    104                 r[i]=inDeriv[m*m+i]; // this is db[i]
    105                 for (int j=0;j<m; j++) {
    106                         r[i]-=inDeriv[i*m+j]*outVal[j]; // this is dA[i][j]*x[j]
    107                 }
    108         }
    109         gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
    110         gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
    111         gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
    112         xMemCpy(outDeriv,gsl_vector_ptr(gsl_dx_p,0),m);
    113         gsl_vector_free(gsl_dx_p);
    114         xDelete(r);
    115         gsl_permutation_free(perm_p);
    116         xDelete(Acopy);
    117 #endif
    118         return 0;
    119 } /*}}}*/
    120 int EDF_fov_forward_for_solverx(int n, IssmPDouble *inVal, int directionCount, IssmPDouble **inDeriv, int m, IssmPDouble *outVal, IssmPDouble **outDeriv) { /*{{{*/
    121 #ifdef _HAVE_GSL_
    122         // the matrix will be modified by LU decomposition. Use gsl_A copy
    123         double* Acopy = xNew<double>(m*m);
    124         xMemCpy(Acopy,inVal,m*m);
    125         /*Initialize gsl matrices and vectors: */
    126         gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
    127         gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
    128         gsl_permutation *perm_p = gsl_permutation_alloc (m);
    129         int  signPerm;
    130         // factorize
    131         gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
    132         gsl_vector *gsl_x_p = gsl_vector_alloc (m);
    133         // solve for the value
    134         gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
    135         /*Copy result*/
    136         xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
    137         gsl_vector_free(gsl_x_p);
    138         // solve for the derivatives acc. to A * dx = r  with r=db - dA * x
    139         double* r=xNew<double>(m);
    140         gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
    141         for (int dir=0;dir<directionCount;++dir) {
    142                 // compute the RHS
    143                 for (int i=0; i<m; i++) {
    144                         r[i]=inDeriv[m*m+i][dir]; // this is db[i]
    145                         for (int j=0;j<m; j++) {
    146                                 r[i]-=inDeriv[i*m+j][dir]*outVal[j]; // this is dA[i][j]*x[j]
    147                         }
    148                 }
    149                 gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
    150                 gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
    151                 // reuse r
    152                 xMemCpy(r,gsl_vector_ptr(gsl_dx_p,0),m);
    153                 for (int i=0; i<m; i++) {
    154                         outDeriv[i][dir]=r[i];
    155                 }
    156         }
    157         gsl_vector_free(gsl_dx_p);
    158         xDelete(r);
    159         gsl_permutation_free(perm_p);
    160         xDelete(Acopy);
    161 #endif
    162         return 0;
    163 }
    164 /*}}}*/
    165 int EDF_fos_reverse_for_solverx(int m, double *dp_U, int n, double *dp_Z, double* dp_x, double* dp_y) { /*{{{*/
    166         // copy to transpose the matrix
    167         double* transposed=xNew<double>(m*m);
    168         for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j];
    169         gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m);
    170         // the adjoint of the solution is our right-hand side
    171         gsl_vector_view x_bar=gsl_vector_view_array(dp_U,m);
    172         // the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute:
    173         gsl_vector_view b_bar=gsl_vector_view_array(dp_Z+m*m,m);
    174         gsl_permutation *perm_p = gsl_permutation_alloc (m);
    175         int permSign;
    176         gsl_linalg_LU_decomp (&aTransposed.matrix, perm_p, &permSign);
    177         gsl_linalg_LU_solve (&aTransposed.matrix, perm_p, &x_bar.vector, &b_bar.vector);
    178         // now do the adjoint of the matrix
    179         for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) dp_Z[i*m+j]-=dp_Z[m*m+i]*dp_y[j];
    180         gsl_permutation_free(perm_p);
    181         xDelete(transposed);
    182         return 0;
    183 }
    184 /*}}}*/
    185 int EDF_fov_reverse_for_solverx(int m, int p, double **dpp_U, int n, double **dpp_Z, double* dp_x, double* dp_y) { /*{{{*/
    186         // copy to transpose the matrix
    187         double* transposed=xNew<double>(m*m);
    188         for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j];
    189         gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m);
    190         gsl_permutation *perm_p = gsl_permutation_alloc (m);
    191         int permSign;
    192         gsl_linalg_LU_decomp (&aTransposed.matrix, perm_p, &permSign);
    193         for (int weightsRowIndex=0;weightsRowIndex<p;++weightsRowIndex) {
    194                 // the adjoint of the solution is our right-hand side
    195                 gsl_vector_view x_bar=gsl_vector_view_array(dpp_U[weightsRowIndex],m);
    196                 // the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute:
    197                 gsl_vector_view b_bar=gsl_vector_view_array(dpp_Z[weightsRowIndex]+m*m,m);
    198                 gsl_linalg_LU_solve (&aTransposed.matrix, perm_p, &x_bar.vector, &b_bar.vector);
    199                 // now do the adjoint of the matrix
    200                 for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) dpp_Z[weightsRowIndex][i*m+j]-=dpp_Z[weightsRowIndex][m*m+i]*dp_y[j];
    201         }
    202         gsl_permutation_free(perm_p);
    203         xDelete(transposed);
    204         return 0;
    205 }
    206 /*}}}*/
    207 void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
    208         // pack inputs to conform to the EDF-prescribed interface
    209         // ensure a contiguous block of locations:
    210         ensureContiguousLocations(n*(n+1));
    211         IssmDouble*  adoubleEDFin=xNew<IssmDouble>(n*(n+1));  // packed inputs, i.e. matrix and right hand side
    212         for(int i=0; i<n*n;i++)adoubleEDFin[i]    =A[i];      // pack matrix
    213         for(int i=0; i<n;  i++)adoubleEDFin[i+n*n]=B[i];      // pack the right hand side
    214         IssmPDouble* pdoubleEDFin=xNew<IssmPDouble>(n*(n+1)); // provide space to transfer inputs during call_ext_fct
    215         IssmPDouble* pdoubleEDFout=xNew<IssmPDouble>(n);           // provide space to transfer outputs during call_ext_fct
    216         // call the wrapped solver through the registry entry we retrieve from parameters
    217         call_ext_fct(dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p,
    218                      n*(n+1), pdoubleEDFin, adoubleEDFin,
    219                      n, pdoubleEDFout,X);
    220         // for(int i=0; i<n;  i++) {ADOLC_DUMP_MACRO(X[i]);}
    221         xDelete(adoubleEDFin);
    222         xDelete(pdoubleEDFin);
    223         xDelete(pdoubleEDFout);
    224 }
    225 /*}}}*/
    226 #endif

  • issm/trunk-jpl/src/c/toolkits/issm/IssmSolver.h

    r14792 r14877  
    2323
    2424void IssmSolve(IssmVec<IssmDouble>** puf,IssmMat<IssmDouble>* Kff, IssmVec<IssmDouble>* pf,Parameters* parameters);
    25 void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n);
    26 void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n);
    27 
    28 #if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_)
    29 void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters);
    30 // call back functions:
    31 ADOLC_ext_fct EDF_for_solverx;
    32 ADOLC_ext_fct_fos_forward EDF_fos_forward_for_solverx;
    33 ADOLC_ext_fct_fos_reverse EDF_fos_reverse_for_solverx;
    34 ADOLC_ext_fct_fov_forward EDF_fov_forward_for_solverx;
    35 ADOLC_ext_fct_fov_reverse EDF_fov_reverse_for_solverx;
    36 #endif
    3725
    3826#endif

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

    r14633 r14877  
    2424#endif
    2525
     26#ifdef _HAVE_GSL_
     27#include "./gsl/gslincludes.h"
     28#endif
     29
    2630#include "./triangle/triangleincludes.h"
    2731#include "./toolkitsenums.h"
Note: See TracChangeset for help on using the changeset viewer.