FGMRES + Diagonal preconditioner. Is this correct?

[mdobossy]

I am trying to put together a FGMRES iterative sparse solver routine, that uses a diagonal preconditioner (the code is being written in c). I am having some trouble with the convergence of the solution, and want to rule out a problem in how I am applying the preconditioner. Most of the code is from the FGMRES with ILU preconditioner examples, with a few small modifications.

My sparse matrix is in a compressed row format with a column index array (_col_ind), row pointer array (_row_ptr), and values (_data) as shown in the MKL documentation. The indices in _col_ind and _row_ptr are in fortran format, so they start with 1. These arrays are STL::vector containers (thus my referencing the memory address of the first element of interest). ivar is the number of rows.

Does this code look like the correct way to apply a diaganal preconditioner during the FGMRES solve with a diagonal preconditioner using the mkl_dcsrsv function? I have bolded the application of the preconditioner and initial guess, but 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.

vector SPMat::iterativeSolve(vector &RHS) {
if(_format!=CRS) //Checks to see if the SPMat has been formatted as a CRS yet
toCRS();
solveError=false;
status statusCode=CONTINUE;
vector solution(_rows);
MKL_INT ipar[128];
double dpar[128], *tmp=new double [_rows*(2*_rows+1)+(_rows*(_rows+9))/2+1];
double *trvec = new double[_rows];
double *b = new double[_rows];
double *residual = new double[_rows];

MKL_INT itercount,ierr=0;
MKL_INT RCI_request, i, ivar;
double dvar;
ivar=_rows;
double alpha=1.0;
i=1;
dcopy(&ivar, &RHS[0], &i, b, &i);

/*---------------------------------------------------------------------------
/* Make an initial guess- inverse of the diagonal * RHS.
/*---------------------------------------------------------------------------*/
mkl_dcsrsv("N", &ivar, α, "DLNF ", &_data[0], &_col_ind[0], &_row_ptr[0], &_row_ptr[1], b, &solution.front());

/*---------------------------------------------------------------------------
/* Initi alize the solver
/*---------------------------------------------------------------------------*/
dfgmres_init(&ivar, &solution.front(), b, &RCI_request, ipar, dpar, tmp);
if (RCI_request!=0) { cout << "dfgmres_init failed." << endl; statusCode=FAILED; }

ipar[30]=1;
dpar[30]=1.E-5;
ipar[14]=2;
ipar[7]=0;
ipar[8]=1;
ipar[9]=0;
ipar[10]=1;
ipar[11]=1;
dpar[0]=1.0E-6;

/*---------------------------------------------------------------------------
/* Check the correctness and consistency of the newly set parameters
/*---------------------------------------------------------------------------*/
dfgmres_check(&ivar, &solution.front(), b, &RCI_request, ipar, dpar, tmp);
if (RCI_request!=0) { cout << "dfgmres_check failed." << endl; statusCode=FAILED; }

while(statusCode==CONTINUE) {
//Compute the solution by RCI (P)FGMRES solver with preconditioning
dfgmres(&ivar, &solution.front(), b, &RCI_request, ipar, dpar, tmp);

/*---------------------------------------------------------------------------
/* If RCI_request=0, then the solution was found with the required precision
/*---------------------------------------------------------------------------*/
if (RCI_request==0)
statusCode=COMPLETE;

/*---------------------------------------------------------------------------
/* If RCI_request=1, then compute the vector A*tmp[ipar[21]-1]
/* and put the result in vector tmp[ipar[22]-1]
/*---------------------------------------------------------------------------
else if (RCI_request==1) {
mkl_dcsrgemv("N", &ivar, &_data[0], &_row_ptr[0], &_col_ind[0],&tmp[ipar[21]-1], &tmp[ipar[22]-1]);
statusCode=CONTINUE;
}
/*---------------------------------------------------------------------------
/* If RCI_request=2, then do the user-defined stopping test
/* The residual stopping test for the computed solution is performed here
/*---------------------------------------------------------------------------
else if (RCI_request==2) {
/* Request to the dfgmres_get routine to put the solution
/* into b via ipar[12]
/*---------------------------------------------------------------------------
ipar[12]=1;
/* Get the current FGMRES solution in the vector b */
dfgmres_get(&ivar, &solution.front(), b, &RCI_request, ipar, dpar,
tmp, &itercount);
/* Compute the current true residual via MKL (Sparse)
/* BLAS routines */
mkl_dcsrgemv("N", &ivar, &_data[0], &_row_ptr[0], &_col_ind[0], b, residual);
dvar=-1.0E0;
i=1;
daxpy(&ivar, &dvar, b, &i, residual, &i);
dvar=dnrm2(&ivar,residual,&i);
if (dvar<1.0E-10)
statusCode=COMPLETE;
else
statusCode=CONTINUE;
}
/*---------------------------------------------------------------------------
/* 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 ILUT routine
/* via standard MKL Sparse Blas solver routine mkl_dcsrtrsv.
/*---------------------------------------------------------------------------*/
& nbsp; else if (RCI_request==3) {
char cvar='N';
double alpha=1.0;
mkl_dcsrsv(&cvar, &ivar, α, "DLNF ", &_data[0], &_col_ind[0], &_row_ptr[0], &_row_ptr[1], &tmp[ipar[21]-1],&tmp[ipar[22]-1]);
statusCode=CONTINUE;
}
/*---------------------------------------------------------------------------
/* 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
/*---------------------------------------------------------------------------*/
else if (RCI_request==4) {
if (dpar[6]<1.0E-12) statusCode=COMPLETE;
else statusCode=CONTINUE;
}
/*---------------------------------------------------------------------------
/* If RCI_request=anything else, then dfgmres subroutine failed
/* to compute the solution vector: computed_solution
/*---------------------------------------------------------------------------*/
else {
statusCode=FAILED;
}
}

if(statusCode==COMPLETE) {
ipar[12]=0;
dfgmres_get(&ivar, &solution.front(), b, &RCI_request, ipar, dpar, tmp, &itercount);

} else if (statusCode==FAILED)
cout << "The solver has returned the ERROR code " << RCI_request << endl;
else
cout << "The FGMRES solver has FAILED" << endl;

delete [] b;
delete [] tmp;
delete [] residual;
delete [] trvec;

return solution;
}

Any thoughts would be greatly appreciated. Thanks!

-Mark

[Sergey_G_Intel]

Hi Mark,

Thank you for using Intel MKL!

The marked part of your code looks OK.

The convergence problem may come from the big condition number of the matrix (roughly speaking, around 10¹² or larger). Diagonal preconditioner cannot always help to decrease the condition number, thus the method may not converge even with the correct program.

If you feel that you need more assistance with your program, please submit a test case through MKL support page (http://www.intel.com/support/performancetools/libraries/mkl/index.htm), register (if you are not registered yet), and describe the problem and submit the test case. Thank you!

[mdobossy]

Sergey-

Thanks for the info!

The problem was indeed in a different portion of the code. Just was trying to rule out improper pre-conditioning. All is working fine now.

-Mark

[Sergey_G_Intel]

Mark,

I'm glad to hear that your problem is resolved.