# Preconditioner Code Examples

This section contains example codes in FORTRAN 77 and C. For description of the preconditioner routines used in these codes, refer to “Preconditioners (Accelerators) based on Incomplete LU Factorization Techniques” in “Sparse Solver Routines”.

## Fortran Example of Using ILU0 Preconditioner with RCI FGMRES Solver

Example results for non-degenerate system. Upon successful execution of the solver, the output contains the following information:

`-------------------------------------------------------------------`
`The FULLY ADVANCED example RCI FGMRES with ILU0 preconditioner`
`to solve the non-degenerate algebraic system of linear equations`
`-------------------------------------------------------------------`
``
`Some info about the current run of RCI FGMRES method:`
``
`As IPAR(8)=0, the automatic test for the maximal number of iterations will be skipped`
`+++`
`As IPAR(9)=0, the automatic residual test will be skipped`
`+++`
`As IPAR(10)=1, the user-defined stopping test will be requested via RCI_REQUEST=2`
`+++`
`As IPAR(11)=1, the Preconditioned FGMRES iterations will be performed, thus,`
`the preconditioner action will be requested via RCI_REQUEST=3`
`+++`
`As IPAR(12)=0, the automatic test for the norm of the next generated vector is`
`not equal to zero up to rounding and computational errors will be skipped,`
`thus, the user-defined test will be requested via RCI_REQUEST=4`
`+++`
``
`The system has been solved`
``
`The following solution has been obtained`
`COMPUTED_SOLUTION(1)= 0.100E+01`
`COMPUTED_SOLUTION(2)= 0.100E+01`
`COMPUTED_SOLUTION(3)= 0.100E+01`
`COMPUTED_SOLUTION(4)= 0.100E+01`
``
`The expected solution is:`
`EXPECTED_SOLUTION(1)= 0.100E+01`
`EXPECTED_SOLUTION(2)= 0.100E+01`
`EXPECTED_SOLUTION(3)= 0.100E+01`
`EXPECTED_SOLUTION(4)= 0.100E+01`
``
`Number of iterations: 2`
``
`--------------------------------------------------------`
`Fortran example of FGMRES with ILU0 preconditioner`
`has successfully PASSED all stages of computations`
`--------------------------------------------------------`

### Note

If the ILU0 preconditioner is not used, then the same code requires 7 iterations.

## Example. FORTRAN 77 Example to Solve Non-symmetric Non-degenerate Linear System

!*******************************************************************************

! INTEL CONFIDENTIAL

! Copyright(C) 2006-2008 Intel Corporation. All Rights Reserved.

! The source code contained or described herein and all documents related to

! the source code ("Material") are owned by Intel Corporation or its suppliers

! or licensors. Title to the Material remains with Intel Corporation or its

! suppliers and licensors. The Material contains trade secrets and proprietary

! and confidential information of Intel or its suppliers and licensors. The

! Material is protected by worldwide copyright and trade secret laws and

! treaty provisions. No part of the Material may be used, copied, reproduced,

! modified, published, uploaded, posted, transmitted, distributed or disclosed

! in any way without Intel's prior express written permission.

! No license under any patent, copyright, trade secret or other intellectual

! property right is granted to or conferred upon you by disclosure or delivery

! of the Materials, either expressly, by implication, inducement, estoppel or

! otherwise. Any license under such intellectual property rights must be

! express and approved by Intel in writing.

C*******************************************************************************

C Content:

C Intel MKL example of RCI Flexible Generalized Minimal RESidual method with

C ILU0 Preconditioner

C*******************************************************************************

C---------------------------------------------------------------------------

C Example program for solving non-degenerate system of equations.

C Full functionality of RCI FGMRES solver is exploited. Example shows how

C ILU0 preconditioner accelerates the solver by reducing the number of

C iterations.

C---------------------------------------------------------------------------

PROGRAM FGMRES_FULL_FUNCT_F

IMPLICIT NONE

INCLUDE "mkl_rci.fi"

INTEGER N

PARAMETER(N=4)

INTEGER SIZE

PARAMETER (SIZE=128)

C---------------------------------------------------------------------------

C Define arrays for the upper triangle of the coefficient matrix

C Compressed sparse row storage is used for sparse representation

C---------------------------------------------------------------------------

INTEGER IA(5)

DATA IA /1,4,7,10,13/

INTEGER JA(12)

DATA JA /1,2,3,1,2,4,1,3,4,2,3,4/

DOUBLE PRECISION A(12),BILU0(12),TRVEC(N)

DATA A / 4.,-1.,-1.,-1.,4.,-1.,-1.,4.,-1.,-1.,-1.,4./

C---------------------------------------------------------------------------

C Allocate storage for the ?par parameters and the solution/rhs/residual vectors

C---------------------------------------------------------------------------

INTEGER IPAR(SIZE),IERR

DOUBLE PRECISION DPAR(SIZE), TMP(N*(2*N+1)+(N*(N+9))/2+1)

DOUBLE PRECISION EXPECTED_SOLUTION(N)

DATA EXPECTED_SOLUTION /1.0,1.0,1.0,1.0/

DOUBLE PRECISION RHS(N), B(N)

DOUBLE PRECISION COMPUTED_SOLUTION(N)

DOUBLE PRECISION RESIDUAL(N)

INTEGER MATSIZE, INCX, REF_NIT

DOUBLE PRECISION REF_NORM2, NRM2

PARAMETER (MATSIZE=12, INCX=1, REF_NIT=2, REF_NORM2=7.772387d0)

C---------------------------------------------------------------------------

C Some additional variables to use with the RCI (P)FGMRES solver

C---------------------------------------------------------------------------

INTEGER ITERCOUNT

INTEGER RCI_REQUEST, I

DOUBLE PRECISION DVAR

C---------------------------------------------------------------------------

C An external BLAS function is taken from MKL BLAS to use

C with the RCI (P)FGMRES solver

C---------------------------------------------------------------------------

DOUBLE PRECISION DNRM2

EXTERNAL DNRM2

WRITE( *,'(A,A)') '---------------------------------------------',

1 '----------------------'

WRITE(*,'(A,A)') 'The FULLY ADVANCED example RCI FGMRES with',

1 ' ILU0 preconditioner'

WRITE(*,'(A,A)') 'to solve the non-degenerate algebraic system',

1 ' of linear equations'

WRITE( *,'(A,A)') '---------------------------------------------',

1 '----------------------'

C---------------------------------------------------------------------------

C Initialize variables and the right hand side through matrix-vector product

C---------------------------------------------------------------------------

CALL MKL_DCSRGEMV('N', N, A, IA, JA, EXPECTED_SOLUTION, RHS)

C---------------------------------------------------------------------------

C Save the right-hand side in vector B for future use

C---------------------------------------------------------------------------

CALL DCOPY(N, RHS, 1, B, 1)

C---------------------------------------------------------------------------

C Initialize the initial guess

C---------------------------------------------------------------------------

DO I=1,N

COMPUTED_SOLUTION(I)=0.d0

ENDDO

COMPUTED_SOLUTION(1)=100.d0

C---------------------------------------------------------------------------

C Initialize the solver

C---------------------------------------------------------------------------

CALL DFGMRES_INIT(N, COMPUTED_SOLUTION, RHS, RCI_REQUEST, IPAR,

1 DPAR, TMP)

IF (RCI_REQUEST.NE.0) GOTO 999

C---------------------------------------------------------------------------

C Calculate ILU0 preconditioner.

C !ATTENTION!

C DCSRILU0 routine uses some IPAR, DPAR set by DFGMRES_INIT routine.

C Important for DCSRILU0 default entries set by DFGMRES_INIT are

C ipar(2) = 6 - output of error messages to the screen,

C ipar(6) = 1 - allow output of error messages,

C ipar(31)= 0 - abort DCSRILU0 calculations if routine meets zero diagonal element.

C

C If ILU0 is going to be used out of MKL FGMRES context, than the values

C of ipar(2), ipar(6), ipar(31), dpar(31), and dpar(32) should be user

C provided before the DCSRILU0 routine call.

C

C In this example, specific for DCSRILU0 entries are set in turn:

C ipar(31)= 1 - change small diagonal value to that given by dpar(32),

C dpar(31)= 1.D-20 instead of the default value set by DFGMRES_INIT.

C It is a small value to compare a diagonal entry with it.

C dpar(32)= 1.D-16 instead of the default value set by DFGMRES_INIT.

C It is the target value of the diagonal value if it is

C small as compared to dpar(31) and the routine should change

C it rather than abort DCSRILU0 calculations.

C---------------------------------------------------------------------------

IPAR(31)=1

DPAR(31)=1.D-20

DPAR(32)=1.D-16

CALL DCSRILU0(N, A, IA, JA, BILU0, IPAR, DPAR, IERR)

NRM2=DNRM2(MATSIZE, BILU0, INCX)

IF(IERR.ne.0) THEN

WRITE(*,'(A,A,I1)') ' Error after calculation of the',

1 ' preconditioner DCSRILU0',IERR

GOTO 998

ENDIF

C---------------------------------------------------------------------------

C Set the desired parameters:

C do the restart after 2 iterations

C LOGICAL parameters:

C do not do the stopping test for the maximal number of iterations

C do the Preconditioned iterations of FGMRES method

C Set parameter IPAR(11) for preconditioner call. For this example,

C it reduces the number of iterations.

C DOUBLE PRECISION parameters

C set the relative tolerance to 1.0D-3 instead of default value 1.0D-6

C NOTE. Preconditioner may increase the number of iterations for an

C arbitrary case of the system and initial guess and even ruin the

C convergence. It is user's responsibility to use a suitable preconditioner

C and to apply it skillfully.

C---------------------------------------------------------------------------

IPAR(15)=2

IPAR(8)=0

IPAR(11)=1

DPAR(1)=1.0D-3

C---------------------------------------------------------------------------

C Check the correctness and consistency of the newly set parameters

C---------------------------------------------------------------------------

CALL DFGMRES_CHECK(N, COMPUTED_SOLUTION, RHS, RCI_REQUEST,

1 IPAR, DPAR, TMP)

IF (RCI_REQUEST.NE.0) GOTO 999

C---------------------------------------------------------------------------

C Print the info about the RCI FGMRES method

C---------------------------------------------------------------------------

WRITE( *,'(A)') ' '

WRITE( *,'(A,A)') 'Some info about the current run of RCI FGMRES',

1 ' method:'

WRITE( *,'(A)') ' '

IF (IPAR(8).NE.0) THEN

WRITE(*,'(A,I1,A,A)') 'As IPAR(8)=',IPAR(8),', the automatic',

1 ' test for the maximal number of iterations will be performed'

ELSE

WRITE(*,'(A,I1,A,A)') 'As IPAR(8)=',IPAR(8),', the automatic',

1 ' test for the maximal number of iterations will be skipped'

ENDIF

WRITE( *,'(A)') '+++'

IF (IPAR(9).NE.0) THEN

WRITE(*,'(A,I1,A,A)') 'As IPAR(9)=',IPAR(9),', the automatic',

1 ' residual test will be performed'

ELSE

WRITE(*,'(A,I1,A,A)') 'As IPAR(9)=',IPAR(9),', the automatic',

1 ' residual test will be skipped'

ENDIF

WRITE( *,'(A)') '+++'

IF (IPAR(10).NE.0) THEN

WRITE(*,'(A,I1,A,A)') 'As IPAR(10)=',IPAR(10),', the',

1 ' user-defined stopping test will be requested via RCI_REQUEST=2'

ELSE

WRITE(*,'(A,I1,A,A,A)') 'As IPAR(10)=',IPAR(10),', the',

1 ' user-defined stopping test will not be requested, thus,',

1 ' RCI_REQUEST will not take the value 2'

ENDIF

WRITE( *,'(A)') '+++'

IF (IPAR(11).NE.0) THEN

WRITE(*,'(A,I1,A,A)') 'As IPAR(11)=',IPAR(11),', the',

1 ' Preconditioned FGMRES iterations will be performed, thus,'

WRITE(*,'(A,A)') 'the preconditioner action will be requested',

1 ' via RCI_REQUEST=3'

ELSE

WRITE(*,'(A,I1,A,A)') 'As IPAR(11)=',IPAR(11),', the',

1 ' Preconditioned FGMRES iterations will not be performed,'

WRITE( *,'(A)') 'thus, RCI_REQUEST will not take the value 3'

ENDIF

WRITE( *,'(A)') '+++'

IF (IPAR(12).NE.0) THEN

WRITE(*,'(A,I1,A,A)')'As IPAR(12)=',IPAR(12),', the automatic',

1 ' test for the norm of the next generated vector is not'

WRITE( *,'(A,A)') ' equal to zero up to rounding and',

1 ' computational errors will be performed,'

WRITE( *,'(A)') 'thus, RCI_REQUEST will not take the value 4'

ELSE

WRITE(*,'(A,I1,A,A)')'As IPAR(12)=',IPAR(12),', the automatic',

1 ' test for the norm of the next generated vector is'

WRITE(*,'(A,A)') 'not equal to zero up to rounding and',

1 ' computational errors will be skipped,'

WRITE(*,'(A,A)') 'thus, the user-defined test will be requested',

1 ' via RCI_REQUEST=4'

ENDIF

WRITE( *,'(A)') '+++'

C---------------------------------------------------------------------------

C Compute the solution by RCI (P)FGMRES solver with preconditioning

C Reverse Communication starts here

C---------------------------------------------------------------------------

1 CALL DFGMRES(N, COMPUTED_SOLUTION, RHS, RCI_REQUEST, IPAR,

1 DPAR, TMP)

C---------------------------------------------------------------------------

C If RCI_REQUEST=0, then the solution was found with the required precision

C---------------------------------------------------------------------------

IF (RCI_REQUEST.EQ.0) GOTO 3

C---------------------------------------------------------------------------

C If RCI_REQUEST=1, then compute the vector A*TMP(IPAR(22))

C and put the result in vector TMP(IPAR(23))

C---------------------------------------------------------------------------

IF (RCI_REQUEST.EQ.1) THEN

CALL MKL_DCSRGEMV('N',N, A, IA, JA, TMP(IPAR(22)), TMP(IPAR(23)))

GOTO 1

ENDIF

C---------------------------------------------------------------------------

C If RCI_request=2, then do the user-defined stopping test

C The residual stopping test for the computed solution is performed here

C---------------------------------------------------------------------------

C NOTE: from this point vector B(N) is no longer containing the right-hand

C side of the problem! It contains the current FGMRES approximation to the

C solution. If you need to keep the right-hand side, save it in some other

C vector before the call to DFGMRES routine. Here we saved it in vector

C RHS(N). The vector B is used instead of RHS to preserve the original

C right-hand side of the problem and guarantee the proper restart of FGMRES

C method. Vector B will be altered when computing the residual stopping

C criterion!

C---------------------------------------------------------------------------

IF (RCI_REQUEST.EQ.2) THEN

C Request to the DFGMRES_GET routine to put the solution into B(N) via IPAR(13)

IPAR(13)=1

C Get the current FGMRES solution in the vector B(N)

CALL DFGMRES_GET(N, COMPUTED_SOLUTION, B, RCI_REQUEST, IPAR,

1 DPAR, TMP, ITERCOUNT)

C Compute the current true residual via MKL (Sparse) BLAS routines

CALL MKL_DCSRGEMV('N', N, A, IA, JA, B, RESIDUAL)

CALL DAXPY(N, -1.0D0, RHS, 1, RESIDUAL, 1)

DVAR=DNRM2(N, RESIDUAL, 1)

IF (DVAR.LT.1.0E-3) THEN

GOTO 3

ELSE

GOTO 1

ENDIF

ENDIF

C---------------------------------------------------------------------------

C If RCI_REQUEST=3, then apply the preconditioner on the vector

C TMP(IPAR(22)) and put the result in vector TMP(IPAR(23))

C Here is the recommended usage of the result produced by ILU0 routine

C via standard MKL Sparse Blas solver routine mkl_dcsrtrsv.

C---------------------------------------------------------------------------

IF (RCI_REQUEST.EQ.3) THEN

CALL MKL_DCSRTRSV('L','N','U',N,BILU0,IA,JA,TMP(IPAR(22)),TRVEC)

CALL MKL_DCSRTRSV('U','N','N',N,BILU0,IA,JA,TRVEC,TMP(IPAR(23)))

GOTO 1

ENDIF

C---------------------------------------------------------------------------

C If RCI_REQUEST=4, then check if the norm of the next generated vector is

C not zero up to rounding and computational errors. The norm is contained

C in DPAR(7) parameter

C---------------------------------------------------------------------------

IF (RCI_REQUEST.EQ.4) THEN

IF (DPAR(7).LT.1.0D-12) THEN

GOTO 3

ELSE

GOTO 1

ENDIF

C---------------------------------------------------------------------------

C If RCI_REQUEST=anything else, then DFGMRES subroutine failed

C to compute the solution vector: COMPUTED_SOLUTION(N)

C---------------------------------------------------------------------------

ELSE

GOTO 999

ENDIF

C---------------------------------------------------------------------------

C Reverse Communication ends here

C Get the current iteration number and the FGMRES solution. (DO NOT FORGET to

C call DFGMRES_GET routine as computed_solution is still containing

C the initial guess!). Request to DFGMRES_GET to put the solution into

C vector COMPUTED_SOLUTION(N) via IPAR(13)

C---------------------------------------------------------------------------

3 IPAR(13)=0

CALL DFGMRES_GET(N, COMPUTED_SOLUTION, RHS, RCI_REQUEST, IPAR,

1 DPAR, TMP, ITERCOUNT)

C---------------------------------------------------------------------------

C Print solution vector: COMPUTED_SOLUTION(N) and

C the number of iterations: ITERCOUNT

C---------------------------------------------------------------------------

WRITE( *,'(A)') ' '

WRITE( *,'(A)') 'The system has been solved'

WRITE( *,'(A)') ' '

WRITE( *,'(A)') 'The following solution has been obtained:'

DO I=1,N

WRITE(*,'(A18,I1,A2,E10.3)') 'COMPUTED_SOLUTION(',I,')=',

1 COMPUTED_SOLUTION(I)

ENDDO

WRITE( *,'(A)') ' '

WRITE( *,'(A)') 'The expected solution is:'

DO I=1,N

WRITE(*,'(A18,I1,A2,E10.3)') 'EXPECTED_SOLUTION(',I,')=',

1 EXPECTED_SOLUTION(I)

ENDDO

WRITE( *,'(A)') ' '

WRITE( *,'(A,I2)') 'Number of iterations: ',ITERCOUNT

WRITE( *,'(A)') ' '

IF(ITERCOUNT.EQ.REF_NIT.AND.DABS(REF_NORM2-NRM2).LT.1.D-6) THEN

WRITE( *,'(A)') ' '

WRITE( *,'(A)') '------------------------------------------------'

WRITE( *,'(A,A)') 'Fortran example of FGMRES with ILU0',

1 ' preconditioner '

WRITE( *,'(A,A)') 'has successfully PASSED all stages of',

1 ' computations'

WRITE( *,'(A)') '------------------------------------------------'

WRITE( *,'(A)') ' '

GOTO 1000

ELSE

WRITE( *,'(A,A)') 'Probably, the preconditioner was computed',

1 ' incorrectly:'

WRITE( *,'(A,F9.6,A,F9.6)')

1 'Either the preconditioner norm',NRM2,

2 ' differs from the expected norm',REF_NORM2

WRITE( *,'(A,I2,A,I2)'),

1 'and/or the number of iterations',ITERCOUNT,

2 ' differs from the expected number ', REF_NIT

WRITE( *,'(A)') ' '

WRITE( *,'(A,A)') '---------------------------------------------',

1 '----------------------'

WRITE( *,'(A,A)') 'Unfortunately, FGMRES+ILU0 Fortran example',

1 ' has FAILED'

WRITE( *,'(A,A)') '---------------------------------------------',

1 '----------------------'

WRITE( *,'(A)') ' '

GOTO 1000

END IF

999 WRITE( *,'(A,I2)') 'The solver has returned the ERROR code ',

1 RCI_REQUEST

998 WRITE( *,'(A)') ' '

WRITE( *,'(A,A)') '---------------------------------------------',

1 '----------------------'

WRITE( *,'(A,A)') 'Unfortunately, FGMRES+ILU0 Fortran example',

1 ' has FAILED'

WRITE( *,'(A,A)') '---------------------------------------------',

1 '----------------------'

WRITE( *,'(A)') ' '

1000 CONTINUE

END

## C Example of Using ILU0 Preconditioner with RCI FGMRES Solver

C example has the same non-degenerate system as in the previous FORTRAN 77 example. The results are the same up to the notational convention between C and Fortran. Please pay special attention to how it is recommended to handle the differences between the Fortran and C arrays. Specifically, in this example we adjust the addresses for input/result for user-defined operations from `IPAR(22)` to `ipar[21]-1` and from `IPAR(23)` to `ipar[22]-1`, respectively. Specific to ILU0 preconditioner parameter entries are referred to as `ipar[30]`, `dpar[30]`, `dpar[31]`. Upon successful execution of the solver, the following information is printed (up to rounding errors that depend on the computer system):

--------------------------------------------------------

The FULLY ADVANCED example RCI FGMRES with ILU0 preconditioner

to solve the non-degenerate algebraic system of linear equations

--------------------------------------------------------

Some info about the current run of RCI FGMRES method:

As ipar[7]=0, the automatic test for the maximal number of iterations will be

skipped

+++

As ipar[8]=0, the automatic residual test will be skipped

+++

As ipar[9]=1 the user-defined stopping test will be requested via

RCI_request=2

+++

As ipar[10]=1, the Preconditioned FGMRES iterations will be performed, thus,

the preconditioner action will be requested via RCI_request=3

+++

As ipar[11]=0, the automatic test for the norm of the next generated vector is

not equal to zero up to rounding and computational errors will be skipped,

thus, the user-defined test will be requested via RCI_request=4

+++

The system has been solved

The following solution has been obtained

computed_solution[0]=1.000000e+00

computed_solution[1]=1.000000e+00

computed_solution[2]=1.000000e+00

computed_solution[3]=1.000000e+00

The expected solution is:

expected_solution[0]=1.000000e+00

expected_solution[1]=1.000000e+00

expected_solution[2]=1.000000e+00

expected_solution[3]=1.000000e+00

Number of iterations: 2

--------------------------------------------------------

C example of FGMRES with ILU0 preconditioner

has successfully PASSED all stages of computations

--------------------------------------------------------

### Note

If the ILU0 preconditioner is not used, then the same code requires 7 iterations.

## Example. C Example to Solve Non-symmetric Non-degenerate Linear System

/********************************************************************************

/* INTEL CONFIDENTIAL

/* Copyright(C) 2006-2008 Intel Corporation. All Rights Reserved.

/* The source code contained or described herein and all documents related to

/* the source code ("Material") are owned by Intel Corporation or its suppliers

/* or licensors. Title to the Material remains with Intel Corporation or its

/* suppliers and licensors. The Material contains trade secrets and proprietary

/* and confidential information of Intel or its suppliers and licensors. The

/* Material is protected by worldwide copyright and trade secret laws and

/* treaty provisions. No part of the Material may be used, copied, reproduced,

/* modified, published, uploaded, posted, transmitted, distributed or disclosed

/* in any way without Intel's prior express written permission.

/* No license under any patent, copyright, trade secret or other intellectual

/* property right is granted to or conferred upon you by disclosure or delivery

/* of the Materials, either expressly, by implication, inducement, estoppel or

/* otherwise. Any license under such intellectual property rights must be

/* express and approved by Intel in writing.

/*

/********************************************************************************

/* Content:

/* Intel MKL example of RCI Flexible Generalized Minimal RESidual method with

/* ILU0 Preconditioner

/********************************************************************************

/*---------------------------------------------------------------------------

/* Example program for solving non-degenerate system of equations.

/* Full functionality of RCI FGMRES solver is exploited. Example shows how

/* ILU0 preconditioner accelerates the solver by reducing the number of

/* iterations.

/*---------------------------------------------------------------------------*/

#include <stdio.h>

#include "math.h"

#include "mkl_blas.h"

#include "mkl_spblas.h"

#include "mkl_rci.h"

#define N 4

#define size 128

int main(void)

{

/*---------------------------------------------------------------------------

/* Define arrays for the upper triangle of the coefficient matrix

/* Compressed sparse row storage is used for sparse representation

/*---------------------------------------------------------------------------*/

MKL_INT ia[5]={1,4,7,10,13};

MKL_INT ja[12]={1,2,3,1,2,4,1,3,4,2,3,4};

double A[12]={4.,-1.,-1.,-1.,4.,-1.,-1.,4.,-1.,-1.,-1.,4.};

/*---------------------------------------------------------------------------

/* Allocate storage for the ?par parameters and the solution/rhs/residual vectors

/*---------------------------------------------------------------------------*/

MKL_INT ipar[size];

double dpar[size], tmp[N*(2*N+1)+(N*(N+9))/2+1];

double trvec[N],bilu0[12];

double expected_solution[N]={1.0,1.0,1.0,1.0};

double rhs[N], b[N];

double computed_solution[N];

double residual[N];

MKL_INT matsize=12, incx=1, ref_nit=2;

double ref_norm2=7.772387E+0, nrm2;

/*---------------------------------------------------------------------------

/* Some additional variables to use with the RCI (P)FGMRES solver

/*---------------------------------------------------------------------------*/

MKL_INT itercount,ierr=0;

MKL_INT RCI_request, i, ivar;

double dvar;

char cvar,cvar1,cvar2;

printf("--------------------------------------------------------\n");

printf("The FULLY ADVANCED example RCI FGMRES with ILU0 preconditioner\n");

printf("to solve the non-degenerate algebraic system of linear equations\n");

printf("--------------------------------------------------------\n\n");

/*---------------------------------------------------------------------------

/* Initialize variables and the right hand side through matrix-vector product

/*---------------------------------------------------------------------------*/

ivar=N;

cvar='N';

mkl_dcsrgemv(&cvar, &ivar, A, ia, ja, expected_solution, rhs);

/*---------------------------------------------------------------------------

/* Save the right-hand side in vector b for future use

/*---------------------------------------------------------------------------*/

i=1;

dcopy(&ivar, rhs, &i, b, &i);

/*---------------------------------------------------------------------------

/* Initialize the initial guess

/*---------------------------------------------------------------------------*/

for(i=0;i<N;i++)

{

computed_solution[i]=0.0;

}

computed_solution[0]=100.0;

/*---------------------------------------------------------------------------

/* Initialize the solver

/*---------------------------------------------------------------------------*/

dfgmres_init(&ivar, computed_solution, rhs, &RCI_request, ipar, dpar, tmp);

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-20;

dpar[31]=1.E-16;

dcsrilu0(&ivar, A, ia, ja, bilu0, ipar, dpar, &ierr);

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]=1.0E-3;

/*---------------------------------------------------------------------------

/* Check the correctness and consistency of the newly set parameters

/*---------------------------------------------------------------------------*/

dfgmres_check(&ivar, computed_solution, rhs, &RCI_request, ipar, dpar, tmp);

if (RCI_request!=0) goto FAILED;

/*---------------------------------------------------------------------------

/* Print the info about the RCI FGMRES method

/*---------------------------------------------------------------------------*/

printf("Some info about the current run of RCI FGMRES method:\n\n");

if (ipar[7])

{

printf("As ipar[7]=%d, the automatic test for the maximal number of iterations will be\n", ipar[7]);

printf("performed\n");

}

else

{

printf("As ipar[7]=%d, the automatic test for the maximal number of iterations will be\n", ipar[7]);

printf("skipped\n");

}

printf("+++\n");

if (ipar[8])

{

printf("As ipar[8]=%d, the automatic residual test will be performed\n", ipar[8]);

}

else

{

printf("As ipar[8]=%d, the automatic residual test will be skipped\n", ipar[8]);

}

printf("+++\n");

if (ipar[9])

{

printf("As ipar[9]=%d the user-defined stopping test will be requested via\n", ipar[9]);

printf("RCI_request=2\n");

}

else

{

printf("As ipar[9]=%d, the user-defined stopping test will not be requested, thus,\n", ipar[9]);

printf("RCI_request will not take the value 2\n");

}

printf("+++\n");

if (ipar[10])

{

printf("As ipar[10]=%d, the Preconditioned FGMRES iterations will be performed, thus,\n", ipar[10]);

printf("the preconditioner action will be requested via RCI_request=3\n");

}

else

{

printf("As ipar[10]=%d, the Preconditioned FGMRES iterations will not be performed,\n", ipar[10]);

printf("thus, RCI_request will not take the value 3\n");

}

printf("+++\n");

if (ipar[11])

{

printf("As ipar[11]=%d, the automatic test for the norm of the next generated vector is\n", ipar[11]);

printf("not equal to zero up to rounding and computational errors will be performed,\n");

printf("thus, RCI_request will not take the value 4\n");

}

else

{

printf("As ipar[11]=%d, the automatic test for the norm of the next generated vector is\n", ipar[11]);

printf("not equal to zero up to rounding and computational errors will be skipped,\n");

printf("thus, the user-defined test will be requested via RCI_request=4\n");

}

printf("+++\n\n");

/*---------------------------------------------------------------------------

/* Compute the solution by RCI (P)FGMRES solver with preconditioning

/* Reverse Communication starts here

/*---------------------------------------------------------------------------*/

ONE: dfgmres(&ivar, computed_solution, rhs, &RCI_request, ipar, dpar, tmp);

/*---------------------------------------------------------------------------

/* 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, A, 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[N] 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

/* rhs[N]. The vector b is used instead of 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[N]

/* 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[N] */

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, A, ia, ja, b, residual);

dvar=-1.0E0;

i=1;

daxpy(&ivar, &dvar, rhs, &i, residual, &i);

dvar=dnrm2(&ivar,residual,&i);

if (dvar<1.0E-3) 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-12) goto COMPLETE;

else goto ONE;

}

/*---------------------------------------------------------------------------

/* If RCI_request=anything else, then dfgmres subroutine failed

/* to compute the solution vector: computed_solution[N]

/*---------------------------------------------------------------------------*/

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[N] via ipar[12]

/*---------------------------------------------------------------------------*/

COMPLETE: ipar[12]=0;

dfgmres_get(&ivar, computed_solution, rhs, &RCI_request, ipar, dpar, tmp,

&itercount);

/*---------------------------------------------------------------------------

/* Print solution vector: computed_solution[N] and the number

/* of iterations: itercount

/*---------------------------------------------------------------------------*/

printf("The system has been solved \n");

printf("\nThe following solution has been obtained: \n");

for (i=0;i<N;i++)

{

printf("computed_solution[%d]=%e\n",i,computed_solution[i]);

}

printf("\nThe expected solution is: \n");

for (i=0;i<N;i++)

{

printf("expected_solution[%d]=%e\n",i,expected_solution[i]);

}

printf("\nNumber of iterations: %d\n" ,itercount);

printf("\n");

if(itercount==ref_nit && fabs(ref_norm2-nrm2)<1.e-6)

{

printf("--------------------------------------------------------\n");

printf("C example of FGMRES with ILU0 preconditioner \n");

printf("has successfully PASSED all stages of computations\n");

printf("--------------------------------------------------------\n");

return 0;

}

else

{

printf("Probably, the preconditioner was computed incorrectly:\n");

printf("Either the preconditioner norm %e differs from the expected norm

%e\n",nrm2,ref_norm2);

printf("and/or the number of iterations %d differs from the expected number

%d\n",itercount,ref_nit);

printf("-------------------------------------------------------------------\n");

printf("Unfortunately, FGMRES+ILU0 C example has FAILED\n");

printf("-------------------------------------------------------------------\n");

return 0;

}

FAILED:

printf("The solver has returned the ERROR code %d \n", RCI_request);

FAILED1:

printf("-------------------------------------------------------------------\n");

printf("Unfortunately, FGMRES+ILU0 C example has FAILED\n");

printf("-------------------------------------------------------------------\n");

return 0;

}

## Fortran Example of Using ILUT Preconditioner with RCI FGMRES Solver

Example results for non-degenerate system. Upon successful execution of the solver, the output contains the following information:

-------------------------------------------------------------------

The FULLY ADVANCED example RCI FGMRES with ILUT preconditioner

to solve the non-degenerate algebraic system of linear equations

-------------------------------------------------------------------

Some info about the current run of RCI FGMRES method:

As IPAR(8)=0, the automatic test for the maximal number of iterations will be skipped

+++

As IPAR(9)=0, the automatic residual test will be skipped

+++

As IPAR(10)=1, the user-defined stopping test will be requested via RCI_REQUEST=2

+++

As IPAR(11)=1, the Preconditioned FGMRES iterations will be performed, thus,

the preconditioner action will be requested via RCI_REQUEST=3

+++

As IPAR(12)=0, the automatic test for the norm of the next generated vector is

not equal to zero up to rounding and computational errors will be skipped,

thus, the user-defined test will be requested via RCI_REQUEST=4

+++

The system has been solved

The following solution has been obtained:

COMPUTED_SOLUTION(1)= 0.1000020E+01

COMPUTED_SOLUTION(2)= 0.1000007E+01

COMPUTED_SOLUTION(3)= 0.1000006E+01

COMPUTED_SOLUTION(4)= 0.9999908E+00

The expected solution is:

EXPECTED_SOLUTION(1)= 0.1000000E+01

EXPECTED_SOLUTION(2)= 0.1000000E+01

EXPECTED_SOLUTION(3)= 0.1000000E+01

EXPECTED_SOLUTION(4)= 0.1000000E+01

Number of iterations: 4

------------------------------------------------

Fortran example of FGMRES with ILUT preconditioner

has successfully PASSED all stages of computations

------------------------------------------------

### Note

If the ILUT preconditioner is not used, then the same code requires 10 iterations.

## Example. Fortran Example to Solve Non-symmetric Non-degenerate Linear System

!*******************************************************************************

! INTEL CONFIDENTIAL

! Copyright(C) 2007-2008 Intel Corporation. All Rights Reserved.

! The source code contained or described herein and all documents related to

! the source code ("Material") are owned by Intel Corporation or its suppliers

! or licensors. Title to the Material remains with Intel Corporation or its

! suppliers and licensors. The Material contains trade secrets and proprietary

! and confidential information of Intel or its suppliers and licensors. The

! Material is protected by worldwide copyright and trade secret laws and

! treaty provisions. No part of the Material may be used, copied, reproduced,

! modified, published, uploaded, posted, transmitted, distributed or disclosed

! in any way without Intel's prior express written permission.

! No license under any patent, copyright, trade secret or other intellectual

! property right is granted to or conferred upon you by disclosure or delivery

! of the Materials, either expressly, by implication, inducement, estoppel or

! otherwise. Any license under such intellectual property rights must be

! express and approved by Intel in writing.

C*******************************************************************************

C Content:

C Intel MKL example of RCI Flexible Generalized Minimal RESidual method with

C ILUT Preconditioner

C*******************************************************************************

C---------------------------------------------------------------------------

C Example program for solving non-degenerate system of equations.

C Full functionality of RCI FGMRES solver is exploited. Example shows how

C ILUT preconditioner accelerates the solver by reducing the number of

C iterations.

C---------------------------------------------------------------------------

PROGRAM FGMRES_FULL_FUNCT_F

IMPLICIT NONE

INCLUDE "mkl_rci.fi"

INTEGER N

PARAMETER(N=4)

INTEGER SIZE

PARAMETER (SIZE=128)

C---------------------------------------------------------------------------

C Define arrays for the upper triangle of the coefficient matrix

C Compressed sparse row storage is used for sparse representation

C---------------------------------------------------------------------------

INTEGER IA(5),IBILUT(5)

DATA IA /1,4,7,10,13/

INTEGER JA(12),JBILUT(16)

DATA JA /1,2,3,1,2,4,1,3,4,2,3,4/

DOUBLE PRECISION A(12),BILUT(16),TRVEC(N)

DATA A / 4.,-1.,-1.,-1.,4.,-1.,-1.,4.,-1.,-1.,-1.,4./

C---------------------------------------------------------------------------

C Allocate storage for the ?par parameters and the solution/rhs/residual vectors

C---------------------------------------------------------------------------

INTEGER IPAR(SIZE),IERR

DOUBLE PRECISION DPAR(SIZE), TMP(N*(2*N+1)+(N*(N+9))/2+1)

DOUBLE PRECISION EXPECTED_SOLUTION(N)

DATA EXPECTED_SOLUTION /1.0,1.0,1.0,1.0/

DOUBLE PRECISION RHS(N), B(N)

DOUBLE PRECISION COMPUTED_SOLUTION(N)

DOUBLE PRECISION RESIDUAL(N)

INTEGER MATSIZE, INCX, REF_NIT

DOUBLE PRECISION REF_NORM2, NRM2

PARAMETER (MATSIZE=10, INCX=1, REF_NIT=4, REF_NORM2=7.836719d0)

C---------------------------------------------------------------------------

C Some additional variables to use with the RCI (P)FGMRES solver

C---------------------------------------------------------------------------

INTEGER ITERCOUNT

INTEGER RCI_REQUEST, I

DOUBLE PRECISION DVAR

C---------------------------------------------------------------------------

C An external BLAS function is taken from MKL BLAS to use

C with the RCI (P)FGMRES solver

C---------------------------------------------------------------------------

DOUBLE PRECISION DNRM2

EXTERNAL DNRM2

C---------------------------------------------------------------------------

C Some additional variables to use for ILUT preconditioner call

C---------------------------------------------------------------------------

INTEGER MAXFIL

DOUBLE PRECISION TOL

WRITE( *,'(A,A)') '---------------------------------------------',

1 '----------------------'

WRITE(*,'(A,A)') 'The FULLY ADVANCED example RCI FGMRES with',

1 ' ILUT preconditioner'

WRITE(*,'(A,A)') 'to solve the non-degenerate algebraic system',

1 ' of linear equations'

WRITE( *,'(A,A)') '---------------------------------------------',

1 '----------------------'

C---------------------------------------------------------------------------

C Initialize variables and the right hand side through matrix-vector product

C---------------------------------------------------------------------------

CALL MKL_DCSRGEMV('N', N, A, IA, JA, EXPECTED_SOLUTION, RHS)

C---------------------------------------------------------------------------

C Save the right-hand side in vector B for future use

C---------------------------------------------------------------------------

CALL DCOPY(N, RHS, 1, B, 1)

C---------------------------------------------------------------------------

C Initialize the initial guess

C---------------------------------------------------------------------------

DO I=1,N

COMPUTED_SOLUTION(I)=0.d0

ENDDO

COMPUTED_SOLUTION(1)=100.d0

C---------------------------------------------------------------------------

C Initialize the solver

C---------------------------------------------------------------------------

CALL DFGMRES_INIT(N, COMPUTED_SOLUTION, RHS, RCI_REQUEST, IPAR,

1 DPAR, TMP)

IF (RCI_REQUEST.NE.0) GOTO 999

C---------------------------------------------------------------------------

C Calculate ILUT preconditioner.

C !ATTENTION!

C DCSRILUT routine uses some IPAR, DPAR set by DFGMRES_INIT routine.

C Important for DCSRILUT default entries set by DFGMRES_INIT are

C ipar(2) = 6 - output of error messages to the screen,

C ipar(6) = 1 - allow output of error messages,

C ipar(31)= 0 - abort DCSRILUT calculations if routine meets zero diagonal element.

C ipar(7) = 1 - output warn messages if any and continue

C

C If ILUT is going to be used out of MKL FGMRES context, than the values

C of ipar(2), ipar(6), ipar(31), and dpar(31), should be user

C provided before the DCSRILUT routine call.

C

C In this example, specific for DCSRILUT entries are set in turn:

C ipar(31)= 1 - change small diagonal value to that given by dpar(31),

C dpar(31)= 1.D-5 instead of the default value set by DFGMRES_INIT.

C It is the target value of the diagonal value if it is

C small as compared to given tolerance multiplied

C by the matrix row norm and the routine should

C change it rather than abort DCSRILUT calculations.

C---------------------------------------------------------------------------

IPAR(31)=1

DPAR(31)=1.D-5

TOL=1.d-6

MAXFIL=1

CALL DCSRILUT(N, A, IA, JA, BILUT, IBILUT, JBILUT,

& TOL, MAXFIL, IPAR, DPAR, IERR)

NRM2=DNRM2(MATSIZE, BILUT, INCX)

IF(IERR.ne.0) THEN

WRITE(*,'(A,A,I1)') ' Error after calculation of the',

1' preconditioner DCSRILUT',IERR

GOTO 998

ENDIF

C---------------------------------------------------------------------------

C Set the desired parameters:

C do the restart after 2 iterations

C LOGICAL parameters:

C do not do the stopping test for the maximal number of iterations

C do the Preconditioned iterations of FGMRES method

C Set parameter IPAR(11) for preconditioner call.

C For this example, it reduces the number of iterations.

C DOUBLE PRECISION parameters

C set the relative tolerance to 1.0D-3 instead of default value 1.0D-6

C NOTE. Preconditioner may increase the number of iterations for an

C arbitrary case of the system and initial guess and even ruin the

C convergence. It is user's responsibility to use a suitable preconditioner

C and to apply it skillfully.

C---------------------------------------------------------------------------

IPAR(15)=2

IPAR(8)=0

IPAR(11)=1

DPAR(1)=1.0D-3

C---------------------------------------------------------------------------

C Check the correctness and consistency of the newly set parameters

C---------------------------------------------------------------------------

CALL DFGMRES_CHECK(N, COMPUTED_SOLUTION, RHS, RCI_REQUEST,

1 IPAR, DPAR, TMP)

IF (RCI_REQUEST.NE.0) GOTO 999

C---------------------------------------------------------------------------

C Print the info about the RCI FGMRES method

C---------------------------------------------------------------------------

WRITE( *,'(A)') ' '

WRITE( *,'(A,A)') 'Some info about the current run of RCI FGMRES',

1 ' method:'

WRITE( *,'(A)') ' '

IF (IPAR(8).NE.0) THEN

WRITE(*,'(A,I1,A,A)') 'As IPAR(8)=',IPAR(8),', the automatic',

1 ' test for the maximal number of iterations will be performed'

ELSE

WRITE(*,'(A,I1,A,A)') 'As IPAR(8)=',IPAR(8),', the automatic',

1 ' test for the maximal number of iterations will be skipped'

ENDIF

WRITE( *,'(A)') '+++'

IF (IPAR(9).NE.0) THEN

WRITE(*,'(A,I1,A,A)') 'As IPAR(9)=',IPAR(9),', the automatic',

1 ' residual test will be performed'

ELSE

WRITE(*,'(A,I1,A,A)') 'As IPAR(9)=',IPAR(9),', the automatic',

1 ' residual test will be skipped'

ENDIF

WRITE( *,'(A)') '+++'

IF (IPAR(10).NE.0) THEN

WRITE(*,'(A,I1,A,A)') 'As IPAR(10)=',IPAR(10),', the',

1 ' user-defined stopping test will be requested via RCI_REQUEST=2'

ELSE

WRITE(*,'(A,I1,A,A,A)') 'As IPAR(10)=',IPAR(10),', the',

1 ' user-defined stopping test will not be requested, thus,',

1 ' RCI_REQUEST will not take the value 2'

ENDIF

WRITE( *,'(A)') '+++'

IF (IPAR(11).NE.0) THEN

WRITE(*,'(A,I1,A,A)') 'As IPAR(11)=',IPAR(11),', the',

1 ' Preconditioned FGMRES iterations will be performed, thus,'

WRITE(*,'(A,A)') 'the preconditioner action will be requested',

1 ' via RCI_REQUEST=3'

ELSE

WRITE(*,'(A,I1,A,A)') 'As IPAR(11)=',IPAR(11),', the',

1 ' Preconditioned FGMRES iterations will not be performed,'

WRITE( *,'(A)') 'thus, RCI_REQUEST will not take the value 3'

ENDIF

WRITE( *,'(A)') '+++'

IF (IPAR(12).NE.0) THEN

WRITE(*,'(A,I1,A,A)')'As IPAR(12)=',IPAR(12),', the automatic',

1 ' test for the norm of the next generated vector is not'

WRITE( *,'(A,A)') ' equal to zero up to rounding and',

1 ' computational errors will be performed,'

WRITE( *,'(A)') 'thus, RCI_REQUEST will not take the value 4'

ELSE

WRITE(*,'(A,I1,A,A)')'As IPAR(12)=',IPAR(12),', the automatic',

1 ' test for the norm of the next generated vector is'

WRITE(*,'(A,A)') 'not equal to zero up to rounding and',

1 ' computational errors will be skipped,'

WRITE(*,'(A,A)') 'thus, the user-defined test will be requested',

1 ' via RCI_REQUEST=4'

ENDIF

WRITE( *,'(A)') '+++'

C---------------------------------------------------------------------------

C Compute the solution by RCI (P)FGMRES solver with preconditioning

C Reverse Communication starts here

C---------------------------------------------------------------------------

1 CALL DFGMRES(N, COMPUTED_SOLUTION, RHS, RCI_REQUEST, IPAR,

1 DPAR, TMP)

C---------------------------------------------------------------------------

C If RCI_REQUEST=0, then the solution was found with the required precision

C---------------------------------------------------------------------------

IF (RCI_REQUEST.EQ.0) GOTO 3

C---------------------------------------------------------------------------

C If RCI_REQUEST=1, then compute the vector A*TMP(IPAR(22))

C and put the result in vector TMP(IPAR(23))

C---------------------------------------------------------------------------

IF (RCI_REQUEST.EQ.1) THEN

CALL MKL_DCSRGEMV('N',N, A, IA, JA, TMP(IPAR(22)), TMP(IPAR(23)))

GOTO 1

ENDIF

C---------------------------------------------------------------------------

C If RCI_request=2, then do the user-defined stopping test

C The residual stopping test for the computed solution is performed here

C---------------------------------------------------------------------------

C NOTE: from this point vector B(N) is no longer containing the right-hand

C side of the problem! It contains the current FGMRES approximation to the

C solution. If you need to keep the right-hand side, save it in some other

C vector before the call to DFGMRES routine. Here we saved it in vector

C RHS(N). The vector B is used instead of RHS to preserve the original

C right-hand side of the problem and guarantee the proper restart of FGMRES

C method. Vector B will be altered when computing the residual stopping

C criterion!

C---------------------------------------------------------------------------

IF (RCI_REQUEST.EQ.2) THEN

C Request to the DFGMRES_GET routine to put the solution into B(N) via IPAR(13)

IPAR(13)=1

C Get the current FGMRES solution in the vector B(N)

CALL DFGMRES_GET(N, COMPUTED_SOLUTION, B, RCI_REQUEST, IPAR,

1 DPAR, TMP, ITERCOUNT)

C Compute the current true residual via MKL (Sparse) BLAS routines

CALL MKL_DCSRGEMV('N', N, A, IA, JA, B, RESIDUAL)

CALL DAXPY(N, -1.0D0, RHS, 1, RESIDUAL, 1)

DVAR=DNRM2(N, RESIDUAL, 1)

IF (DVAR.LT.1.0E-3) THEN

GOTO 3

ELSE

GOTO 1

ENDIF

ENDIF

C---------------------------------------------------------------------------

C If RCI_REQUEST=3, then apply the preconditioner on the vector

C TMP(IPAR(22)) and put the result in vector TMP(IPAR(23))

C Here is the recommended usage of the result produced by ILUT routine

C via standard MKL Sparse Blas solver routine mkl_dcsrtrsv.

C---------------------------------------------------------------------------

IF (RCI_REQUEST.EQ.3) THEN

CALL MKL_DCSRTRSV('L','N','U',N,BILUT,IBILUT,JBILUT,

& TMP(IPAR(22)),TRVEC)

CALL MKL_DCSRTRSV('U','N','N',N,BILUT,IBILUT,JBILUT,

& TRVEC,TMP(IPAR(23)))

GOTO 1

ENDIF

C---------------------------------------------------------------------------

C If RCI_REQUEST=4, then check if the norm of the next generated vector is

C not zero up to rounding and computational errors. The norm is contained

C in DPAR(7) parameter

C---------------------------------------------------------------------------

IF (RCI_REQUEST.EQ.4) THEN

IF (DPAR(7).LT.1.0D-12) THEN

GOTO 3

ELSE

GOTO 1

ENDIF

C---------------------------------------------------------------------------

C If RCI_REQUEST=anything else, then DFGMRES subroutine failed

C to compute the solution vector: COMPUTED_SOLUTION(N)

C---------------------------------------------------------------------------

ELSE

GOTO 999

ENDIF

C---------------------------------------------------------------------------

C Reverse Communication ends here

C Get the current iteration number and the FGMRES solution. (DO NOT FORGET to

C call DFGMRES_GET routine as computed_solution is still containing

C the initial guess!). Request to DFGMRES_GET to put the solution into

C vector COMPUTED_SOLUTION(N) via IPAR(13)

C---------------------------------------------------------------------------

3 IPAR(13)=0

CALL DFGMRES_GET(N, COMPUTED_SOLUTION, RHS, RCI_REQUEST, IPAR,

1 DPAR, TMP, ITERCOUNT)

C---------------------------------------------------------------------------

C Print solution vector: COMPUTED_SOLUTION(N) and

C the number of iterations: ITERCOUNT

C---------------------------------------------------------------------------

WRITE( *,'(A)') ' '

WRITE( *,'(A)') 'The system has been solved'

WRITE( *,'(A)') ' '

WRITE( *,'(A)') 'The following solution has been obtained:'

DO I=1,N

WRITE(*,'(A18,I1,A2,E14.7)') 'COMPUTED_SOLUTION(',I,')=',

1 COMPUTED_SOLUTION(I)

ENDDO

WRITE( *,'(A)') ' '

WRITE( *,'(A)') 'The expected solution is:'

DO I=1,N

WRITE(*,'(A18,I1,A2,E14.7)') 'EXPECTED_SOLUTION(',I,')=',

1 EXPECTED_SOLUTION(I)

ENDDO

WRITE( *,'(A)') ' '

WRITE( *,'(A,I2)') 'Number of iterations: ',ITERCOUNT

WRITE( *,'(A)') ' '

IF(ITERCOUNT.EQ.REF_NIT.AND.DABS(REF_NORM2-NRM2).LT.1.D-6) THEN

WRITE( *,'(A)') ' '

WRITE( *,'(A)') '------------------------------------------------'

WRITE( *,'(A,A)') 'Fortran example of FGMRES with ILUT',

1 ' preconditioner '

WRITE( *,'(A,A)') 'has successfully PASSED all stages of',

1 ' computations'

WRITE( *,'(A)') '------------------------------------------------'

WRITE( *,'(A)') ' '

GOTO 1000

ELSE

WRITE( *,'(A,A)') 'Probably, the preconditioner was computed',

1 ' incorrectly:'

WRITE( *,'(A,F9.6,A,F9.6)')

1 'Either the preconditioner norm',NRM2,

2 ' differs from the expected norm',REF_NORM2

WRITE( *,'(A,I2,A,I2)'),

1 'and/or the number of iterations ', ITERCOUNT,

2 ' differs from the expected number ', REF_NIT

WRITE( *,'(A)') ' '

WRITE( *,'(A,A)') '---------------------------------------------',

1 '----------------------'

WRITE( *,'(A,A)') 'Unfortunately, FGMRES+ILUT Fortran example',

1 ' has FAILED'

WRITE( *,'(A,A)') '---------------------------------------------',

1 '----------------------'

WRITE( *,'(A)') ' '

GOTO 1000

END IF

999 WRITE( *,'(A,I2)') 'The solver has returned the ERROR code ',

1 RCI_REQUEST

998 WRITE( *,'(A)') ' '

WRITE( *,'(A,A)') '---------------------------------------------',

1 '----------------------'

WRITE( *,'(A,A)') 'Unfortunately, FGMRES+ILUT Fortran example',

1 ' has FAILED'

WRITE( *,'(A,A)') '---------------------------------------------',

1 '----------------------'

WRITE( *,'(A)') ' '

1000 CONTINUE

END

## C Example of Using ILUT Preconditioner with RCI FGMRES Solver

C example has the same non-degenerate system as in the previous FORTRAN 77 example. The results are the same up to the notational convention between C and Fortran. Please pay special attention to how it is recommended to handle the differences between the Fortran and C arrays. Specifically, in this example we adjust the addresses for input/result for user-defined operations from `IPAR(22)` to `ipar[21]-1` and from `IPAR(23)` to `ipar[22]-1`, respectively. Specific to ILUT preconditioner parameter entries are referred to as `ipar[30]`, `dpar[30]`, `dpar[31]`. Upon successful execution of the solver, the following information is printed (up to rounding errors that depend on the computer system):

--------------------------------------------------------

The FULLY ADVANCED example RCI FGMRES with ILUT

preconditioner to solve

the non-degenerate algebraic system of linear equations

--------------------------------------------------------

Some info about the current run of RCI FGMRES method:

As ipar[7]=0, the automatic test for the maximal number of iterations will be

skipped

+++

As ipar[8]=0, the automatic residual test will be skipped

+++

As ipar[9]=1, the user-defined stopping test will be requested via

RCI_request=2

+++

As ipar[10]=1, the Preconditioned FGMRES iterations will be performed, thus,

the preconditioner action will be requested via RCI_request=3

+++

As ipar[11]=0, the automatic test for the norm of the next generated vector is

not equal to zero up to rounding and computational errors will be skipped,

thus, the user-defined test will be requested via RCI_request=4

+++

The system has been solved

The following solution has been obtained:

computed_solution[0]=1.000020e+00

computed_solution[1]=1.000007e+00

computed_solution[2]=1.000006e+00

computed_solution[3]=9.999908e-01

The expected solution is:

expected_solution[0]=1.000000e+00

expected_solution[1]=1.000000e+00

expected_solution[2]=1.000000e+00

expected_solution[3]=1.000000e+00

Number of iterations: 4

--------------------------------------------------------

C example of FGMRES with ILUT preconditioner

has successfully PASSED all stages of computations

--------------------------------------------------------

### Note

If the ILUT preconditioner is not used, then the same code requires 10 iterations.

## Example. C Example to Solve Non-symmetric Non-degenerate Linear System

/********************************************************************************

/* INTEL CONFIDENTIAL

/* Copyright(C) 2007-2008 Intel Corporation. All Rights Reserved.

/* The source code contained or described herein and all documents related to

/* the source code ("Material") are owned by Intel Corporation or its suppliers

/* or licensors. Title to the Material remains with Intel Corporation or its

/* suppliers and licensors. The Material contains trade secrets and proprietary

/* and confidential information of Intel or its suppliers and licensors. The

/* Material is protected by worldwide copyright and trade secret laws and

/* treaty provisions. No part of the Material may be used, copied, reproduced,

/* modified, published, uploaded, posted, transmitted, distributed or disclosed

/* in any way without Intel's prior express written permission.

/* No license under any patent, copyright, trade secret or other intellectual

/* property right is granted to or conferred upon you by disclosure or delivery

/* of the Materials, either expressly, by implication, inducement, estoppel or

/* otherwise. Any license under such intellectual property rights must be

/* express and approved by Intel in writing.

/*

/********************************************************************************

/* Content:

/* Intel MKL example of RCI Flexible Generalized Minimal RESidual method with

/* ILUT Preconditioner

/********************************************************************************

/*---------------------------------------------------------------------------

/* Example program for solving non-degenerate system of equations.

/* Full functionality of RCI FGMRES solver is exploited. Example shows how

/* ILUT preconditioner accelerates the solver by reducing the number of

/* iterations.

/*---------------------------------------------------------------------------*/

#include <stdio.h>

#include "math.h"

#include "mkl_blas.h"

#include "mkl_spblas.h"

#include "mkl_rci.h"

#define N 4

#define size 128

int main(void)

{

/*---------------------------------------------------------------------------

/* Define arrays for the upper triangle of the coefficient matrix

/* Compressed sparse row storage is used for sparse representation

/*---------------------------------------------------------------------------*/

MKL_INT ia[5]={1,4,7,10,13}, ibilut[5];

MKL_INT ja[12]={1,2,3,1,2,4,1,3,4,2,3,4}, jbilut[16];

double A[12]={4.,-1.,-1.,-1.,4.,-1.,-1.,4.,-1.,-1.,-1.,4.};

/*---------------------------------------------------------------------------

/* Allocate storage for the ?par parameters and the solution/rhs/residual vectors

/*---------------------------------------------------------------------------*/

MKL_INT ipar[size];

double dpar[size], tmp[N*(2*N+1)+(N*(N+9))/2+1];

double trvec[N],bilut[12];

double expected_solution[N]={1.0,1.0,1.0,1.0};

double rhs[N], b[N];

double computed_solution[N];

double residual[N];

/*---------------------------------------------------------------------------

/* Some additional variables to use with the RCI (P)FGMRES solver

/*---------------------------------------------------------------------------*/

MKL_INT itercount,ierr=0;

MKL_INT RCI_request, i, ivar;

double dvar;

char cvar,cvar1,cvar2;

MKL_INT matsize=10, incx=1, ref_nit=4;

double ref_norm2=7.836719E+0, nrm2;

/*---------------------------------------------------------------------------

/* Some additional variables to use for ILUT preconditioner call

/*---------------------------------------------------------------------------*/

MKL_INT maxfil;

double tol;

printf("--------------------------------------------------------\n");

printf("The FULLY ADVANCED example RCI FGMRES with ILUT\n");

printf("preconditioner to solve \n");

printf("the non-degenerate algebraic system of linear equations\n");

printf("--------------------------------------------------------\n\n");

/*---------------------------------------------------------------------------

/* Initialize variables and the right hand side through matrix-vector product

/*---------------------------------------------------------------------------*/

ivar=N;

cvar='N';

mkl_dcsrgemv(&cvar, &ivar, A, ia, ja, expected_solution, rhs);

/*---------------------------------------------------------------------------

/* Save the right-hand side in vector b for future use

/*---------------------------------------------------------------------------*/

i=1;

dcopy(&ivar, rhs, &i, b, &i);

/*---------------------------------------------------------------------------

/* Initialize the initial guess

/*---------------------------------------------------------------------------*/

for(i=0;i<N;i++)

{

computed_solution[i]=0.0;

}

computed_solution[0]=100.0;

/*---------------------------------------------------------------------------

/* Initialize the solver

/*---------------------------------------------------------------------------*/

dfgmres_init(&ivar, computed_solution, rhs, &RCI_request, ipar, dpar, tmp);

if (RCI_request!=0) goto FAILED;

/*---------------------------------------------------------------------------

/* Calculate ILUT preconditioner.

/* !ATTENTION!

/* DCSRILUT routine uses some IPAR, DPAR set by DFGMRES_INIT routine.

/* Important for DCSRILUT 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[6] = 1 - output warn messages if any and continue

/* ipar[30]= 0 - abort DCSRILUT calculations if routine meets zero diagonal element.

/*

/* If ILUT is going to be used out of MKL FGMRES context, than the values

/* of ipar[1], ipar[5], ipar[6], ipar[30], and dpar[30] should be user

/* provided before the DCSRILUT routine call.

/*

/* In this example, specific for DCSRILUT entries are set in turn:

/* ipar[30]= 1 - change small diagonal value to that given by dpar[31],

/* dpar[30]= 1.E-5 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 the given tolerance multiplied

/* by the matrix row norm and the routine should

/* change it rather than abort DCSRILUT calculations.

/*---------------------------------------------------------------------------*/

ipar[30]=1;

dpar[30]=1.E-5;

tol=1.E-6;

maxfil=1;

dcsrilut(&ivar, A, ia, ja, bilut, ibilut, jbilut, &tol, &maxfil, ipar, dpar, &ierr);

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

if (ierr!=0)

{

printf("Preconditioner dcsrilut 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]=1.0E-3;

/*---------------------------------------------------------------------------

/* Check the correctness and consistency of the newly set parameters

/*---------------------------------------------------------------------------*/

dfgmres_check(&ivar, computed_solution, rhs, &RCI_request, ipar, dpar, tmp);

if (RCI_request!=0) goto FAILED;

/*---------------------------------------------------------------------------

/* Print the info about the RCI FGMRES method

/*---------------------------------------------------------------------------*/

printf("Some info about the current run of RCI FGMRES method:\n\n");

if (ipar[7])

{

printf("As ipar[7]=%d, the automatic test for the maximal number of

iterations will be\n", ipar[7]);

printf("performed\n");

}

else

{

printf("As ipar[7]=%d, the automatic test for the maximal number of

iterations will be\n", ipar[7]);

printf("skipped\n");

}

printf("+++\n");

if (ipar[8])

{

printf("As ipar[8]=%d, the automatic residual test will be

performed\n", ipar[8]);

}

else

{

printf("As ipar[8]=%d, the automatic residual test will be skipped\n",

ipar[8]);

}

printf("+++\n");

if (ipar[9])

{

printf("As ipar[9]=%d, the user-defined stopping test will be requested

via\n", ipar[9]);

printf("RCI_request=2\n");

}

else

{

printf("As ipar[9]=%d, the user-defined stopping test will not be

requested, thus,\n", ipar[9]);

printf("RCI_request will not take the value 2\n");

}

printf("+++\n");

if (ipar[10])

{

printf("As ipar[10]=%d, the Preconditioned FGMRES iterations will

be performed, thus,\n", ipar[10]);

printf("the preconditioner action will be requested via

RCI_request=3\n");

}

else

{

printf("As ipar[10]=%d, the Preconditioned FGMRES iterations will

not be performed,\n", ipar[10]);

printf("thus, RCI_request will not take the value 3\n");

}

printf("+++\n");

if (ipar[11])

{

printf("As ipar[11]=%d, the automatic test for the norm of the next

generated vector is\n", ipar[11]);

printf("not equal to zero up to rounding and computational errors

will be performed,\n");

printf("thus, RCI_request will not take the value 4\n");

}

else

{

printf("As ipar[11]=%d, the automatic test for the norm of the next

generated vector is\n", ipar[11]);

printf("not equal to zero up to rounding and computational errors

will be skipped,\n");

printf("thus, the user-defined test will be requested via

RCI_request=4\n");

}

printf("+++\n\n");

/*---------------------------------------------------------------------------

/* Compute the solution by RCI (P)FGMRES solver with preconditioning

/* Reverse Communication starts here

/*---------------------------------------------------------------------------*/

ONE: dfgmres(&ivar, computed_solution, rhs, &RCI_request, ipar, dpar, tmp);

/*---------------------------------------------------------------------------

/* 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, A, 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[N] 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

/* rhs[N]. The vector b is used instead of 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[N] 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[N] */

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, A, ia, ja, b, residual);

dvar=-1.0E0;

i=1;

daxpy(&ivar, &dvar, rhs, &i, residual, &i);

dvar=dnrm2(&ivar,residual,&i);

if (dvar<1.0E-3) 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 ILUT 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,bilut,ibilut,jbilut,&tmp[ipar[21]-1],trvec);

cvar1='U';

cvar='N';

cvar2='N';

mkl_dcsrtrsv(&cvar1,&cvar,&cvar2,&ivar,bilut,ibilut,jbilut,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-12) goto COMPLETE;

else goto ONE;

}

/*---------------------------------------------------------------------------

/* If RCI_request=anything else, then dfgmres subroutine failed

/* to compute the solution vector: computed_solution[N]

/*---------------------------------------------------------------------------*/

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[N] via ipar[12]

/*---------------------------------------------------------------------------*/

COMPLETE: ipar[12]=0;

dfgmres_get(&ivar, computed_solution, rhs, &RCI_request, ipar, dpar, tmp,

&itercount);

/*---------------------------------------------------------------------------

/* Print solution vector: computed_solution[N] and the number

/* of iterations: itercount

/*---------------------------------------------------------------------------*/

printf("The system has been solved \n");

printf("\nThe following solution has been obtained: \n");

for (i=0;i<N;i++)

{

printf("computed_solution[%d]=",i);

printf("%e\n",computed_solution[i]);

}

printf("\nThe expected solution is: \n");

for (i=0;i<N;i++)

{

printf("expected_solution[%d]=",i);

printf("%e\n",expected_solution[i]);

}

printf("\nNumber of iterations: %d\n" ,itercount);

printf("\n");

if(itercount==ref_nit && fabs(ref_norm2-nrm2)<1.e-6)

{

printf("--------------------------------------------------------\n");

printf("C example of FGMRES with ILUT preconditioner \n");

printf("has successfully PASSED all stages of computations\n");

printf("--------------------------------------------------------\n");

return 0;

}

else

{

printf("Probably, the preconditioner was computed incorrectly:\n");

printf("Either the preconditioner norm %e differs from the expected norm

%e\n",nrm2,ref_norm2);

printf("and/or the number of iterations %d differs from the expected number

%d\n",itercount,ref_nit);

printf("-------------------------------------------------------------------\n");

printf("Unfortunately, FGMRES+ILUT C example has FAILED\n");

printf("-------------------------------------------------------------------\n");

return 0;

}

FAILED:

printf("The solver has returned the ERROR code %d", RCI_request);

FAILED1:

printf("-------------------------------------------------------------------\n");

printf("Unfortunately, FGMRES+ILUT C example has FAILED\n");

printf("-------------------------------------------------------------------\n");

return 0;

}

Parent topic: Code Examples