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