[12850] | 1 | /*!\file SolverxSeq
|
---|
| 2 | * \brief implementation of sequential solver using the GSL librarie
|
---|
[11726] | 3 | */
|
---|
| 4 |
|
---|
[12445] | 5 | #ifdef HAVE_CONFIG_H
|
---|
| 6 | #include <config.h>
|
---|
| 7 | #else
|
---|
| 8 | #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
|
---|
| 9 | #endif
|
---|
[12444] | 10 | #include <cstring>
|
---|
| 11 |
|
---|
[11726] | 12 | #include "./Solverx.h"
|
---|
| 13 | #include "../../shared/shared.h"
|
---|
| 14 | #include "../../include/include.h"
|
---|
| 15 | #include "../../io/io.h"
|
---|
| 16 |
|
---|
[12850] | 17 | #ifdef _HAVE_GSL_
|
---|
| 18 | #include <gsl/gsl_linalg.h>
|
---|
| 19 | #endif
|
---|
[11726] | 20 |
|
---|
[13195] | 21 | #ifdef _HAVE_ADOLC_
|
---|
| 22 | #include "../../shared/Numerics/adolc_edf.h"
|
---|
| 23 | #endif
|
---|
[12850] | 24 |
|
---|
[13216] | 25 | void SolverxSeq(SeqVec<IssmDouble>** puf,SeqMat<IssmDouble>* Kff, SeqVec<IssmDouble>* pf, Parameters* parameters){/*{{{*/
|
---|
[13195] | 26 |
|
---|
[13375] | 27 | #ifdef _HAVE_GSL_
|
---|
[11726] | 28 | /*Intermediary: */
|
---|
| 29 | int M,N,N2,s;
|
---|
[13216] | 30 | SeqVec<IssmDouble> *uf = NULL;
|
---|
[11726] | 31 |
|
---|
| 32 | Kff->GetSize(&M,&N);
|
---|
| 33 | pf->GetSize(&N2);
|
---|
| 34 |
|
---|
[13056] | 35 | if(N!=N2)_error_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !");
|
---|
| 36 | if(M!=N)_error_("Stiffness matrix should be square!");
|
---|
[13375] | 37 | IssmDouble *x = xNew<IssmDouble>(N);
|
---|
| 38 |
|
---|
[13195] | 39 | #ifdef _HAVE_ADOLC_
|
---|
[13196] | 40 | SolverxSeq(x,Kff->matrix,pf->vector,N,parameters);
|
---|
[13195] | 41 | #else
|
---|
[13196] | 42 | SolverxSeq(x,Kff->matrix,pf->vector,N);
|
---|
[13195] | 43 | #endif
|
---|
[13375] | 44 |
|
---|
[13216] | 45 | uf=new SeqVec<IssmDouble>(x,N);
|
---|
[13196] | 46 | xDelete(x);
|
---|
[11726] | 47 |
|
---|
| 48 | /*Assign output pointers:*/
|
---|
| 49 | *puf=uf;
|
---|
[12850] | 50 |
|
---|
[13375] | 51 | #else
|
---|
| 52 | _error_("GSL support not compiled in!");
|
---|
| 53 | #endif
|
---|
[12850] | 54 |
|
---|
[12417] | 55 | }/*}}}*/
|
---|
[13261] | 56 | void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
|
---|
| 57 |
|
---|
| 58 | /*Allocate output*/
|
---|
| 59 | double* X = xNew<double>(n);
|
---|
| 60 |
|
---|
| 61 | /*Solve*/
|
---|
| 62 | SolverxSeq(X,A,B,n);
|
---|
| 63 |
|
---|
| 64 | /*Assign output pointer*/
|
---|
| 65 | *pX=X;
|
---|
| 66 | }
|
---|
| 67 | /*}}}*/
|
---|
| 68 |
|
---|
[12988] | 69 | #ifdef _HAVE_ADOLC_
|
---|
[13295] | 70 | int 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
|
---|
[13185] | 72 | return 0;
|
---|
[13334] | 73 | } /*}}}*/
|
---|
[13295] | 74 | int EDF_fos_forward_for_solverx(int n, IssmPDouble *inVal, IssmPDouble *inDeriv, int m, IssmPDouble *outVal, IssmPDouble *outDeriv) { /*{{{*/
|
---|
| 75 | #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++) {
|
---|
[13364] | 98 | r[i]-=inDeriv[i*m+j]*outVal[j]; // this is dA[i][j]*x[j]
|
---|
[13295] | 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;
|
---|
[13334] | 111 | } /*}}}*/
|
---|
[13297] | 112 | int EDF_fov_forward_for_solverx(int n, IssmPDouble *inVal, int directionCount, IssmPDouble **inDeriv, int m, IssmPDouble *outVal, IssmPDouble **outDeriv) { /*{{{*/
|
---|
| 113 | #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++) {
|
---|
[13364] | 138 | r[i]-=inDeriv[i*m+j][dir]*outVal[j]; // this is dA[i][j]*x[j]
|
---|
[13297] | 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 | /*}}}*/
|
---|
[13313] | 157 | int EDF_fos_reverse_for_solverx(int m, double *dp_U, int n, double *dp_Z) { /*{{{*/
|
---|
| 158 | return 0;
|
---|
| 159 | }
|
---|
| 160 | /*}}}*/
|
---|
[13261] | 161 | void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
|
---|
[13195] | 162 | // pack inputs to conform to the EDF-prescribed interface
|
---|
[13261] | 163 | IssmDouble* adoubleEDFin=xNew<IssmDouble>(n*(n+1)); // packed inputs, i.e. matrix and right hand side
|
---|
| 164 | for(int i=0; i<n*n;i++)adoubleEDFin[i] =A[i]; // pack matrix
|
---|
| 165 | for(int i=0; i<n; i++)adoubleEDFin[i+n*n]=B[i]; // pack the right hand side
|
---|
[13196] | 166 | IssmPDouble* pdoubleEDFin=xNew<IssmPDouble>(n*(n+1)); // provide space to transfer inputs during call_ext_fct
|
---|
[13261] | 167 | IssmPDouble* pdoubleEDFout=xNew<IssmPDouble>(n); // provide space to transfer outputs during call_ext_fct
|
---|
[13195] | 168 | // call the wrapped solver through the registry entry we retrieve from parameters
|
---|
| 169 | call_ext_fct(dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p,
|
---|
[13196] | 170 | n*(n+1), pdoubleEDFin, adoubleEDFin,
|
---|
| 171 | n, pdoubleEDFout,X);
|
---|
| 172 | xDelete(adoubleEDFin);
|
---|
| 173 | xDelete(pdoubleEDFin);
|
---|
| 174 | xDelete(pdoubleEDFout);
|
---|
[12988] | 175 | }
|
---|
| 176 | /*}}}*/
|
---|
| 177 | #endif
|
---|
[13261] | 178 | void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
|
---|
| 179 | #ifdef _HAVE_GSL_
|
---|
[12988] | 180 | /*GSL Matrices and vectors: */
|
---|
| 181 | int s;
|
---|
| 182 | gsl_matrix_view a;
|
---|
| 183 | gsl_vector_view b;
|
---|
| 184 | gsl_vector *x = NULL;
|
---|
| 185 | gsl_permutation *p = NULL;
|
---|
| 186 | /*A will be modified by LU decomposition. Use copy*/
|
---|
| 187 | double* Acopy = xNew<double>(n*n);
|
---|
[13196] | 188 | xMemCpy(Acopy,A,n*n);
|
---|
[12417] | 189 |
|
---|
[12988] | 190 | /*Initialize gsl matrices and vectors: */
|
---|
| 191 | a = gsl_matrix_view_array (Acopy,n,n);
|
---|
| 192 | b = gsl_vector_view_array (B,n);
|
---|
| 193 | x = gsl_vector_alloc (n);
|
---|
[12417] | 194 |
|
---|
[12988] | 195 | /*Run LU and solve: */
|
---|
| 196 | p = gsl_permutation_alloc (n);
|
---|
| 197 | gsl_linalg_LU_decomp (&a.matrix, p, &s);
|
---|
| 198 | gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x);
|
---|
[12417] | 199 |
|
---|
[12988] | 200 | //printf ("x = \n");
|
---|
| 201 | //gsl_vector_fprintf (stdout, x, "%g");
|
---|
[12417] | 202 |
|
---|
[12988] | 203 | /*Copy result*/
|
---|
[13196] | 204 | xMemCpy(X,gsl_vector_ptr(x,0),n);
|
---|
[12417] | 205 |
|
---|
[12988] | 206 | /*Clean up and assign output pointer*/
|
---|
[13196] | 207 | xDelete(Acopy);
|
---|
[12988] | 208 | gsl_permutation_free(p);
|
---|
| 209 | gsl_vector_free(x);
|
---|
[13261] | 210 | #endif
|
---|
[12988] | 211 | }
|
---|
| 212 | /*}}}*/
|
---|