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

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

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

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

}/*}}}*/
void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* 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;
		#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_("SolverxSeq: 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
	#endif
}/*}}}*/
