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

/*!\file SolverxSeq
 * \brief implementation of sequential solver using the GSL librarie
 */

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

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

#ifdef _HAVE_GSL_
#include <gsl/gsl_linalg.h>
#endif

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

        #ifdef _HAVE_GSL_
        /*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!");

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

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

        #else
                _error2_("GSL support not compiled in!");
        #endif

}/*}}}*/
void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n){/*{{{*/

        #ifdef _HAVE_GSL_
                /*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_("SolverxSeq: 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
        #endif
}/*}}}*/