source: issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp@ 13375

Last change on this file since 13375 was 13375, checked in by Mathieu Morlighem, 13 years ago

NEW: enable to Fetch matlab single vectors

File size: 6.8 KB
RevLine 
[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]25void 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]56void 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]70int 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]74int 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]112int 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]157int EDF_fos_reverse_for_solverx(int m, double *dp_U, int n, double *dp_Z) { /*{{{*/
158 return 0;
159}
160/*}}}*/
[13261]161void 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]178void 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/*}}}*/
Note: See TracBrowser for help on using the repository browser.