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

/*!\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

}/*}}}*/
#ifdef _HAVE_ADOLC_
void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n){//{{{
        /* if we use Adol-C then the IssmDouble will be an adouble
           and the calls to gsl_... will not work. 
           We therefore call a wrapped solver instead. 
        */

        /*Output: */
        IssmDouble* X=NULL;

        /*Intermediary: */
        int     i;
        double* doubleA=NULL;
        double* doubleB=NULL;
        double* doubleX=NULL;

        /*First, transfer from IssmDouble to double all our matrices and vectors: */
        doubleA=xNew<double>(n*n);
        doubleB=xNew<double>(n);
        for(i=0;i<n*n;i++)A[i]>>=doubleA[i];
        for(i=0;i<n;i++)B[i]>>=doubleB[i];
        
        /*Call wrapped solver: */
        SolverxSeq(&doubleX,doubleA, doubleB, n);

        /*Transfer solution vector from double to IssmDouble: */
        X = xNew<IssmDouble>(n);
        for(i=0;i<n;i++)X[i]<<=doubleX[i];

        /*Free ressources:*/
        xDelete<double>(doubleA);
        xDelete<double>(doubleB);

        /*Assign output pointers: */
        *pX=X;
}
/*}}}*/
#endif
//void SolverxSeq(double** pX,double* A,double* B,int n){ //{{{
#ifdef _HAVE_ADOLC_
void SolverxSeq(double** pX,double* A,double* B,int n){ 
#else
void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n){
#endif
        #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;
        /*A will be modified by LU decomposition. Use copy*/
        double* Acopy = xNew<double>(n*n);
        xMemCpy<double>(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*/
        double* X = xNew<double>(n);
        memcpy(X,gsl_vector_ptr(x,0),n*sizeof(double));

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