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 |
|
---|
18 | void 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 | }
|
---|