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

/*!\file SolverxGsl
 * \brief Gsl implementation of solver
 */

#include "./Solverx.h"
#include "../../shared/shared.h"
#include "../../include/include.h"
#include "../../io/io.h"

#ifdef HAVE_CONFIG_H
        #include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif
 
#include <gsl/gsl_linalg.h>

void    SolverxGsl(SeqVec** puf,SeqMat* Kff, SeqVec* pf){

        /*intermediary: */
        SeqMat* KffCopy=NULL;

        /*Output: */
        SeqVec*        uf               = NULL;

        
        /*GSL Matrices and vectors: */
        gsl_matrix_view m;
        gsl_vector_view b;
        gsl_vector* x=NULL;
        gsl_permutation* p=NULL;
        
        /*We are going to do an in place LU decomp, so we need to save the matrix with its original structure: */
        KffCopy=Kff->Duplicate();


        /*Intermediary: */
        int M,N,N2,s;

        Kff->GetSize(&M,&N);
        pf->GetSize(&N2);

        if(N!=N2)_error_("Right hand side vector of size %i, when matrix is of size %i-%i !",N2,M,N);
        if(M!=N)_error_("Stiffness matrix should be square!");

        /*Initialize gsl matrices and vectors: */
        m = gsl_matrix_view_array (KffCopy->matrix, M, N);
        b = gsl_vector_view_array (pf->vector, N);
        x = gsl_vector_alloc (N);

        /*Run LU and solve: */
        p = gsl_permutation_alloc (N);
        gsl_linalg_LU_decomp (&m.matrix, p, &s);
        gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);

        //printf ("x = \n");
        //gsl_vector_fprintf (stdout, x, "%g");

        /*Get uf initialized with the results: */
        uf=new SeqVec(gsl_vector_ptr(x,0),M);

        /*Free resources:*/
        gsl_permutation_free (p);
        gsl_vector_free (x);

        delete KffCopy;


        /*Assign output pointers:*/
        *puf=uf;
}