/***********************************************************
 * test_scal.c                                             *
 *  Solve the linear system Ax=b with scalapack and checks *
 *  the solution.                                          *
 *  The entries of A are random numbers and the solution   *
 *  is the vector with all components equal to one.        *
 *                                                         *
 ***********************************************************/

#include <stdio.h>
#include <math.h>
#include <mpi.h>
#include <stdlib.h>
#include <time.h>
#include <mkl_scalapack.h>

int bl=4, nb=2, zero=0, one=1;

/** Header files for Blacs and Pblas are not available from Intel Mkl */
void Cblacs_pinfo(int*, int*);
void Cblacs_get(int, int, int*);
void Cblacs_gridinit(int*, char*, int, int);
void Cblacs_gridinfo(int, int*, int*, int*, int*);
void Cblacs_barrier(int , char*);
void Cblacs_gridexit(int);
void descinit_(int*, int*, int*, int*, int*, int*, int*, int*, int*, int*);  
  
void pdgemv_(char*,int*,int*,double*,double*,int*,int*,int*,double*,
	int*,int*,int*,int*,double*,double*,int*,int*,int*,int*);
void pdaxpy_(int*, double*, double*, int*, int*, int*, int*, double*, int*, int*, int*, int*);
void pdnrm2_(int*, double*, double*, int*, int*, int*, int*);

/** Local functions */

void block_A (int bl, double* a)
/* Fills a square matrix block of order bl with random entries */
{  int i,j;

   for (i=0; i<bl; i++)
      for (j=0; j<bl; j++ ) a[bl*i+j] = (rand()/((double)RAND_MAX+1));
}; 

void block_x (int bl, double* x)
/* Fills the vector block of length bl with ones*/
{  int j;
   double one=1;
   for (j=0; j<bl; j++ ) x[j]=one;
}; 

int main(int argc, char* argv[]) 
{   int ctxt, n=nb*bl, nprow=nb, npcol=nb, myrow, mycol;  /** BLACS stuff: Context, grid */
    int descA[9], descx[9], descy[9], info;  /** BLACS stuff: descriptors */
    
    double A[bl*bl], x[bl], y[bl]; /* local size of the arrays */
    int ipiv[bl+nb];
    
    double alpha=1, beta=0, norm, gflops;
    int  iam,nprocs,i,j,ib,jb,iad,jad;
	
    clock_t clock0, clock1;

    
    /* Start MPI */
    
    MPI_Init(&argc, &argv);
    
    /* Initialize blacs and the process grid */
    Cblacs_pinfo(&iam,&nprocs);
    if (iam==0 && nprocs<nb*nb)
    {  printf("\n Number of processors should be at least %d\n",nb*nb);
       exit(1);
    }
    

/**   define the process grid */

    Cblacs_get(0,0,&ctxt);
    Cblacs_gridinit(&ctxt,"Row-major",nprow,npcol);
    
    /* Find who I am in the grid */
    Cblacs_gridinfo(ctxt,&nprow,&npcol,&myrow,&mycol);

    /* Initialize BLACS descriptors */
    descinit_(descA,&n,&n,&nb,&nb,&zero,&zero,&ctxt,&bl,&info);
    descinit_(descx,&one,&n,&one,&nb,&zero,&zero,&ctxt,&one,&info);
    descinit_(descy,&n,&one,&nb,&one,&zero,&zero,&ctxt,&bl,&info);

    /** Initialize matrix A and right-hand-side */
    unsigned int  seed=(float) 100*(myrow+2)*200*(mycol+8); srand(seed);
    block_A(bl,A);
    if (myrow==0) block_x (bl, x);
    Cblacs_barrier(ctxt,"A");

    pdgemv_("No Transpose",&n,&n,&alpha,A,&one,&one,descA,
	         x,&one,&one,descx,&one,&beta,y,&one,&one,descy,&one);
	
	/** Solve the linear system */   
 
    if (iam==0) clock0=clock();
		
    pdgesv_(........);
	
	if (iam==0)
	   if (info!=0) 
	   {  printf("\n error in pdgesv: info=%d\n",info);
	      exit(0);
	   }
	   else 
	   {  clock1=clock(); printf("\n cputime: %e\n",(double)(clock1 - clock0)/CLOCKS_PER_SEC);
	      gflops=1e-9*n*n*n/3/(((double)(clock1 - clock0)/CLOCKS_PER_SEC));
	   }

    /** Compute  ||x-y|| */
	alpha = ...;
    pdaxpy_(&n,&alpha,x,&one,&one,descx,&one,y,&one,&one,descy,&one);
	pdnrm2_(&n,&norm,y,&one,&one,descy,&one);
	
    /** Print the results */
	if (iam==0)
	{  printf("\n number of GFLOPS: %e\n",gflops);
	   printf("\n ||x-y||/||x||=%e\n\n",norm/=...); 
	}

    /* Release process grid */
    Cblacs_gridexit(ctxt);

    /* Shut down MPI */
    MPI_Finalize();
    
    return 0;
}