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

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

BUG: fixed Kriging with new solver prototype
File size: 3.5 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 IssmDouble *x = xNew<IssmDouble>(N);
38#ifdef _HAVE_ADOLC_
39 SolverxSeq(x,Kff->matrix,pf->vector,N,parameters);
40#else
41 SolverxSeq(x,Kff->matrix,pf->vector,N);
42#endif
43 uf=new SeqVec<IssmDouble>(x,N);
44 xDelete(x);
45
46 /*Assign output pointers:*/
47 *puf=uf;
48
49 #else
50 _error_("GSL support not compiled in!");
51 #endif
52
53}/*}}}*/
54void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
55
56 /*Allocate output*/
57 double* X = xNew<double>(n);
58
59 /*Solve*/
60 SolverxSeq(X,A,B,n);
61
62 /*Assign output pointer*/
63 *pX=X;
64}
65/*}}}*/
66
67#ifdef _HAVE_ADOLC_
68int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y){/*{{{*/
69 if(m*(m+1)!=n)_error_("Stiffness matrix should be square!");
70 SolverxSeq(y,x, x+m*m, m);
71 return 0;
72}/*}}}*/
73void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
74 // pack inputs to conform to the EDF-prescribed interface
75 IssmDouble* adoubleEDFin=xNew<IssmDouble>(n*(n+1)); // packed inputs, i.e. matrix and right hand side
76 for(int i=0; i<n*n;i++)adoubleEDFin[i] =A[i]; // pack matrix
77 for(int i=0; i<n; i++)adoubleEDFin[i+n*n]=B[i]; // pack the right hand side
78 IssmPDouble* pdoubleEDFin=xNew<IssmPDouble>(n*(n+1)); // provide space to transfer inputs during call_ext_fct
79 IssmPDouble* pdoubleEDFout=xNew<IssmPDouble>(n); // provide space to transfer outputs during call_ext_fct
80 // call the wrapped solver through the registry entry we retrieve from parameters
81 call_ext_fct(dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p,
82 n*(n+1), pdoubleEDFin, adoubleEDFin,
83 n, pdoubleEDFout,X);
84 xDelete(adoubleEDFin);
85 xDelete(pdoubleEDFin);
86 xDelete(pdoubleEDFout);
87}
88/*}}}*/
89#endif
90void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
91#ifdef _HAVE_GSL_
92 /*GSL Matrices and vectors: */
93 int s;
94 gsl_matrix_view a;
95 gsl_vector_view b;
96 gsl_vector *x = NULL;
97 gsl_permutation *p = NULL;
98 /*A will be modified by LU decomposition. Use copy*/
99 double* Acopy = xNew<double>(n*n);
100 xMemCpy(Acopy,A,n*n);
101
102 /*Initialize gsl matrices and vectors: */
103 a = gsl_matrix_view_array (Acopy,n,n);
104 b = gsl_vector_view_array (B,n);
105 x = gsl_vector_alloc (n);
106
107 /*Run LU and solve: */
108 p = gsl_permutation_alloc (n);
109 gsl_linalg_LU_decomp (&a.matrix, p, &s);
110 gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x);
111
112 //printf ("x = \n");
113 //gsl_vector_fprintf (stdout, x, "%g");
114
115 /*Copy result*/
116 xMemCpy(X,gsl_vector_ptr(x,0),n);
117
118 /*Clean up and assign output pointer*/
119 xDelete(Acopy);
120 gsl_permutation_free(p);
121 gsl_vector_free(x);
122#endif
123}
124/*}}}*/
Note: See TracBrowser for help on using the repository browser.