Need Help : FGMRES ILU0, Diagonal and RCI iterative solvers

[Naoufel_E_]

I am trying to put together an FGMRES iterative sparse solver routine, that uses one of the examples given by INTEL, (I have tried many technics; the ILU, Diagonal preconditioner and the RCI; Reverse Communication Interface - the code is being written in c++).

I am having some trouble with the results of all the solver types above, the given solution seems does contain the correct results, and want to rule out a problem in how I am applying the preconditioner or the RCI.

My sparse matrix is in a compressed row format with a column index array (ja), row pointer array (ia), and values (ar) and the right-hand side array (r_rhs) as shown in the MKL documentation. The indices in ja and ia are in the Fortran format, so they start with 1. m_dResidualValue is the tolerance (1.0e-7) in the default case, n is the number of rows, nz is the number of the non-zeros.

Please find the code used in the ILU preconditioner, most of the code is from the FGMRES with ILU preconditioner examples, with a few small modifications.  Does this code look like the correct way to apply a diagonal preconditioner during the FGMRES solve with a preconditioner? I have included the entire iterative solve method I'm attempting for reference (again, most of the code here is directly from the FGMRES solver examples, so much of it should be correct). Obviously, portions of this code need to be cleaned up, but please bear with me as this is a first attempt at getting the solver to work.

int FGMRES_Iterative_Solver_ILU_Precond(CString path)
{

	int bSucceeded = 1;


	/*---------------------------------------------------------------------------
	/* Allocate storage for the ?par parameters and the solution/r_rhs/residual vectors
	/*---------------------------------------------------------------------------*/
	double *b, *computed_solution, *residual;
	double *trvec, *bilu0, *dpar, *tmp;; //trvec[*n],bilu0[*nz];
	MKL_INT *ipar, i;
	// allcoate memory 

	b = (double *)calloc(n, sizeof(double));
	computed_solution = (double *)calloc(n, sizeof(double));
	residual = (double *)calloc(n, sizeof(double));

	trvec = (double *)calloc(n, sizeof(double));
	bilu0 = (double *)calloc(nz, sizeof(double));
	ipar = (MKL_INT *)calloc(128, sizeof(MKL_INT));

	unsigned int iSize = n;
	MKL_INT sizetmp = (n)*(2 * (n)+1) + ((n)*((n)+9)) / 2 + 1;
	printf("sizetmpi= %d, sizeof(ia)= %d\n", sizetmp, sizeof(ia));

	//double dpar[size];
	// double tmp[sizetmp];
	dpar = (double *)calloc(128, sizeof(double));
	tmp = (double *)calloc(sizetmp, sizeof(double));
	printf("n=%d, nz= %d\n", n, nz);
	double tol = 1.0e-6;
	MKL_INT matsize = nz, incx = 1;
	double nrm2;

	//for(i=0; i<=*n; i++) --ia;
	//for(i=0; i<*nz; i++) --ja;

	/*---------------------------------------------------------------------------
	/* Some additional variables to use with the (P)FGMRES solver
	/*---------------------------------------------------------------------------*/
	MKL_INT itercount, ierr = 0;
	MKL_INT RCI_request, ivar;
	printf("LinearSolver: Line 48\n");
	double dvar;
	char cvar, cvar1, cvar2;


	/*---------------------------------------------------------------------------
	/* Initialize variables and the right hand side through matrix-vector product
	/*---------------------------------------------------------------------------*/
	ivar = n;
	cvar = 'n';
	/*---------------------------------------------------------------------------
	/* Save the right-hand side in vector b for future use
	/*---------------------------------------------------------------------------*/
	i = 1;
	dcopy(&n, r_rhs, &i, b, &i);

	//clbas_dcopy(ivar, r_rhs, i, b, i);
	/*---------------------------------------------------------------------------
	/* Initialize the initial guess
	/*---------------------------------------------------------------------------*/
	for (i = 0; i < n; i++)
		computed_solution = 0.0;
	//computed_solution[0] = 1.0;


	/*---------------------------------------------------------------------------
	/* Initialize the solver
	/*---------------------------------------------------------------------------*/
	printf("dfgmres_init\n");
	dfgmres_init(&ivar, computed_solution, r_rhs, &RCI_request, ipar, dpar, tmp);
	printf("dfgmres_end\n");
	if (RCI_request != 0)
		goto FAILED;

	/*---------------------------------------------------------------------------
	/* Calculate ILU0 preconditioner.
	/*                      !ATTENTION!
	/* DCSRILU0 routine uses some IPAR, DPAR set by DFGMRES_INIT routine.
	/* Important for DCSRILU0 default entries set by DFGMRES_INIT are
	/* ipar[1] = 6 - output of error messages to the screen,
	/* ipar[5] = 1 - allow output of errors,
	/* ipar[30]= 0 - abort DCSRILU0 calculations if routine meets zero diagonal element.
	/*
	/* If ILU0 is going to be used out of MKL FGMRES context, than the values
	/* of ipar[1], ipar[5], ipar[30], dpar[30], and dpar[31] should be user
	/* provided before the DCSRILU0 routine call.
	/*
	/* In this example, specific for DCSRILU0 entries are set in turn:
	/* ipar[30]= 1 - change small diagonal value to that given by dpar[31],
	/* dpar[30]= 1.E-20 instead of the default value set by DFGMRES_INIT.
	/*                  It is a small value to compare a diagonal entry with it.
	/* dpar[31]= 1.E-16 instead of the default value set by DFGMRES_INIT.
	/*                  It is the target value of the diagonal value if it is
	/*                  small as compared to dpar[30] and the routine should change
	/*                  it rather than abort DCSRILU0 calculations.
	/*---------------------------------------------------------------------------*/

	//ipar[30] = 1;
	//dpar[30] = 1.E-50;
	//dpar[31] = 1.E-50;

	printf("ilu begin\n");
	dcsrilu0(&n, &ar[0], &ia[0], &ja[0], &bilu0[0], &ipar[0], &dpar[0], &ierr);
	printf("ilu end\n");

	nrm2 = dnrm2(&matsize, bilu0, &incx);

	if (ierr != 0)
	{
		printf("Preconditioner dcsrilu0 has returned the ERROR code %d", ierr);
		goto FAILED1;
	}

	/*---------------------------------------------------------------------------
	/* Set the desired parameters:
	/* do the restart after 2 iterations
	/* LOGICAL parameters:
	/* do not do the stopping test for the maximal number of iterations
	/* do the Preconditioned iterations of FGMRES method
	/* Set parameter ipar[10] for preconditioner call. For this example,
	/* it reduces the number of iterations.
	/* DOUBLE PRECISION parameters
	/* set the relative tolerance to 1.0D-3 instead of default value 1.0D-6
	/* NOTE. Preconditioner may increase the number of iterations for an
	/* arbitrary case of the system and initial guess and even ruin the
	/* convergence. It is user's responsibility to use a suitable preconditioner
	/* and to apply it skillfully.
	/*---------------------------------------------------------------------------*/
	ipar[14] = 2;
	ipar[7] = 0;
	ipar[10] = 1;
	dpar[0] = m_dResidualValue;

	/*---------------------------------------------------------------------------
	/* Check the correctness and consistency of the newly set parameters
	/*---------------------------------------------------------------------------*/
	printf("dfgmres check begin\n");
	dfgmres_check(&ivar, computed_solution, r_rhs, &RCI_request, ipar, dpar, tmp);
	if (RCI_request != 0) goto FAILED;
	printf("dfgmres check end\n");

	/*---------------------------------------------------------------------------
	/* Compute the solution by RCI (P)FGMRES solver with preconditioning
	/* Reverse Communication starts here
	/*---------------------------------------------------------------------------*/
	printf("dfgmres begin to solve\n");
ONE:  dfgmres(&ivar, computed_solution, r_rhs, &RCI_request, ipar, dpar, tmp);

	printf("dfgmres end of solve\n");

	/*---------------------------------------------------------------------------
	/* If RCI_request=0, then the solution was found with the required precision
	/*---------------------------------------------------------------------------*/
	if (RCI_request == 0) goto COMPLETE;
	/*---------------------------------------------------------------------------
	/* If RCI_request=1, then compute the vector A*tmp[ipar[21]-1]
	/* and put the result in vector tmp[ipar[22]-1]
	/*---------------------------------------------------------------------------
	/* NOTE that ipar[21] and ipar[22] contain FORTRAN style addresses,
	/* therefore, in C code it is required to subtract 1 from them to get C style
	/* addresses
	/*---------------------------------------------------------------------------*/
	if (RCI_request == 1)
	{
		mkl_dcsrgemv(&cvar, &ivar, ar, ia, ja, &tmp[ipar[21] - 1], &tmp[ipar[22] - 1]);
		goto ONE;
	}
	/*---------------------------------------------------------------------------
	/* If RCI_request=2, then do the user-defined stopping test
	/* The residual stopping test for the computed solution is performed here
	/*---------------------------------------------------------------------------
	/* NOTE: from this point vector b is no longer containing the right-hand
	/* side of the problem! It contains the current FGMRES approximation to the
	/* solution. If you need to keep the right-hand side, save it in some other
	/* vector before the call to dfgmres routine. Here we saved it in vector
	/* r_rhs. The vector b is used instead of r_rhs to preserve the
	/* original right-hand side of the problem and guarantee the proper
	/* restart of FGMRES method. Vector b will be altered when computing the
	/* residual stopping criterion!
	/*---------------------------------------------------------------------------*/
	if (RCI_request == 2)
	{
		/* Request to the dfgmres_get routine to put the solution into b via ipar[12]
		/*---------------------------------------------------------------------------
		/* WARNING: beware that the call to dfgmres_get routine with ipar[12]=0 at this stage may
		/* destroy the convergence of the FGMRES method, therefore, only advanced users should
		/* exploit this option with care */
		ipar[12] = 1;
		/* Get the current FGMRES solution in the vector b */
		dfgmres_get(&ivar, computed_solution, b, &RCI_request, ipar, dpar, tmp, &itercount);
		/* Compute the current true residual via MKL (Sparse) BLAS routines */
		mkl_dcsrgemv(&cvar, &ivar, ar, ia, ja, b, residual);
		dvar = -1.0E0;
		i = 1;
		daxpy(&ivar, &dvar, r_rhs, &i, residual, &i);
		dvar = dnrm2(&ivar, residual, &i);
		if (dvar < m_dResidualValue)
			goto COMPLETE;

		else goto ONE;
	}
	/*---------------------------------------------------------------------------
	/* If RCI_request=3, then apply the preconditioner on the vector
	/* tmp[ipar[21]-1] and put the result in vector tmp[ipar[22]-1]
	/*---------------------------------------------------------------------------
	/* NOTE that ipar[21] and ipar[22] contain FORTRAN style addresses,
	/* therefore, in C code it is required to subtract 1 from them to get C style
	/* addresses
	/* Here is the recommended usage of the result produced by ILU0 routine
	/* via standard MKL Sparse Blas solver routine mkl_dcsrtrsv.
	/*---------------------------------------------------------------------------*/
	if (RCI_request == 3)
	{
		cvar1 = 'L';
		cvar = 'n';
		cvar2 = 'U';
		mkl_dcsrtrsv(&cvar1, &cvar, &cvar2, &ivar, bilu0, ia, ja, &tmp[ipar[21] - 1], trvec);
		cvar1 = 'U';
		cvar = 'n';
		cvar2 = 'n';
		mkl_dcsrtrsv(&cvar1, &cvar, &cvar2, &ivar, bilu0, ia, ja, trvec, &tmp[ipar[22] - 1]);
		goto ONE;
	}

	/*---------------------------------------------------------------------------
	/* If RCI_request=4, then check if the norm of the next generated vector is
	/* not zero up to rounding and computational errors. The norm is contained
	/* in dpar[6] parameter
	/*---------------------------------------------------------------------------*/
	if (RCI_request == 4)
	{
		if (dpar[6] < 1.0E-50)
			goto COMPLETE;
		else goto ONE;
	}
	/*---------------------------------------------------------------------------
	/* If RCI_request=anything else, then dfgmres subroutine failed
	/* to compute the solution vector: computed_solution
	/*---------------------------------------------------------------------------*/
	else
	{
		goto FAILED;
	}
	/*---------------------------------------------------------------------------
	Reverse Communication ends here
	Get the current iteration number and the FGMRES solution (DO NOT FORGET to
	call dfgmres_get routine as computed_solution is still containing
	the initial guess!). Request to dfgmres_get to put the solution
	into vector computed_solution via ipar[12]
	---------------------------------------------------------------------------*/
COMPLETE:   ipar[12] = 0;
	dfgmres_get(&ivar, computed_solution, r_rhs, &RCI_request, ipar, dpar, tmp, &itercount);
	/*---------------------------------------------------------------------------
	Print solution vector: computed_solution and the number of iterations: itercount
	---------------------------------------------------------------------------*/
	printf("The system has been solved \n");
	/*for (i = 0; i < n; i++)
	r_rhs[i + 1] = computed_solution;*/
	writeSolution(path, computed_solution, 1, n);
	printf("\nNumber of iterations: %d\n", itercount);
	printf("\n");
	goto CLEANMEMORY;

FAILED:
	bSucceeded = 0;
	printf("The solver has returned the ERROR code %d \n", RCI_request);
FAILED1:
	bSucceeded = 0;
	printf("-------------------------------------------------------------------\n");
	printf("Unfortunately, FGMRES+ILU0 C example has FAILED\n");
	printf("-------------------------------------------------------------------\n");

CLEANMEMORY:
	free(b);
	free(computed_solution);
	free(residual);
	free(trvec);
	free(bilu0);
	free(tmp);
	free(dpar);
	free(ipar);

	return bSucceeded;
}