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

}/*}}}*/
void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/

	/*Allocate output*/
	double* X = xNew<double>(n); 

	/*Solve*/
	SolverxSeq(X,A,B,n);

	/*Assign output pointer*/
	*pX=X;
}
/*}}}*/

#ifdef _HAVE_ADOLC_
int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y){ /*{{{*/
    SolverxSeq(y,x, x+m*m, m); // x is where the matrix starts, x+m*m is where the right-hand side starts
    return 0;
} /*}}}*/
int EDF_fos_forward_for_solverx(int n, IssmPDouble *inVal, IssmPDouble *inDeriv, int m, IssmPDouble *outVal, IssmPDouble *outDeriv) { /*{{{*/
#ifdef _HAVE_GSL_
  // the matrix will be modified by LU decomposition. Use gsl_A copy
  double* Acopy = xNew<double>(m*m);
  xMemCpy(Acopy,inVal,m*m);
  /*Initialize gsl matrices and vectors: */
  gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
  gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
  gsl_permutation *perm_p = gsl_permutation_alloc (m);
  int  signPerm;
  // factorize
  gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
  gsl_vector *gsl_x_p = gsl_vector_alloc (m);
  // solve for the value
  gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
  /*Copy result*/
  xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
  gsl_vector_free(gsl_x_p);
  // solve for the derivatives acc. to A * dx = r  with r=db - dA * x
  // compute the RHS
  double* r=xNew<double>(m);
  for (int i=0; i<m; i++) {
    r[i]=inDeriv[m*m+i]; // this is db[i]
    for (int j=0;j<m; j++) {
      r[i]-=inDeriv[i*m+j]*outVal[j]; // this is dA[i][j]*x[j]
    }
  }
  gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
  gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
  gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
  xMemCpy(outDeriv,gsl_vector_ptr(gsl_dx_p,0),m);
  gsl_vector_free(gsl_dx_p);
  xDelete(r);
  gsl_permutation_free(perm_p);
  xDelete(Acopy);
  #endif
  return 0;
} /*}}}*/
int EDF_fov_forward_for_solverx(int n, IssmPDouble *inVal, int directionCount, IssmPDouble **inDeriv, int m, IssmPDouble *outVal, IssmPDouble **outDeriv) { /*{{{*/
#ifdef _HAVE_GSL_
  // the matrix will be modified by LU decomposition. Use gsl_A copy
  double* Acopy = xNew<double>(m*m);
  xMemCpy(Acopy,inVal,m*m);
  /*Initialize gsl matrices and vectors: */
  gsl_matrix_view gsl_A = gsl_matrix_view_array (Acopy,m,m);
  gsl_vector_view gsl_b = gsl_vector_view_array (inVal+m*m,m); // the right hand side starts at address inVal+m*m
  gsl_permutation *perm_p = gsl_permutation_alloc (m);
  int  signPerm;
  // factorize
  gsl_linalg_LU_decomp (&gsl_A.matrix, perm_p, &signPerm);
  gsl_vector *gsl_x_p = gsl_vector_alloc (m);
  // solve for the value
  gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_b.vector, gsl_x_p);
  /*Copy result*/
  xMemCpy(outVal,gsl_vector_ptr(gsl_x_p,0),m);
  gsl_vector_free(gsl_x_p);
  // solve for the derivatives acc. to A * dx = r  with r=db - dA * x
  double* r=xNew<double>(m);
  gsl_vector *gsl_dx_p = gsl_vector_alloc(m);
  for (int dir=0;dir<directionCount;++dir) {
    // compute the RHS
    for (int i=0; i<m; i++) {
      r[i]=inDeriv[m*m+i][dir]; // this is db[i]
      for (int j=0;j<m; j++) {
        r[i]-=inDeriv[i*m+j][dir]*outVal[j]; // this is dA[i][j]*x[j]
      }
    }
    gsl_vector_view gsl_r=gsl_vector_view_array(r,m);
    gsl_linalg_LU_solve (&gsl_A.matrix, perm_p, &gsl_r.vector, gsl_dx_p);
    // reuse r
    xMemCpy(r,gsl_vector_ptr(gsl_dx_p,0),m);
    for (int i=0; i<m; i++) {
      outDeriv[i][dir]=r[i];
    }
  }
  gsl_vector_free(gsl_dx_p);
  xDelete(r);
  gsl_permutation_free(perm_p);
  xDelete(Acopy);
  #endif
  return 0;
}
/*}}}*/
int EDF_fos_reverse_for_solverx(int m, double *dp_U, int n, double *dp_Z) { /*{{{*/
  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
}
/*}}}*/