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

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

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

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

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

        /*Intermediary: */
        int M,N,N2,s;
        SeqVec *uf = NULL;
        IssmDouble *x  = NULL;

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

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

        SolverxGsl(&x,Kff->matrix,pf->vector,N);
        uf=new SeqVec(x,N);

        /*Assign output pointers:*/
        *puf=uf;
}/*}}}*/
void SolverxGsl(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n){/*{{{*/

        /*GSL Matrices and vectors: */
        int              s;
        gsl_matrix_view  a;
        gsl_vector_view  b;
        gsl_vector      *x = NULL;
        gsl_permutation *p = NULL;
#ifdef _HAVE_ADOLC_
        // if we use Adol-C then the IssmDouble will be an adouble
        // and the calls to gsl_... will not work
        // and we should call  a suitable wrapped solve instead
        _error2_("SolverxGsl: should not be here with Adol-C");
#else
        /*A will be modified by LU decomposition. Use copy*/
        IssmDouble* Acopy = xNew<IssmDouble>(n*n);
        xMemCpy<IssmDouble>(Acopy,A,n*n);

        /*Initialize gsl matrices and vectors: */
        a = gsl_matrix_view_array (Acopy,n,n);
        b = gsl_vector_view_array (B,n);
        x = gsl_vector_alloc (n);

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

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

        /*Copy result*/
        IssmDouble* X = xNew<IssmDouble>(n);
        memcpy(X,gsl_vector_ptr(x,0),n*sizeof(IssmDouble));

        /*Clean up and assign output pointer*/
        xDelete<IssmDouble>(Acopy);
        gsl_permutation_free(p);
        gsl_vector_free(x);
        *pX=X;
#endif
}/*}}}*/