/* Test file for MKL CG RCI method */
/*
   Calling:
   
   export LD_LIBRARY_PATH=/Components64/MKL/11.2.0/Linux-x86-64/vendor
   gcc -o test-mkl test-mkl-cg-rci-ilu.c -I /Components64/MKL/11.2.0/Linux-x86-64/vendor -L /Components64/MKL/11.2.0/Linux-x86-64/vendor -lmkl_def -lmkl_intel_ilp64 -lmkl_core -lmkl_intel_thread -liomp5 -lm
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#define MKL_ILP64
#include "mkl.h"

#define RCI_PAR_LEN 128

static int MKL_CG_method(MKL_INT n, double *x, double *b, double *ad, 
			 MKL_INT *ia, MKL_INT *ja, MKL_INT a_len, 
			 MKL_INT *itercount)
{
   int i, fail = 0;
   MKL_INT len_ilut;
   int err;
   MKL_INT RCI_request;
   MKL_INT ilu_err;
   MKL_INT ipar[RCI_PAR_LEN];
   MKL_INT ilut_maxfil, *iilut, *jilut;
   MKL_INT maxit = 1000;
   double dpar[RCI_PAR_LEN];
   double *tmp, *trvec;
   double tol = 1e-4, ilut_tol = 1e-12;
   double *ilu;
   double alpha = 1.0, beta = 0.0;
   char matdescra[4] = {'G',' ',' ','F'};
   int rel_resQ = 1;
   int use_precond = 1;
   int use_ilut = 0;
   
   if (use_precond) {
      trvec = malloc(sizeof(double)*n);

      if (use_ilut) {
	 iilut = malloc(sizeof(MKL_INT)*(n+1));
	 ilut_maxfil = (MKL_INT) (((double)a_len)/((double)n));
	 len_ilut = (2*ilut_maxfil+1)*n-ilut_maxfil*(ilut_maxfil+1)+1;
	 ilu = malloc(sizeof(double)*len_ilut);
	 jilut = malloc(sizeof(MKL_INT)*len_ilut);
      } else {
	 ilu = malloc(sizeof(double)*a_len);
      }
   }
   
   /* create 4 packed arrays corresponding to tmp[n,1;4] required for 
      A.x and preconditioning. */
   tmp = (double*) malloc(sizeof(double)*4*n);

   /* initialize CG */
   dcg_init(&n, x, b, &RCI_request, ipar, dpar, tmp);
   if (RCI_request == -10000) {
      fail = -1;
      goto cleanup;
   }
   /* change parameters */
   ipar[4] = maxit;  /* maximum number of iterations */
   ipar[5] = 0; /* do not output error messages */
   ipar[6] = 0; /* do not output warning messages */
   ipar[7] = 1; /* stopping test for the maximum # of iterations */
   ipar[8] = 1; /* perform residual stopping test */
   ipar[9] = 0; /* do not perform the user-defined stopping test */
   ipar[10] = use_precond;  /* perform preconditioning if non-NULL */
   ipar[30] = 1; /* zero diagonal element of ILU is replaced */
   
   if (rel_resQ) {
      dpar[0] = tol;  /* relative tolerance */
      dpar[1] = 0.0;  /* absolute tolerance */
   } else {
      dpar[0] = 0.0;  /* relative tolerance */
      dpar[1] = tol;  /* absolute tolerance */      
   }
   /* def. small value for ILU diagonal elements */
   if (use_precond && use_ilut)
      dpar[30] = ilut_tol; 
   else
      dpar[30] = 1e-14;
   
   dcg_check(&n, x, b, &RCI_request, ipar, dpar, tmp);
   if (RCI_request == -1100) {
      fail = -1;
      goto cleanup;
   }
   
   /* initialize ILU */
   if (use_precond && !use_ilut) {
      dcsrilu0(&n, ad, ia, (MKL_INT*)ja, ilu, ipar, dpar, &ilu_err);
      if (ilu_err < 0) {
	 fail = -1;
	 goto cleanup;
      }
   } else if (use_precond && use_ilut) {
      dcsrilut(&n, ad, ia, (MKL_INT*)ja, ilu, iilut, jilut, 
	       &ilut_tol, &ilut_maxfil, ipar, dpar, &ilu_err);
      if (ilu_err < 0) {
	 fail = -1;
	 goto cleanup;
      }
      if (ilu_err > 0) {
	 fail = 1;
      }
      
   } else { /* no other ILU method is supported */
      fail = -2;
      goto cleanup;
   }
   
   do { /* do-while loop */
      dcg(&n, x, b, &RCI_request, ipar, dpar, tmp);
      
      if (RCI_request < 0) {
	 /* error */
	 if (RCI_request == -1) {
	    fail = 1;
	    break;
	 } else if (RCI_request == -2) {
	    fail = 2;
	    break;
	 } else {
	    fail = 1;
	    goto cleanup;
	 }
      } 
      switch(RCI_request) {
	 case 1:
	    /* multiply A by tmp[0]; save result in tmp[n] */
	    mkl_dcsrgemv("N", &n, ad, ia, ja, &tmp[0], &tmp[n]);
	    break;
	 case 2:
	    /* stopping test is currently not implemented; */
	    break;
	 case 3:
	    /* apply the preconditioner to tmp[2*n]; save the result in tmp[3*n] */
	    if (use_precond && !use_ilut) {
	       mkl_dcsrtrsv("L", "N", "U", &n, ilu, ia, (MKL_INT*)ja,
			    &tmp[2*n], trvec);
	       mkl_dcsrtrsv("U", "N", "N", &n, ilu, ia, (MKL_INT*)ja,
			    trvec, &tmp[3*n]);
	    } else if (use_precond && use_ilut) { 
	       mkl_dcsrtrsv("L", "N", "U", &n, ilu, iilut, jilut,
			    &tmp[2*n], trvec);
	       mkl_dcsrtrsv("U", "N", "N", &n, ilu, iilut, jilut,
			    trvec, &tmp[3*n]);
	    } else { /* no ILU */
	       memcpy(&tmp[3*n], &tmp[2*n], n*sizeof(double));
	    }
	    break;
      }
      
   } while(RCI_request != 0);
   
   dcg_get(&n, x, b, &RCI_request, ipar, dpar, tmp, itercount);
   
   printf("tmp:\n");
   for (i = 0; i < n; i++)
      printf("tmp[%d]=%g, ", i, tmp[i]);
   printf("\n");
   
 cleanup:
   free(tmp);
   if (use_precond) {
      free(ilu);
      free(trvec);
      if (use_ilut) {
	 free(iilut);
	 free(jilut);
      }
   }
   
   return fail;
}

int main()
{
   int i, fail;
   MKL_INT n = 10;
   MKL_INT a_len = 28;
   MKL_INT niter = 0;
   double x[10] = {0.0}; /* all zeros */
   double b[10] = {.0, .0, .0, .0, .0, .0, .0, .0, .0, 1.0};
   double ad[28] = {-2., 1., 1., -2., 1., 1., -2., 1., 1., -2., 1., 1., -2., 1.,
      1., -2., 1., 1., -2., 1., 1., -2., 1., 1., -2., 1., 1., -2.};
   MKL_INT ia[11] = {1, 3, 6, 9, 12, 15, 18, 21, 24, 27, 29};
   MKL_INT ja[28] = {1, 2, 1, 2, 3, 2, 3, 4, 3, 4, 5, 4, 5, 6, 5, 6, 7, 6,
      7, 8, 7, 8, 9, 8, 9, 10, 9, 10};
   
   fail = MKL_CG_method(n, x, b, ad, ia, ja, a_len, &niter);
   if (fail != 0) {
      printf("fail = %d\n", fail);
      return 1;
   }
   printf("# of iterations = %d\n", niter);
   printf("solution x:\n");
   for (i = 0; i < n; i++) {
      printf("x[%d]=%g, ", i, x[i]);
      if (n == 4) printf("\n");
   }
   printf ("\n");
   
   return 0;
}