Ignore:
Timestamp:
11/16/12 08:10:16 (12 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 13974

Location:
issm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/c/modules/Solverx/SolverxSeq.cpp

    r13395 r13975  
    1212#include "./Solverx.h"
    1313#include "../../shared/shared.h"
     14#include "../../classes/classes.h"
    1415#include "../../include/include.h"
    1516#include "../../io/io.h"
    … …  
    1920#endif
    2021
    21 #ifdef _HAVE_ADOLC_
    22 #include "../../shared/Numerics/adolc_edf.h"
    23 #endif
    24 
    2522void SolverxSeq(SeqVec<IssmDouble>** puf,SeqMat<IssmDouble>* Kff, SeqVec<IssmDouble>* pf, Parameters* parameters){/*{{{*/
    2623
    2724#ifdef _HAVE_GSL_
    2825        /*Intermediary: */
    29         int M,N,N2,s;
     26        int M,N,N2;
    3027        SeqVec<IssmDouble> *uf = NULL;
    3128
    … …  
    3532        if(N!=N2)_error_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !");
    3633        if(M!=N)_error_("Stiffness matrix should be square!");
     34#ifdef _HAVE_ADOLC_
     35        ensureContiguousLocations(N);
     36#endif
    3737        IssmDouble *x  = xNew<IssmDouble>(N);
    38 
    3938#ifdef _HAVE_ADOLC_
    4039        SolverxSeq(x,Kff->matrix,pf->vector,N,parameters);
    … …  
    6968#ifdef _HAVE_ADOLC_
    7069int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y){ /*{{{*/
    71     SolverxSeq(y,x, x+m*m, m); // x is where the matrix starts, x+m*m is where the right-hand side starts
    72     return 0;
     70        SolverxSeq(y,x, x+m*m, m); // x is where the matrix starts, x+m*m is where the right-hand side starts
     71        return 0;
    7372} /*}}}*/
    7473int EDF_fos_forward_for_solverx(int n, IssmPDouble *inVal, IssmPDouble *inDeriv, int m, IssmPDouble *outVal, IssmPDouble *outDeriv) { /*{{{*/
    7574#ifdef _HAVE_GSL_
    76   // the matrix will be modified by LU decomposition. Use gsl_A copy
    77   double* Acopy = xNew<double>(m*m);
    78   xMemCpy(Acopy,inVal,m*m);
    79   /*Initialize gsl matrices and vectors: */
    80   gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
    81   gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
    82   gsl_permutation *perm_p = gsl_permutation_alloc (m);
    83   int  signPerm;
    84   // factorize
    85   gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
    86   gsl_vector *gsl_x_p = gsl_vector_alloc (m);
    87   // solve for the value
    88   gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
    89   /*Copy result*/
    90   xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
    91   gsl_vector_free(gsl_x_p);
    92   // solve for the derivatives acc. to A * dx = r  with r=db - dA * x
    93   // compute the RHS
    94   double* r=xNew<double>(m);
    95   for (int i=0; i<m; i++) {
    96     r[i]=inDeriv[m*m+i]; // this is db[i]
    97     for (int j=0;j<m; j++) {
    98       r[i]-=inDeriv[i*m+j]*outVal[j]; // this is dA[i][j]*x[j]
    99     }
    100   }
    101   gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
    102   gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
    103   gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
    104   xMemCpy(outDeriv,gsl_vector_ptr(gsl_dx_p,0),m);
    105   gsl_vector_free(gsl_dx_p);
    106   xDelete(r);
    107   gsl_permutation_free(perm_p);
    108   xDelete(Acopy);
    109   #endif
    110   return 0;
     75        //  for (int i=0; i<m*m; ++i) std::cout << "EDF_fos_forward_for_solverx A["<< i << "]=" << inVal[i] << std::endl;
     76        //  for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx b["<< i << "]=" << inVal[i+m*m] << std::endl;
     77        // the matrix will be modified by LU decomposition. Use gsl_A copy
     78        double* Acopy = xNew<double>(m*m);
     79        xMemCpy(Acopy,inVal,m*m);
     80        /*Initialize gsl matrices and vectors: */
     81        gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
     82        gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
     83        gsl_permutation *perm_p = gsl_permutation_alloc (m);
     84        int  signPerm;
     85        // factorize
     86        gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
     87        gsl_vector *gsl_x_p = gsl_vector_alloc (m);
     88        // solve for the value
     89        gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
     90        /*Copy result*/
     91        xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
     92        gsl_vector_free(gsl_x_p);
     93        //  for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx x["<< i << "]=" << outVal[i] << std::endl;
     94        // solve for the derivatives acc. to A * dx = r  with r=db - dA * x
     95        // compute the RHS
     96        double* r=xNew<double>(m);
     97        for (int i=0; i<m; i++) {
     98                r[i]=inDeriv[m*m+i]; // this is db[i]
     99                for (int j=0;j<m; j++) {
     100                        r[i]-=inDeriv[i*m+j]*outVal[j]; // this is dA[i][j]*x[j]
     101                }
     102        }
     103        gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
     104        gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
     105        gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
     106        xMemCpy(outDeriv,gsl_vector_ptr(gsl_dx_p,0),m);
     107        gsl_vector_free(gsl_dx_p);
     108        xDelete(r);
     109        gsl_permutation_free(perm_p);
     110        xDelete(Acopy);
     111#endif
     112        return 0;
    111113} /*}}}*/
    112114int EDF_fov_forward_for_solverx(int n, IssmPDouble *inVal, int directionCount, IssmPDouble **inDeriv, int m, IssmPDouble *outVal, IssmPDouble **outDeriv) { /*{{{*/
    113115#ifdef _HAVE_GSL_
    114   // the matrix will be modified by LU decomposition. Use gsl_A copy
    115   double* Acopy = xNew<double>(m*m);
    116   xMemCpy(Acopy,inVal,m*m);
    117   /*Initialize gsl matrices and vectors: */
    118   gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
    119   gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
    120   gsl_permutation *perm_p = gsl_permutation_alloc (m);
    121   int  signPerm;
    122   // factorize
    123   gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
    124   gsl_vector *gsl_x_p = gsl_vector_alloc (m);
    125   // solve for the value
    126   gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
    127   /*Copy result*/
    128   xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
    129   gsl_vector_free(gsl_x_p);
    130   // solve for the derivatives acc. to A * dx = r  with r=db - dA * x
    131   double* r=xNew<double>(m);
    132   gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
    133   for (int dir=0;dir<directionCount;++dir) {
    134     // compute the RHS
    135     for (int i=0; i<m; i++) {
    136       r[i]=inDeriv[m*m+i][dir]; // this is db[i]
    137       for (int j=0;j<m; j++) {
    138         r[i]-=inDeriv[i*m+j][dir]*outVal[j]; // this is dA[i][j]*x[j]
    139       }
    140     }
    141     gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
    142     gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
    143     // reuse r
    144     xMemCpy(r,gsl_vector_ptr(gsl_dx_p,0),m);
    145     for (int i=0; i<m; i++) {
    146       outDeriv[i][dir]=r[i];
    147     }
    148   }
    149   gsl_vector_free(gsl_dx_p);
    150   xDelete(r);
    151   gsl_permutation_free(perm_p);
    152   xDelete(Acopy);
    153   #endif
    154   return 0;
    155 }
    156 /*}}}*/
    157 int EDF_fos_reverse_for_solverx(int m, double *dp_U, int n, double *dp_Z) { /*{{{*/
    158   return 0;
     116        // the matrix will be modified by LU decomposition. Use gsl_A copy
     117        double* Acopy = xNew<double>(m*m);
     118        xMemCpy(Acopy,inVal,m*m);
     119        /*Initialize gsl matrices and vectors: */
     120        gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
     121        gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
     122        gsl_permutation *perm_p = gsl_permutation_alloc (m);
     123        int  signPerm;
     124        // factorize
     125        gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
     126        gsl_vector *gsl_x_p = gsl_vector_alloc (m);
     127        // solve for the value
     128        gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
     129        /*Copy result*/
     130        xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
     131        gsl_vector_free(gsl_x_p);
     132        // solve for the derivatives acc. to A * dx = r  with r=db - dA * x
     133        double* r=xNew<double>(m);
     134        gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
     135        for (int dir=0;dir<directionCount;++dir) {
     136                // compute the RHS
     137                for (int i=0; i<m; i++) {
     138                        r[i]=inDeriv[m*m+i][dir]; // this is db[i]
     139                        for (int j=0;j<m; j++) {
     140                                r[i]-=inDeriv[i*m+j][dir]*outVal[j]; // this is dA[i][j]*x[j]
     141                        }
     142                }
     143                gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
     144                gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
     145                // reuse r
     146                xMemCpy(r,gsl_vector_ptr(gsl_dx_p,0),m);
     147                for (int i=0; i<m; i++) {
     148                        outDeriv[i][dir]=r[i];
     149                }
     150        }
     151        gsl_vector_free(gsl_dx_p);
     152        xDelete(r);
     153        gsl_permutation_free(perm_p);
     154        xDelete(Acopy);
     155#endif
     156        return 0;
     157}
     158/*}}}*/
     159int EDF_fos_reverse_for_solverx(int m, double *dp_U, int n, double *dp_Z, double* dp_x, double* dp_y) { /*{{{*/
     160        // copy to transpose the matrix
     161        double* transposed=xNew<double>(m*m);
     162        for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j];
     163        gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m);
     164        // the adjoint of the solution is our right-hand side
     165        gsl_vector_view x_bar=gsl_vector_view_array(dp_U,m);
     166        // the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute:
     167        gsl_vector_view b_bar=gsl_vector_view_array(dp_Z+m*m,m);
     168        gsl_permutation *perm_p = gsl_permutation_alloc (m);
     169        int permSign;
     170        gsl_linalg_LU_decomp (&aTransposed.matrix, perm_p, &permSign);
     171        gsl_linalg_LU_solve (&aTransposed.matrix, perm_p, &x_bar.vector, &b_bar.vector);
     172        // now do the adjoint of the matrix
     173        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];
     174        gsl_permutation_free(perm_p);
     175        xDelete(transposed);
     176        return 0;
     177}
     178/*}}}*/
     179int EDF_fov_reverse_for_solverx(int m, int p, double **dpp_U, int n, double **dpp_Z, double* dp_x, double* dp_y) { /*{{{*/
     180        // copy to transpose the matrix
     181        double* transposed=xNew<double>(m*m);
     182        for (int i=0; i<m; ++i) for (int j=0; j<m; ++j) transposed[j*m+i]=dp_x[i*m+j];
     183        gsl_matrix_view aTransposed = gsl_matrix_view_array (transposed,m,m);
     184        gsl_permutation *perm_p = gsl_permutation_alloc (m);
     185        int permSign;
     186        gsl_linalg_LU_decomp (&aTransposed.matrix, perm_p, &permSign);
     187        for (int weightsRowIndex=0;weightsRowIndex<p;++weightsRowIndex) {
     188                // the adjoint of the solution is our right-hand side
     189                gsl_vector_view x_bar=gsl_vector_view_array(dpp_U[weightsRowIndex],m);
     190                // the last m elements of dp_Z representing the adjoint of the right-hand side we want to compute:
     191                gsl_vector_view b_bar=gsl_vector_view_array(dpp_Z[weightsRowIndex]+m*m,m);
     192                gsl_linalg_LU_solve (&aTransposed.matrix, perm_p, &x_bar.vector, &b_bar.vector);
     193                // now do the adjoint of the matrix
     194                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];
     195        }
     196        gsl_permutation_free(perm_p);
     197        xDelete(transposed);
     198        return 0;
    159199}
    160200/*}}}*/
    161201void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
    162202        // pack inputs to conform to the EDF-prescribed interface
     203        // ensure a contiguous block of locations:
     204        ensureContiguousLocations(n*(n+1));
    163205        IssmDouble*  adoubleEDFin=xNew<IssmDouble>(n*(n+1));  // packed inputs, i.e. matrix and right hand side
    164206        for(int i=0; i<n*n;i++)adoubleEDFin[i]    =A[i];      // pack matrix
    … …  
    170212                     n*(n+1), pdoubleEDFin, adoubleEDFin,
    171213                     n, pdoubleEDFout,X);
     214        // for(int i=0; i<n;  i++) {ADOLC_DUMP_MACRO(X[i]);}
    172215        xDelete(adoubleEDFin);
    173216        xDelete(pdoubleEDFin);
    … …  
    181224        int              s;
    182225        gsl_matrix_view  a;
    183         gsl_vector_view  b;
    184         gsl_vector      *x = NULL;
     226        gsl_vector_view  b,x;
    185227        gsl_permutation *p = NULL;
     228//      for (int i=0; i<n*n; ++i) std::cout << "SolverxSeq A["<< i << "]=" << A[i] << std::endl;
     229//      for (int i=0; i<n; ++i) std::cout << "SolverxSeq b["<< i << "]=" << B[i] << std::endl;
    186230        /*A will be modified by LU decomposition. Use copy*/
    187231        double* Acopy = xNew<double>(n*n);
    … …  
    191235        a = gsl_matrix_view_array (Acopy,n,n);
    192236        b = gsl_vector_view_array (B,n);
    193         x = gsl_vector_alloc (n);
     237        x = gsl_vector_view_array (X,n);
    194238
    195239        /*Run LU and solve: */
    196240        p = gsl_permutation_alloc (n);
    197241        gsl_linalg_LU_decomp (&a.matrix, p, &s);
    198         gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x);
    199 
    200         //printf ("x = \n");
    201         //gsl_vector_fprintf (stdout, x, "%g");
    202 
    203         /*Copy result*/
    204         xMemCpy(X,gsl_vector_ptr(x,0),n);
     242        gsl_linalg_LU_solve (&a.matrix, p, &b.vector, &x.vector);
    205243
    206244        /*Clean up and assign output pointer*/
    207245        xDelete(Acopy);
    208246        gsl_permutation_free(p);
    209         gsl_vector_free(x);
    210 #endif
    211 }
    212 /*}}}*/
     247#endif
     248}
     249/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.