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