/*!\file SolverxGsl * \brief Gsl implementation of solver */ #ifdef HAVE_CONFIG_H #include #else #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" #endif #include #include #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(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*/ IssmDouble* X = xNew(n); memcpy(X,gsl_vector_ptr(x,0),n*sizeof(IssmDouble)); /*Clean up and assign output pointer*/ xDelete(Acopy); gsl_permutation_free(p); gsl_vector_free(x); *pX=X; #endif }/*}}}*/