/*!\file PlapackInvertMatrix.cpp
 * \brief invert petsc matrix using Plapack
 */

/* petsc: */
#include "../../petsc/petscincludes.h"

/*plapack: */
#include "../plapackincludes.h"

/* Some fortran routines: */
#include "../../scalapack/FortranMapping.h"

void PlapackInvertMatrixLocalCleanup(PLA_Obj* pa,PLA_Template* ptempl,double** parrayA,int** pidxnA,ISSM_MPI_Comm* pcomm_2d);

int PlapackInvertMatrix(Mat* A,Mat* inv_A,int status,int con,COMM comm){ 
	/*inv_A does not yet exist, inv_A was just allocated, that's all*/

	/*Error management*/
	int i;

	/*input*/
	int mA,nA;
	int local_mA,local_nA;
	int lower_row,upper_row;
	MatType type;

	/*Plapack: */
	ISSM_MPI_Datatype   datatype;
	ISSM_MPI_Comm       comm_2d;
	PLA_Obj a=NULL;
	PLA_Template   templ;	
	double one=1.0;		
	int ierror;
	int nb,nb_alg;
	int nprows,npcols;

	/*Petsc to Plapack: */
	double    *arrayA=NULL;
	int* idxnA=NULL;

	/*Verify that A is square*/
	MatGetSize(*A,&mA,&nA);
	MatGetLocalSize(*A,&local_mA,&local_nA);

	/*Some dimensions checks: */
	if (mA!=nA) _error_("trying to take the invert of a non-square matrix!");

	/* Set default Plapack parameters */
	//First find nprows*npcols=num_procs;
	CyclicalFactorization(&nprows,&npcols,num_procs); 
	//nprows=num_procs;
	//npcols=1;
	ierror = 0;
	nb     = nA/num_procs;
	if(nA - nb*num_procs) nb++; /* without cyclic distribution */

	if (ierror){
		PLA_Set_error_checking(ierror,PETSC_TRUE,PETSC_TRUE,PETSC_FALSE );
	}
	else {
		PLA_Set_error_checking(ierror,PETSC_FALSE,PETSC_FALSE,PETSC_FALSE );
	}
	nb_alg = 0;
	if (nb_alg){
		pla_Environ_set_nb_alg (PLA_OP_ALL_ALG,nb_alg);
	}

	/*Verify that plapack is not already initialized: */
	if(PLA_Initialized(NULL)==TRUE)PLA_Finalize();
	/* Create a 2D communicator */
	PLA_Comm_1D_to_2D(comm,nprows,npcols,&comm_2d); 

	/*Initlialize plapack: */
	PLA_Init(comm_2d);

	templ = NULL;
	PLA_Temp_create(nb, 0, &templ);

	/* Use suggested nb_alg if it is not provided by user */
	if (nb_alg == 0){
		PLA_Environ_nb_alg(PLA_OP_PAN_PAN,templ,&nb_alg);
		pla_Environ_set_nb_alg(PLA_OP_ALL_ALG,nb_alg);
	}

	/* Set the datatype */
	datatype = ISSM_MPI_DOUBLE;

	/* Copy A into a*/
	PLA_Matrix_create(datatype,mA,nA,templ,PLA_ALIGN_FIRST,PLA_ALIGN_FIRST,&a);  
	PLA_Obj_set_to_zero(a);
	/*Take array from A: use MatGetValues, because we are sure this routine works with 
	any matrix type.*/
	MatGetOwnershipRange(*A,&lower_row,&upper_row);
	upper_row--; 
	arrayA = xNew<double>(nA);
	idxnA  = xNew<int>(nA);
	for (i=0;i<nA;i++){
		*(idxnA+i)=i;
	}
	PLA_API_begin();
	PLA_Obj_API_open(a);  	
	for (i=lower_row;i<=upper_row;i++){
		MatGetValues(*A,1,&i,nA,idxnA,arrayA); 
		PLA_API_axpy_matrix_to_global(1,nA, &one,(void *)arrayA,1,a,i,0); 
	}
	PLA_Obj_API_close(a); 
	PLA_API_end(); 

	/*Call the plapack invert routine*/
	PLA_General_invert(PLA_METHOD_INV,a);

	/*Translate Plapack a into Petsc invA*/
	MatGetType(*A,&type);
	PlapackToPetsc(inv_A,local_mA,local_nA,mA,nA,type,a,templ,nprows,npcols,nb);

	/*Free ressources:*/
	PLA_Obj_free(&a);
	PLA_Temp_free(&templ);
	xDelete<double>(arrayA);
	xDelete<int>(idxnA);

	/*Finalize PLAPACK*/
	PLA_Finalize();
	ISSM_MPI_Comm_free(&comm_2d);
}
