RCI - ISS

marcias_martinez · ‎05-01-2007

Hi Guys,

I have been looking at this simple example that I obtained from Intel Support on RCI - ISS. The issue that I am having is the following:

1. We are trying to solve using conjugate gradient solver Kx=b.

K is represented in the program by vector a(i). x(i) is the expected solution and rhs is equal to b.

In Conjugate gradient solver one needs to take an initial guess (double expected_sol[8]={1.0e0, 0.0e0, 1.0e0, 0.0e0, 1.0e0, 0.0e0, 1.0e0, 0.0e0};).

The program some how runs through a series of iterations and congerce to this initial solution double expected_sol[8]={1.0e0, 0.0e0, 1.0e0, 0.0e0, 1.0e0, 0.0e0, 1.0e0, 0.0e0};. So I try chaging K or b and it still provides me with the same initial guess.

Any ideas what is going on or what I am missing?
Thanks,

Marcias

*******************************************************************************

! Content: Intel MKL RCI (P)CG C example

!

!******************************************************************************/

#include

#define

length 128

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

/* Example program for solving symmetric positive definite system */

/* of equations. Full case: full functionality of RCI (P)CG is used. */

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

extern

"C" {

void mkl_dcsrsymv(char *uplo, int *m, double *a, int *ia, int *ja, double *x, double *y);

void DCG_INIT(int *n, double *x,double *b, int *rci_request, int *ipar, double *dpar, double *tmp);

void DCG_CHECK(int *n, double *x,double *b, int *rci_request, int *ipar, double *dpar, double *tmp);

void DCG(int *n, double *x,double *b, int *rci_request, int *ipar, double *dpar, double *tmp);

void DCG_GET(int *n, double *x, double *b, int *rci_request, int *ipar, double *dpar, double *tmp, int *itercount);

void mkl_dcsrsv(char *transa, int *m, double *alpha, char *matdescra, double *val, int *indx, int *pntrb, int *pntre, double *x, double *y);

void daxpy(int *n,double *alpha,double *x,int *incx,double *y,int *incy);

double dnrm2(int *n,double *x,int *incx);

}

int

main(void)

{

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

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

/* preconditioner as well as an array for rhs vector */

/* Compressed sparse row storage is used for sparse representation */

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

int RCI_request, itercount, i;

double rhs[8];

/* Fill all arrays containing matrix data.*/

int ia[9]={1,5,8,10,12,15,17,18,19};

int ja[18]={1, 3, 6,7,

2,3, 5,

3, 8,

4, 7,

5,6,7,

6, 8,

7,

8};

double a[18]={ 7.0e0, 1.0e0, 2.0e0, 7.0e0,

-4.0e0, 8.0e0, 2.0e0,

1.0e0, 5.0e0,

7.0e0, 9.0e0,

5.0e0, 1.0e0, 5.0e0,

-1.0e0, 5.0e0,

11.0e0,

5.0e0};

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

/* Allocate storage for the solver ?par and the initial solution vector */

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

int ipar[length];

double dpar[length],tmp[4*8];

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

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

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

char matdes[3], U='U', NT='N';

int N=8, ione=1;

double mdone=-1.0e0, one=1.0e0, solution[8];

double expected_sol[8]={1.0e0, 0.0e0, 1.0e0, 0.0e0, 1.0e0, 0.0e0, 1.0e0, 0.0e0};

double euclidean_norm, temp[8];

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

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

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

mkl_dcsrsymv(&U, &N, a, ia, ja, expected_sol, rhs);

/*---------------------------------------------------------------------------*/< /P>

/* Initialize the initial guess */

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

for(i=0;i=0.0e0;

matdes[0]='D';

matdes[1]='L';

matdes[2]='N';

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

/* Initialize the solver */

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

DCG_INIT(&N, solution,rhs, &RCI_request,ipar,dpar,tmp);

if (RCI_request!=0) goto FAILED;

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

/* Set the desired parameters: */

/* INTEGER parameters: */

/* set the maximal number of iterations to 100 */

/* LOGICAL parameters: */

/* run the Preconditioned version of RCI (P)CG with preconditioner C_inverse */

/* DOUBLE PRECISION parameters */

/* - */

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

ipar[4]=100;

ipar[10]=1;

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

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

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

DCG_CHECK(&N,solution,rhs,&RCI_request,ipar,dpar,tmp);

if (RCI_request!=0) goto FAILED;

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

/* Compute the solution by RCI (P)CG solver */

/* Reverse Communications starts here */

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

RCI: DCG(&N,solution,rhs,&RCI_request,ipar,dpar,tmp);

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

/* If RCI_request=0, then the solution was found according to the reque sted */

/* stopping tests. In this case, this means that it was found after 100 */

/* iterations. */

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

if (RCI_request==0) goto STOP_RCI;

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

/* If RCI_request=1, then compute the vector A*TMP[0] */

/* and put the result in vector TMP */

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

if (RCI_request==1)

{

mkl_dcsrsymv(&U, &N, a, ia, ja, tmp, &tmp);

goto RCI;

}

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

/* If RCI_request=2, then do the user-defined stopping test: compute the */

/* Euclidean norm of the actual residual using MKL routines and check if */

/* it is less than 1.0e-8 */

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

if (RCI_request==2)

{

mkl_dcsrsymv(&U, &N, a, ia, ja, solution, temp);

daxpy(&N,&mdone,rhs,&ione,temp,&ione);

euclidean_norm = dnrm2(&N,temp,&ione);

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

/* The solution has not been found yet according to the user-defined stopping*/

/* test. Continue RCI (P)CG iterations. */

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

if (euclidean_norm>1.0e-8) goto RCI;

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

/* The solution has been found according to the user-defined stopping test */

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

else goto STOP_RCI;

}

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

/* If RCI_request=3, then compute apply the preconditioner matrix C_inverse */

/* on vector TMP[2*N] and put the result in vector TMP[3*N] */

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

if (RCI_request==3)

{

mkl_dcsrsv(&NT, &N, &one, matdes, a, ja, ia, &ia[2], &tmp[2*N], &tmp[3*N]);

goto RCI;

}

else

{

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

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

/* to compute the solution vector: solution */

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

goto FAILED;

}

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

/* Reverse Communication ends here */

/* Get the current iteration number */

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

STOP_RCI: DCG_GET(&N,solution,rhs,&RCI_request,ipar,dpar,tmp,&itercount);

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

/* Print solution vector: solution and number of iterations: itercount */

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

printf(" The system is successfully solved ");

printf(" The following solution obtained ");

for(i=0;i);

printf(" Expected solution: ");

for(i=0;i);

printf(" Number of iterations: %d",itercount);

goto JUMP;

FAILED: printf(" Solver returned error code %d ", RCI_request);

return -1;

JUMP:

return 0;

}

Sergey_G_Intel · ‎05-24-2007

Hi Marcias,

The initial guess is actually in the vector solution[]. According to the code,

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

/* Initialize the initial guess */

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

for(i=0;i=0.0e0;

the initial guess is vector of zeros. As far as I can see from your message (excuse meif I'm mistaken), you are looking at theexpected solution vector expected_solution[]. This vector contains the exact solution to the problem and it is created in the line

double expected_sol[8]={1.0e0, 0.0e0, 1.0e0, 0.0e0, 1.0e0, 0.0e0, 1.0e0, 0.0e0};

It does not change through the code at all independently on the changes you made in K and b.

If you change matrix K and/or vector b, you will see changes in the computed solution vector solution[]. You may alsochange the initial guess by changing the loop above.

With kind regards,

Sergey