/*!\file SolverxGsl
 * \brief Gsl implementation of solver
 */

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

#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<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
}/*}}}*/