#include "pargrid.h"

/*+ main - the main routine for this grid program

     Input Parameters:
     argc, argv - the usual
     argv[1] - the number of processors in the x direction
     argv[2] - the number of processors in the y direction
     argv[3] - the number of processors in the z direction
     argv[4] - the number of points in the x direction on each processor
               do not use less than 3 points
     argv[5] - the number of points in the y direction on each processor
               do not use less than 3 points
     argv[6] - the number of points in the z direction on each processor
               do not use less than 3 points
     argv[7] - symmetric matrix = 0, nonsymmetric matrix = 1
     argv[8] - number of components per grid point
     argv[9] - 0 = do not use, 1 = use inode/clique stuff
     argv[10] - 0 = ICC, 1 = ILU 
     argv[11] - 0 => method = CG, 1 => method = GMRES 

     Notes: The grid program solves a linear system associated with
            a 3-D grid distributed across the processors.  The 3-D
            grid is partitioned in all three dimensions among the
            processors.  A 7pt stencil is used.

 +*/

int main(int argc, char **argv)
{
	par_grid grid;
	BSprocinfo *procinfo;
	int	inode, method, scaling, num_rhs;
	int	nonsym_storage, nonsymmetric;

	/* Call BSinit() to initialize BlocklSolve and MPI */
    BSinit(&argc,&argv);

	/* set up the context for BlockSolve */
	procinfo = BScreate_ctx(); CHKERRN(0);
	/* tell it to print out some information on the reordering */
	BSctx_set_pr(procinfo,TRUE); CHKERRN(0);

	if(procinfo->my_id==0) {
		printf("\n");
		printf("************** Blocksolve Example Grid5 *******************\n");
	}

	/* read in grid parameters */
	if (argc < 14) {
		SETERRC(ARG_ERR,"Argument list too small\n");
		if(procinfo->my_id==0) {
		printf("Example: mpirun -np P grid5.ARCH PX PY PZ NX NY NZ SYM NC IN PRE METHOD SCALE NUM_RHS\n");
		printf("ARCH = Machine architecture\n");
		printf("P = number of processors\n");
		printf("PX = the number of processors in the x direction\n");
		printf("PY = the number of processors in the y direction\n");
		printf("PZ = the number of processors in the z direction\n");
		printf("NX = number of points in the x direction on each processor\n");
		printf("NY = number of points in the y direction on each processor\n");
		printf("NZ = number of points in the z direction on each processor\n");
		printf("SYM = 0, generate symmetric matrix, = 1, generate nonsymmetric\n");
		printf("NC = number of components per grid point\n");
		printf("IN = 0, do not use, = 1, use inode/clique stuff\n");
		printf("PRE = 0, compute ICC factor, = 1, compute ILU factor\n");
		printf("METHOD = 0, use PCG, = 1, use GMRES\n");
		printf("SCALE = 0, do not scale matrix, = 1, scale matrix\n");
		printf("NUM_RHS = number of RHS to solve for\n");
		}
		return(-1);
	}
	sscanf(argv[1],"%d",&grid.worker_x);
	sscanf(argv[2],"%d",&grid.worker_y);
	sscanf(argv[3],"%d",&grid.worker_z);
	if (procinfo->my_id == 0) {
		printf("o  ");
		printf("Number of workers (x,y,z): %d %d %d\n",
			grid.worker_x,grid.worker_y,grid.worker_z);
	}
	if (procinfo->nprocs != grid.worker_x*grid.worker_y*grid.worker_z) {
		SETERRC(ARG_ERR,"Number of processors is not correct\n");
		return(-1);
	}
	sscanf(argv[4],"%d",&grid.l_num_x);
	sscanf(argv[5],"%d",&grid.l_num_y);
	sscanf(argv[6],"%d",&grid.l_num_z);
	sscanf(argv[7],"%d",&nonsymmetric);
	sscanf(argv[8],"%d",&grid.ncomp);
	sscanf(argv[9],"%d",&inode);
	sscanf(argv[10],"%d",&nonsym_storage);
	sscanf(argv[11],"%d",&method);
	sscanf(argv[12],"%d",&scaling);
	sscanf(argv[13],"%d",&num_rhs);
	if(grid.ncomp<1) grid.ncomp = 1;
	if(grid.ncomp>20) grid.ncomp = 20;
	grid.symmetric = (!nonsymmetric);
	grid.icc_storage = (!nonsym_storage);
	grid.positive = FALSE;

	/* tell it to use the nonsymmetric matrix data structure */
	if (procinfo->my_id == 0) {
		printf("o  ");
		if(nonsymmetric) 
			printf("Generating nonsymmetric matrix\n");
		else 
			printf("Generating symmetric matrix\n");
	}

	if (procinfo->my_id == 0) {
		printf("o  ");
		printf("Multi component grid, %d components per point\n",grid.ncomp);
	}

	/* decide whether to look for i-nodes or cliques */
	if(inode) {
		BSctx_set_si(procinfo,FALSE); CHKERRN(0);
	} else {
		BSctx_set_si(procinfo,TRUE); CHKERRN(0);
	}

	/* decide on preconditioner */
	if(nonsym_storage) {
		BSctx_set_pre(procinfo,PRE_ILU); CHKERRN(0);
	} else {
		BSctx_set_pre(procinfo,PRE_ICC); CHKERRN(0);
	}

	/* decide on iterative method */
	if(method) {
		BSctx_set_method(procinfo,GMRES); CHKERRN(0);
	} else {
		BSctx_set_method(procinfo,CG); CHKERRN(0);
	}

	BSctx_set_scaling(procinfo,scaling); CHKERRN(0);
	BSctx_set_num_rhs(procinfo,num_rhs); CHKERRN(0);

	/* local grid size and type */
	grid.type = 7;
	if (procinfo->my_id == 0) {
		printf("o  ");
		printf("Local discretizations (x,y,z): %d %d %d\n",
			grid.l_num_x,grid.l_num_y,grid.l_num_z);
	}

	/* call the worker */
	worker(&grid,procinfo); 

	/* print out the log info (if desired) */
	/* MLOG_PRINT(procinfo); */
	
	if(procinfo->my_id==0) {
		printf("************ End Blocksolve Example Grid5 *****************\n");
		printf("\n");
	}

	/* print logging if enabled */
	BSprint_log(procinfo); CHKERRN(0);

	/* free the context */
	BSfree_ctx(procinfo); CHKERRN(0);

	/* finalize BlockSolve and MPI */
	BSfinalize();

	exit(0);
}