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

Last change on this file since 13407 was 13407, checked in by utke, 12 years ago

CHG use call to ensure contiguous locations for arguments passed to EDF ; commented out debug prints - to be removed when regression is successfull

File size: 7.6 KB
Line 
1/*!\file SolverxSeq
2 * \brief implementation of sequential solver using the GSL librarie
3 */
4
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
10#include <cstring>
11
12#include "./Solverx.h"
13#include "../../shared/shared.h"
14#include "../../include/include.h"
15#include "../../io/io.h"
16
17#ifdef _HAVE_GSL_
18#include <gsl/gsl_linalg.h>
19#endif
20
21#ifdef _HAVE_ADOLC_
22#include "../../shared/Numerics/adolc_edf.h"
23#endif
24
25void SolverxSeq(SeqVec<IssmDouble>** puf,SeqMat<IssmDouble>* Kff, SeqVec<IssmDouble>* pf, Parameters* parameters){/*{{{*/
26
27#ifdef _HAVE_GSL_
28 /*Intermediary: */
29 int M,N,N2,s;
30 SeqVec<IssmDouble> *uf = NULL;
31
32 Kff->GetSize(&M,&N);
33 pf->GetSize(&N2);
34
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!");
37#ifdef _HAVE_ADOLC_
38 ensureContiguousLocations(N);
39#endif
40 IssmDouble *x = xNew<IssmDouble>(N);
41#ifdef _HAVE_ADOLC_
42 SolverxSeq(x,Kff->matrix,pf->vector,N,parameters);
43#else
44 SolverxSeq(x,Kff->matrix,pf->vector,N);
45#endif
46
47 uf=new SeqVec<IssmDouble>(x,N);
48 xDelete(x);
49
50 /*Assign output pointers:*/
51 *puf=uf;
52
53#else
54 _error_("GSL support not compiled in!");
55#endif
56
57}/*}}}*/
58void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
59
60 /*Allocate output*/
61 double* X = xNew<double>(n);
62
63 /*Solve*/
64 SolverxSeq(X,A,B,n);
65
66 /*Assign output pointer*/
67 *pX=X;
68}
69/*}}}*/
70
71#ifdef _HAVE_ADOLC_
72int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y){ /*{{{*/
73 SolverxSeq(y,x, x+m*m, m); // x is where the matrix starts, x+m*m is where the right-hand side starts
74 return 0;
75} /*}}}*/
76int EDF_fos_forward_for_solverx(int n, IssmPDouble *inVal, IssmPDouble *inDeriv, int m, IssmPDouble *outVal, IssmPDouble *outDeriv) { /*{{{*/
77#ifdef _HAVE_GSL_
78// for (int i=0; i<m*m; ++i) std::cout << "EDF_fos_forward_for_solverx A["<< i << "]=" << inVal[i] << std::endl;
79// for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx b["<< i << "]=" << inVal[i+m*m] << std::endl;
80 // the matrix will be modified by LU decomposition. Use gsl_A copy
81 double* Acopy = xNew<double>(m*m);
82 xMemCpy(Acopy,inVal,m*m);
83 /*Initialize gsl matrices and vectors: */
84 gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
85 gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
86 gsl_permutation *perm_p = gsl_permutation_alloc (m);
87 int signPerm;
88 // factorize
89 gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
90 gsl_vector *gsl_x_p = gsl_vector_alloc (m);
91 // solve for the value
92 gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
93 /*Copy result*/
94 xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
95 gsl_vector_free(gsl_x_p);
96// for (int i=0; i<m; ++i) std::cout << "EDF_fos_forward_for_solverx x["<< i << "]=" << outVal[i] << std::endl;
97 // solve for the derivatives acc. to A * dx = r with r=db - dA * x
98 // compute the RHS
99 double* r=xNew<double>(m);
100 for (int i=0; i<m; i++) {
101 r[i]=inDeriv[m*m+i]; // this is db[i]
102 for (int j=0;j<m; j++) {
103 r[i]-=inDeriv[i*m+j]*outVal[j]; // this is dA[i][j]*x[j]
104 }
105 }
106 gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
107 gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
108 gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
109 xMemCpy(outDeriv,gsl_vector_ptr(gsl_dx_p,0),m);
110 gsl_vector_free(gsl_dx_p);
111 xDelete(r);
112 gsl_permutation_free(perm_p);
113 xDelete(Acopy);
114 #endif
115 return 0;
116} /*}}}*/
117int EDF_fov_forward_for_solverx(int n, IssmPDouble *inVal, int directionCount, IssmPDouble **inDeriv, int m, IssmPDouble *outVal, IssmPDouble **outDeriv) { /*{{{*/
118#ifdef _HAVE_GSL_
119 // the matrix will be modified by LU decomposition. Use gsl_A copy
120 double* Acopy = xNew<double>(m*m);
121 xMemCpy(Acopy,inVal,m*m);
122 /*Initialize gsl matrices and vectors: */
123 gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
124 gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
125 gsl_permutation *perm_p = gsl_permutation_alloc (m);
126 int signPerm;
127 // factorize
128 gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
129 gsl_vector *gsl_x_p = gsl_vector_alloc (m);
130 // solve for the value
131 gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
132 /*Copy result*/
133 xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
134 gsl_vector_free(gsl_x_p);
135 // solve for the derivatives acc. to A * dx = r with r=db - dA * x
136 double* r=xNew<double>(m);
137 gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
138 for (int dir=0;dir<directionCount;++dir) {
139 // compute the RHS
140 for (int i=0; i<m; i++) {
141 r[i]=inDeriv[m*m+i][dir]; // this is db[i]
142 for (int j=0;j<m; j++) {
143 r[i]-=inDeriv[i*m+j][dir]*outVal[j]; // this is dA[i][j]*x[j]
144 }
145 }
146 gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
147 gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
148 // reuse r
149 xMemCpy(r,gsl_vector_ptr(gsl_dx_p,0),m);
150 for (int i=0; i<m; i++) {
151 outDeriv[i][dir]=r[i];
152 }
153 }
154 gsl_vector_free(gsl_dx_p);
155 xDelete(r);
156 gsl_permutation_free(perm_p);
157 xDelete(Acopy);
158 #endif
159 return 0;
160}
161/*}}}*/
162int EDF_fos_reverse_for_solverx(int m, double *dp_U, int n, double *dp_Z) { /*{{{*/
163 return 0;
164}
165/*}}}*/
166void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
167 // pack inputs to conform to the EDF-prescribed interface
168 // ensure a contiguous block of locations:
169 ensureContiguousLocations(n*(n+1));
170 IssmDouble* adoubleEDFin=xNew<IssmDouble>(n*(n+1)); // packed inputs, i.e. matrix and right hand side
171 for(int i=0; i<n*n;i++)adoubleEDFin[i] =A[i]; // pack matrix
172 for(int i=0; i<n; i++)adoubleEDFin[i+n*n]=B[i]; // pack the right hand side
173 IssmPDouble* pdoubleEDFin=xNew<IssmPDouble>(n*(n+1)); // provide space to transfer inputs during call_ext_fct
174 IssmPDouble* pdoubleEDFout=xNew<IssmPDouble>(n); // provide space to transfer outputs during call_ext_fct
175 // call the wrapped solver through the registry entry we retrieve from parameters
176 call_ext_fct(dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p,
177 n*(n+1), pdoubleEDFin, adoubleEDFin,
178 n, pdoubleEDFout,X);
179 // for(int i=0; i<n; i++) {ADOLC_DUMP_MACRO(X[i]);}
180 xDelete(adoubleEDFin);
181 xDelete(pdoubleEDFin);
182 xDelete(pdoubleEDFout);
183}
184/*}}}*/
185#endif
186void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
187#ifdef _HAVE_GSL_
188 /*GSL Matrices and vectors: */
189 int s;
190 gsl_matrix_view a;
191 gsl_vector_view b;
192 gsl_vector *x = NULL;
193 gsl_permutation *p = NULL;
194// for (int i=0; i<n*n; ++i) std::cout << "SolverxSeq A["<< i << "]=" << A[i] << std::endl;
195// for (int i=0; i<n; ++i) std::cout << "SolverxSeq b["<< i << "]=" << B[i] << std::endl;
196 /*A will be modified by LU decomposition. Use copy*/
197 double* Acopy = xNew<double>(n*n);
198 xMemCpy(Acopy,A,n*n);
199
200 /*Initialize gsl matrices and vectors: */
201 a = gsl_matrix_view_array (Acopy,n,n);
202 b = gsl_vector_view_array (B,n);
203 x = gsl_vector_alloc (n);
204
205 /*Run LU and solve: */
206 p = gsl_permutation_alloc (n);
207 gsl_linalg_LU_decomp (&a.matrix, p, &s);
208 gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x);
209
210 //printf ("x = \n");
211 //gsl_vector_fprintf (stdout, x, "%g");
212
213 /*Copy result*/
214 xMemCpy(X,gsl_vector_ptr(x,0),n);
215// for (int i=0; i<n; ++i) std::cout << "SolverxSeq x["<< i << "]=" << X[i] << std::endl;
216
217 /*Clean up and assign output pointer*/
218 xDelete(Acopy);
219 gsl_permutation_free(p);
220 gsl_vector_free(x);
221#endif
222}
223/*}}}*/
Note: See TracBrowser for help on using the repository browser.