/***********************************************************
 *  *
 * Parallel Basic Linear Algebra Subroutines (PBLAS)
 
 mpicc.openmpi TS.c -lblacs-openmpi -lscalapack-openmpi -o TS  
 mpirun.openmpi -np 4  ./TS   *
 *                                         *
 ***********************************************************/

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

/***********************************************************
 *                      Globals                            *
 ***********************************************************/

int N=10, ZERO=0, ONE=1;;

main(int argc, char* argv[]) {
    int myid;                              /* process rank */
    int p;                         /* number of processes  */
    int n, nb, i,j, i0;
    int ctxt, prow, pcol, myrow, mycol;    /* BLACS stuff */
    int descx[9], desca[9], info;
    int size = N * sizeof(double);
    double s, *x, **a;
    
printf("avant size=%d\n", size);
    /* Start the MPI engine */
    MPI_Init(&argc, &argv);

    /* Find out number of processes */
    MPI_Comm_size(MPI_COMM_WORLD, &p);

    size = (size + p - 1)/p;
    printf("size=%d\n", size);
    x = (double*) malloc(size);
    a = (double**) malloc(size);
   

    /* Find out process rank  */
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);

    /* Get a BLACS context */
    Cblacs_get(0, 0, &ctxt);

    /* Initialize a 1-dimensional grid */
    Cblacs_gridinit(&ctxt, "R", p, 1);
    Cblacs_gridinfo(ctxt, &prow, &pcol, &myrow, &mycol);

    nb = (N + p - 1)/p;     /* nb = ceil(N/p)  */
    i0 = myid*nb;
    n = nb;
    if (myid == p-1) n = N - i0;
     for(i=0;i<n;i++)
    a[i] = (double*) malloc((i) * sizeof(double));
    printf("P=%d nb=%d i0=%d n=%dprow=%d pcol=%d myrow=%d mycol=%d\n", myid,nb,i0,n,prow,pcol,myrow,mycol);

    descinit_(descx,&ONE,&N,&nb,&nb,&ZERO,&ZERO,&ctxt,&ONE,&info);
    
    descinit_(desca,&N,&N,&nb,&nb,&ZERO,&ZERO,&ctxt,&n,&info);
    for(i=0;i<=9;i++)
    printf("qqqqqq P%d descx[%d]=%d desca[i]=%d\n",myid,i,descx[i],desca[i]);

    /* Initialize vectors */
    for (i = 0; i < n; i++){
      //x[i] = 1.0;
       x[i] = (double)i0 + (double)(i+1);
      for(j=0;j<x[i]-1;j++)
      {
      a[i][j] = 1.0;
      printf("p=%d a[%d][%d]=%f\n",myid,i,j,a[i][j]);
      }
   
      printf("p=%d i=%d x[i]=%f\n",myid,i,x[i]);
    }

    /* Check result */
   
    
    //pstrsv_(UPLO, TRANS, DIAG, &N,A, &IA, &JA, DESCA,X, &IX, &JX, DESCX, &INCX);
    pstrsv_("Lower", "No Transpose", 'U', &N,a, &ONE, &ONE, desca,x, &ONE, &ONE, descx, &info);

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

    /* Shut down MPI */
    MPI_Finalize();

} /* main */