source: issm/trunk/src/c/modules/Solverx/SolverxGsl.cpp@ 11995

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

merged trunk-jpl and trunk for revision 11994M

File size: 1.5 KB
Line 
1/*!\file SolverxGsl
2 * \brief Gsl implementation of solver
3 */
4
5#include "./Solverx.h"
6#include "../../shared/shared.h"
7#include "../../include/include.h"
8#include "../../io/io.h"
9
10#ifdef HAVE_CONFIG_H
11 #include <config.h>
12#else
13#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
14#endif
15
16#include <gsl/gsl_linalg.h>
17
18void SolverxGsl(SeqVec** puf,SeqMat* Kff, SeqVec* pf){
19
20 /*intermediary: */
21 SeqMat* KffCopy=NULL;
22
23 /*Output: */
24 SeqVec* uf = NULL;
25
26
27 /*GSL Matrices and vectors: */
28 gsl_matrix_view m;
29 gsl_vector_view b;
30 gsl_vector* x=NULL;
31 gsl_permutation* p=NULL;
32
33 /*We are going to do an in place LU decomp, so we need to save the matrix with its original structure: */
34 KffCopy=Kff->Duplicate();
35
36
37 /*Intermediary: */
38 int M,N,N2,s;
39
40 Kff->GetSize(&M,&N);
41 pf->GetSize(&N2);
42
43 if(N!=N2)_error_("Right hand side vector of size %i, when matrix is of size %i-%i !",N2,M,N);
44 if(M!=N)_error_("Stiffness matrix should be square!");
45
46 /*Initialize gsl matrices and vectors: */
47 m = gsl_matrix_view_array (KffCopy->matrix, M, N);
48 b = gsl_vector_view_array (pf->vector, N);
49 x = gsl_vector_alloc (N);
50
51 /*Run LU and solve: */
52 p = gsl_permutation_alloc (N);
53 gsl_linalg_LU_decomp (&m.matrix, p, &s);
54 gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);
55
56 //printf ("x = \n");
57 //gsl_vector_fprintf (stdout, x, "%g");
58
59 /*Get uf initialized with the results: */
60 uf=new SeqVec(gsl_vector_ptr(x,0),M);
61
62 /*Free resources:*/
63 gsl_permutation_free (p);
64 gsl_vector_free (x);
65
66 delete KffCopy;
67
68
69 /*Assign output pointers:*/
70 *puf=uf;
71}
Note: See TracBrowser for help on using the repository browser.