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

/*!\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)_error_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !");
        if(M!=N)_error_("Stiffness matrix should be square!");

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

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

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

}/*}}}*/
#ifdef _HAVE_ADOLC_

int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y) {
    if(m*(m+1)!=n)_error_("Stiffness matrix should be square!");
    SolverxSeq(&y,x, x+m*m, m);
    return 0;
}

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;
        IssmPDouble* pdoubleA=NULL;
        IssmPDouble* pdoubleB=NULL;
        IssmPDouble* pdoubleX=NULL;

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

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

        /*Free ressources:*/
        xDelete<IssmPDouble>(pdoubleA);
        xDelete<IssmPDouble>(pdoubleB);

        /*Assign output pointers: */
        *pX=X;
}
/*}}}*/
#endif
void SolverxSeq(IssmPDouble** pX,IssmPDouble* A,IssmPDouble* 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;
        /*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
}
/*}}}*/