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

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

CHG: renaming _error2_ _error_
File size: 2.9 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
21void SolverxSeq(SeqVec** puf,SeqMat* Kff, SeqVec* pf){/*{{{*/
22
23 #ifdef _HAVE_GSL_
24 /*Intermediary: */
25 int M,N,N2,s;
26 SeqVec *uf = NULL;
27 IssmDouble *x = NULL;
28
29 Kff->GetSize(&M,&N);
30 pf->GetSize(&N2);
31
32 if(N!=N2)_error_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !");
33 if(M!=N)_error_("Stiffness matrix should be square!");
34
35 SolverxSeq(&x,Kff->matrix,pf->vector,N);
36 uf=new SeqVec(x,N);
37
38 /*Assign output pointers:*/
39 *puf=uf;
40
41 #else
42 _error_("GSL support not compiled in!");
43 #endif
44
45}/*}}}*/
46#ifdef _HAVE_ADOLC_
47void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n){//{{{
48 /* if we use Adol-C then the IssmDouble will be an adouble
49 and the calls to gsl_... will not work.
50 We therefore call a wrapped solver instead.
51 */
52
53 /*Output: */
54 IssmDouble* X=NULL;
55
56 /*Intermediary: */
57 int i;
58 double* doubleA=NULL;
59 double* doubleB=NULL;
60 double* doubleX=NULL;
61
62 /*First, transfer from IssmDouble to double all our matrices and vectors: */
63 doubleA=xNew<double>(n*n);
64 doubleB=xNew<double>(n);
65 for(i=0;i<n*n;i++)A[i]>>=doubleA[i];
66 for(i=0;i<n;i++)B[i]>>=doubleB[i];
67
68 /*Call wrapped solver: */
69 SolverxSeq(&doubleX,doubleA, doubleB, n);
70
71 /*Transfer solution vector from double to IssmDouble: */
72 X = xNew<IssmDouble>(n);
73 for(i=0;i<n;i++)X[i]<<=doubleX[i];
74
75 /*Free ressources:*/
76 xDelete<double>(doubleA);
77 xDelete<double>(doubleB);
78
79 /*Assign output pointers: */
80 *pX=X;
81}
82/*}}}*/
83#endif
84//void SolverxSeq(double** pX,double* A,double* B,int n){ //{{{
85#ifdef _HAVE_ADOLC_
86void SolverxSeq(double** pX,double* A,double* B,int n){
87#else
88void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n){
89#endif
90 #ifdef _HAVE_GSL_
91 /*GSL Matrices and vectors: */
92 int s;
93 gsl_matrix_view a;
94 gsl_vector_view b;
95 gsl_vector *x = NULL;
96 gsl_permutation *p = NULL;
97 /*A will be modified by LU decomposition. Use copy*/
98 double* Acopy = xNew<double>(n*n);
99 xMemCpy<double>(Acopy,A,n*n);
100
101 /*Initialize gsl matrices and vectors: */
102 a = gsl_matrix_view_array (Acopy,n,n);
103 b = gsl_vector_view_array (B,n);
104 x = gsl_vector_alloc (n);
105
106 /*Run LU and solve: */
107 p = gsl_permutation_alloc (n);
108 gsl_linalg_LU_decomp (&a.matrix, p, &s);
109 gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x);
110
111 //printf ("x = \n");
112 //gsl_vector_fprintf (stdout, x, "%g");
113
114 /*Copy result*/
115 double* X = xNew<double>(n);
116 memcpy(X,gsl_vector_ptr(x,0),n*sizeof(double));
117
118 /*Clean up and assign output pointer*/
119 xDelete<double>(Acopy);
120 gsl_permutation_free(p);
121 gsl_vector_free(x);
122 *pX=X;
123 #endif
124}
125/*}}}*/
Note: See TracBrowser for help on using the repository browser.