#include <stdlib.h>
#include "stack-c.h"
#include "api_scilab.h"
#include "Scierror.h"
#include "sciprint.h"

#include "Kriging.h"
#include "Erreur.h"

/**********************************
 * SCILAB function : GaussHermite *
 **********************************/

int gtw_GaussHermite ( char * fname )
{
	int m, n, l;
    	int N;
	int un = 1;
	double * X , * W ;
	
	CheckRhs(1,1);
	CheckLhs(2,2);

	/* Variable N */

	GetRhsVar(1,"i",&m,&n,&l);
	CheckScalar(1,m,n);
    	N = *istk(l);
	
    	/* Appel a l'API */

    	if ( GaussHermite ( & X , & W , N ) )
        	Scierror(999,GetErreur());

	/* Creation des vecteurs (lignes) d'abscisses et de poids sur le stack */	
		
    	CreateVarFromPtr ( Rhs + 1 , "d", & un , & N , & X );
	LhsVar(1) = Rhs + 1;

    	CreateVarFromPtr ( Rhs + 2 , "d", & un , & N , & W );
	LhsVar(2) = Rhs + 2;
	
	free(X);
	free(W);
	
    	return(0);
}

/********************************************
 * SCILAB function : KrigCholesky           *
 ********************************************/

int gtw_KrigCholesky ( char * fname )
{
	int m, n, l;
	int lda,iuplo,info ;
	double * A;

	CheckRhs(1,1);
	CheckLhs(1,1);

	/* Symmetric definite positive matrix */	
	
	GetRhsVar(1,"d",&m,&n,&l);
	CheckDimProp(1, 1, m != n );
	A = stk(l);
	lda = m;

	/* Compute Cholesky factorization */

	iuplo = 0;

	(void) mydpotrf_ ( & iuplo , & lda , A , & lda , & info ) ;
	
	/* Return factorized matrix */

	CreateVarFromPtr ( Rhs + 1 , "d", & lda , & lda , & A );
	LhsVar(1) = Rhs + 1;
	
	return(0);	
}

/********************************************
 * SCILAB function : KrigSolve              *
 ********************************************/

int gtw_KrigSolve ( char * fname )
{
	int m, n, l;
	int lda,nrhs,iuplo,info ;
	double * A;
	double * B;

	CheckRhs(2,2);
	CheckLhs(1,1);

	/* Cholesky-factorized covariance matrix */	
	
	GetRhsVar(1,"d",&m,&n,&l);
	CheckDimProp(1, 1, m != n );
	A = stk(l);
	lda = m;

	/* Right-Hand Side */

	GetRhsVar(2,"d",&m,&n,&l);
	CheckDimProp(1, 2, m != lda );
	B = stk(l);
	nrhs = n;

	/* Solve sparse system */

	iuplo = 0;

	(void) mydpotrs_ ( & iuplo , & lda , & nrhs , A , & lda , B , & lda , & info ) ;
	
	/* Return solution on the stack */

	CreateVarFromPtr ( Rhs + 1 , "d", & lda , & nrhs , & B );
	LhsVar(1) = Rhs + 1;
	
	return(0);	
}

/********************************************
 * SCILAB function : KrigBunchKaufman       *
 ********************************************/

int gtw_KrigBunchKaufman ( char * fname )
{
	int m, n, l;
	int un = 1;
	int lda,iuplo,info,lwork;
	int * IPIV ;
	double * WORK ;
	double * A;

	CheckRhs(1,1);
	CheckLhs(2,2);

	/* Symmetric matrix */	
	
	GetRhsVar(1,"d",&m,&n,&l);
	CheckDimProp(1, 1, m != n );
	A = stk(l);
	lda = m;

    /* Vector IPIV */
	
	IPIV = (int *) malloc ( lda * sizeof(int) );
	
	lwork = lda * lda ;
	WORK = (double *) malloc ( lda * lda * sizeof(double) );
			
	/* Compute Bunch-Kaufman factorization */

	iuplo = 0;

	(void) mydsytrf_ ( & iuplo , & lda , A , & lda , IPIV , WORK , & lwork , & info ) ;
	
	/* Return factorized matrix and permutation vector */

	CreateVarFromPtr ( Rhs + 1 , "d", & lda , & lda , & A );
	LhsVar(1) = Rhs + 1;
	
	CreateVarFromPtr ( Rhs + 2 , "i", & lda , & un , & IPIV );
	LhsVar(2) = Rhs + 2;
	
	free(IPIV);
	free(WORK);
	
	return(0);	
}

/********************************************
 * SCILAB function : KrigSolve              *
 ********************************************/

int gtw_KrigSolveOK ( char * fname )
{
	int m, n, l;
	int lda,nrhs,iuplo,info ;
	int * IPIV ;
	double * A;
	double * B;

	CheckRhs(3,3);
	CheckLhs(1,1);

	/* LU factorized kriging matrix */	
	
	GetRhsVar(1,"d",&m,&n,&l);
	CheckDimProp(1, 1, m != n );
	A = stk(l);
	lda = m;

	/* Permutation vector */
	
	GetRhsVar(2,"i",&m,&n,&l);
	CheckVector(2,m,n);
	CheckDimProp(1, 2, ( m * n )!= lda );
	IPIV = istk(l);
		
	/* Right-Hand Side */

	GetRhsVar(3,"d",&m,&n,&l);
	CheckDimProp(1, 3, m != lda );
	B = stk(l);
	nrhs = n;

	/* Solve system */

	iuplo = 0;

	(void) mydsytrs_ ( & iuplo , & lda , & nrhs , A , & lda , IPIV , B , & lda , & info ) ;
	
	/* Return solution on the stack */

	CreateVarFromPtr ( Rhs + 1 , "d", & lda , & nrhs , & B );
	LhsVar(1) = Rhs + 1;
	
	return(0);	
}


