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

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

#ifdef _HAVE_ADOLC_
#include "../../shared/Numerics/adolc_edf.h"
#endif

void SolverxSeq(SeqVec<IssmDouble>** puf,SeqMat<IssmDouble>* Kff, SeqVec<IssmDouble>* pf, Parameters* parameters){/*{{{*/

        #ifdef _HAVE_GSL_
        /*Intermediary: */
        int M,N,N2,s;
        SeqVec<IssmDouble> *uf = 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!");
        IssmDouble *x  = xNew<IssmDouble>(N);
#ifdef _HAVE_ADOLC_
        SolverxSeq(x,Kff->matrix,pf->vector,N,parameters);
#else
        SolverxSeq(x,Kff->matrix,pf->vector,N);
#endif
        uf=new SeqVec<IssmDouble>(x,N);
        xDelete(x);

        /*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 *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){//{{{
        // pack inputs to conform to the EDF-prescribed interface
        IssmDouble*  adoubleEDFin=xNew<IssmDouble>(n*(n+1)); // packed inputs, i.e. matrix and right hand side
        for(int i=0; i<n*n;i++)adoubleEDFin[i]    =A[i]; // pack matrix
        for(int i=0; i<n;  i++)adoubleEDFin[i+n*n]=B[i]; // pack the right hand side
        IssmPDouble* pdoubleEDFin=xNew<IssmPDouble>(n*(n+1)); // provide space to transfer inputs during call_ext_fct
        IssmPDouble* pdoubleEDFout=xNew<IssmPDouble>(n);       // provide space to transfer outputs during call_ext_fct
        // call the wrapped solver through the registry entry we retrieve from parameters
        call_ext_fct(dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p,
                     n*(n+1), pdoubleEDFin, adoubleEDFin,
                     n, pdoubleEDFout,X);
        xDelete(adoubleEDFin);
        xDelete(pdoubleEDFin);
        xDelete(pdoubleEDFout);
}
/*}}}*/
#endif
void SolverxSeq(IssmPDouble * X, 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(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*/
        xMemCpy(X,gsl_vector_ptr(x,0),n);

        /*Clean up and assign output pointer*/
        xDelete(Acopy);
        gsl_permutation_free(p);
        gsl_vector_free(x);
        #endif
}
/*}}}*/