/***********************************************************
 * dot_pblas.c -- parallel dot product                     *
 *                                                         *
 * Compute the dot product of two N-vectors using the      *
 * Parallel Basic Linear Algebra Subroutines (PBLAS)       *
 * 1.11.2002 PA                                            *
 ***********************************************************/

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

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

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

main(int argc, char* argv[]) {
    int myid;                              /* process rank */
    int p;                         /* number of processes  */
    int n, nb, i, i0;
    int ctxt, prow, pcol, myrow, mycol;    /* BLACS stuff */
    int descx[9], descy[9], info;
    int size = N * sizeof(double);
    double s, *x, *y;

    /* 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;
    x = (double*) malloc(size);
    y = (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;

    descinit_(descx,&N,&ONE,&nb,&ONE,&ZERO,&ZERO,&ctxt,&N,&info);
    descinit_(descy,&N,&ONE,&nb,&ONE,&ZERO,&ZERO,&ctxt,&N,&info);

    /* Initialize vectors */
    for (i = 0; i < n; i++){
      x[i] = 1.0;
      y[i] = (double)i0 + (double)i;
    }

    /* Check result */
    s = 0.0;
    pddot_(&N,&s,x,&ONE,&ONE,descx,&ONE,y,&ONE,&ONE,descy,&ONE);
    if (myid == 0) printf("s = %20.2f  (should be %d)\n", s, N*(N-1)/2);

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

    /* Shut down MPI */
    MPI_Finalize();

} /* main */