/* SMP:
export MPI=...
export fig=...
$MPIEXEC -np 2 ./solve.bin -f0 $fig/opciones.ejecu -f1 $fig/temporal/femsii_a -f2 $fig/temporal/femsii_b -f3 $fig/temporal/femsii_x
*/

/* 
  Include "petscksp.h" so that we can use KSP solvers.  Note that this file
  automatically includes:
     petsc.h       - base PETSc routines   petscvec.h - vectors
     petscsys.h    - system routines       petscmat.h - matrices
     petscis.h     - index sets            petscksp.h - Krylov subspace methods
     petscviewer.h - viewers               petscpc.h  - preconditioners


*/

/*!\file femsii_solve.cpp
\brief
 * Este programa acepta un sistema de ecuaciones, lo resuelve segun las
 * indicaciones que se le dan, y escribe un ventor de resultados.
 *
 * Los archivos de entrada son 4:
 * -f0 -> opciones de ejecucion, aplicables a la resolucion del sistema.
 * -f1 -> matriz de coheficientes.
 * -f2 -> vector de terminos independientes.
 * -f3 -> archivo de salida para femsii_x con la resolucion del sistema.
 *
 * DESCRIPCION DE main(int argc, char **args): Se encarga de entender los parametros de entrada, rescatar la matriz de rigidez y el vector de terminos independientes, y resolver el sistema de ecuaciones. solve() lee el archivo opciones.ejecu para cargar las opcions con las que se lanzara el resolutor de subespacios de Krilov. Tras crear un vector de resultados lo escribe en la carpeta temporal del proyecto. La estructura de de la funcion y de el codigo de este archivo en general ha sido creado pensado en crear un codigo siguiendo las directrices generales de la documentacion de PETSc para el caso de sistemas paralelizables. De esta manera solbe.bin puede ser ejecutado en varios procesos concurrentes.
 */
static char help[] = "Solves [A][x]=[b] system.\n\n";

#include "petscksp.h"

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
	KSP ksp;
	PetscErrorCode ierr;
	Mat A;
	Vec x,b;
	PetscViewer fd;
	char file[4][PETSC_MAX_PATH_LEN];
	PetscInt iteraciones, m, n;
	PetscMPIInt rank, size;
	PetscBool flg;
//	PetscLogDouble tsolve,tsolve1,tsolve2;

	PetscInitialize(&argc,&args,(char *)0,help);
	MPI_Comm_rank(PETSC_COMM_WORLD,&rank); // cada procesador recibira un valor distinto
	MPI_Comm_size(PETSC_COMM_WORLD,&size);
	PetscPrintf(PETSC_COMM_WORLD,"Number of processors = %d, rank = %d\n",size,rank);

	PetscOptionsGetString(PETSC_NULL,"-f0",file[0],PETSC_MAX_PATH_LEN-1,&flg); //options

	PetscOptionsGetString(PETSC_NULL,"-f1",file[1],PETSC_MAX_PATH_LEN-1,&flg); //Matriz
	PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&fd);
	MatCreate(PETSC_COMM_WORLD,&A);
	MatSetType(A,MATSBAIJ);
	MatLoad(A,fd); // MatLoad(fd,MATSBAIJ,&A);
	PetscViewerDestroy(&fd);
	MatSetFromOptions(A);

	PetscOptionsGetString(PETSC_NULL,"-f2",file[2],PETSC_MAX_PATH_LEN-1,&flg); //Vector
	PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&fd);
	VecCreate(PETSC_COMM_WORLD,&b);
	VecLoad(b,fd);
	PetscViewerDestroy(&fd);

	PetscOptionsGetString(PETSC_NULL,"-f3",file[3],PETSC_MAX_PATH_LEN-1,&flg); //destiny

	MatGetSize(A,&m,&n);
//	if (m != n) {SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_SIZ, "Not for rectangular matrices (%d, %d)\n\n", m, n);}
	VecGetSize(b,&m);
//	if (m != n) {SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_ARG_SIZ, "Vector b incoherente %d, %d\n\n", m, n);}
	PetscPrintf(PETSC_COMM_WORLD,"Tamanho de x: %d\n", n);

	VecCreate(PETSC_COMM_WORLD,&x);
	VecSetSizes(x,PETSC_DECIDE,m);
	VecSetFromOptions(x);
	KSPCreate(PETSC_COMM_WORLD,&ksp);
	KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
	PetscOptionsInsertFile(PETSC_COMM_WORLD,file[0],PETSC_TRUE);
	KSPSetFromOptions(ksp);
	ierr = KSPSolve(ksp,b,x); CHKERRQ(ierr);
	VecAssemblyBegin(x);
	VecAssemblyEnd(x);

	KSPGetIterationNumber(ksp,&iteraciones);
	PetscPrintf(PETSC_COMM_WORLD,"  Numero de iteraciones = %3D\n",iteraciones);
//	PetscGetTime(&tsolve2);
//	tsolve = tsolve2 - tsolve1;
//	PetscPrintf(PETSC_COMM_WORLD,"Tiempo KSP: %f\n", tsolve2-tsolve1);

	PetscPrintf(PETSC_COMM_WORLD,"Escribiendo binario del vector de resultados.\n");
	PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[3],FILE_MODE_WRITE,&fd);
	VecView(x,fd);
	PetscViewerDestroy(&fd);

	MatDestroy(&A);
	VecDestroy(&b);
	VecDestroy(&x);
	KSPDestroy(&ksp);
	PetscFinalize();
	return 0;
}
