/*

   Benchmark #8 -- Dynamic Program

   Large sparse matrices (Part b)

   Assumption: All values in the A matrices are nonnegative
               (Since these are state transition probabilities, this
                is a reasonable assumption.)

  The Problem:

  Given: A[1],A[2],...A[K],  N x N floating point matrices
          D, an 8-bit integer vector of length T where
              1 <= D(t) <= K for t = 1,2,...T

  Calculate: Y[0],Y[1],...Y[T] 64-bit floating point vectors of length N
            P[1],P[2],...P[T] 32-bit integer vectors of length N
            B, a T+1 long 32-bit integer vector where
               1 <= B(t) <= N for t = 0,1,2,...T
  according to the following algorithm:

     Y[0](i) = 1 for i = 1,2,...N
     For t = 1,2,...,T
         k = D(t)
         For j = 1,2,...,N
             Y[t](j) = max Y[t-1](i)*A[k](i,j)
                        i
             P[t](j) = any i which maximized Y[t](j)
         Next j
     Next t

     Let Z = max Y[T](i), and let I be any i which maximized Z.
              i

     Set B(T) = I
     For t = T-1,T-2,...,1,0
         Set B(t) = P[t+1](B(t+1))
     Next t

  Output Z and B(t) for t = 0,1,2,...T.

  Parameters:    N = 10**6   K = 2   T = 2000

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

   Main Program for the large, sparse matrix version of Benchmark #8

   Call time parameters:

     N    =  Size of the A matrices
                Maximum = Default = 10**6

     K    =  Number of A matrices
                Maximum = 5
                Default = 2

     L    =  The number of non-zero entries per row and column.
                Maximum = 5
                Default = 4

     T    =  Length of D array
                Maximum = Default = 2000

     D2   =  Partitions of T.  D1 and D2 are greater than 1, and
             their product is T.  (D1 is computed from D2)
                Maximum = Default = 50

************************************************************************* */
#include "bench8sc.h"
#include <unistd.h>
#include <time.h>
#include <sys/time.h>
#include <stdio.h>
#include <mpi.h>

static char cvs_info[] = "BMARKGRP $Date: 2005/01/10 21:04:37 $ $Revision: 1.2 $ $RCSfile: m8sc.c,v $ $Name: rel_5 $"; 


/* Define these matrices global to keep the stack from getting too large. */

  int32 ia[NONZERO][MATRIXSIZE];
  long d[DLEN];			/* Which matrix to use in each step */
  int32 b[DLEN+1];		/* Best path */

int main( int argc, char *argv[] ) 
{

  int ier;			/* Error flag */
  int n = -1;			/* The size of the A matrices */
  int k = -1;			/* The number of A matrices */
  int l = -1;			/* The number of non-zero entries */
  int d1, d2=-1, t=-1;                /* d1 * d2 = t */
  double z,zz;
  int i,j,e_code = 0;
  double cputime(), wall();
  double setcputime;               /* CPU time for setup */
  double nsetcputime;               /* CPU time for setup */
  double setrealtime;              /* Real time for setup */
  double runcputime;               /* CPU time for benchmark */
  double runrealtime;              /* Real time for benchmark */


  /* The array to hold the matrices, and the sparse index */
  double *a;
  extern int32 ia[NONZERO][MATRIXSIZE];
  extern long d[DLEN];			/* Which matrix to use in each step */
  extern int32 b[DLEN+1];			/* Best path */
  int iii,ii;

  int mype = 0;			/* rank of process	*/
  int npes;			/* number of processes	*/
  int ionode = 0;		/* rank of pe doing our IO*/
  int s8s();
  void c8s(), p8s(), r8s();


  /* Get the input parameters */
  /*processargs(argc,argv);  */	/* Need to write this */
  /* Get N, the matrix size */
  if (n <= 0) n = MATRIXSIZE;
  if (n > MATRIXSIZE)
    {
      printf("Maximum allowed value for n is %d\n",MATRIXSIZE);
      exit(-1);
    }
  /* Get K, the number of matrices */
  if (k <= 0) k = NMATRICES;
  if (k > NMATRICES)
    {
      printf("Maximum allowed value for k is %d\n",NMATRICES);
      exit(-1);
    }
  /* Get l, the number of non-zero entries per row and column */
  if (l <= 0) l = NONZERO;
  if (l > NONZERO)
    {
      printf("Maximum allowed value for l is %d\n",NONZERO);
      exit(-1);
    }
  /* Get t, the length of the D array */
  if (t <= 0) t = DLEN;
  if (t > DLEN)
    {
      printf("Maximum allowed value for t is %d\n",DLEN);
      exit(-1);
    }
  /* Get d2,  the partitions of T */
  if (d2 <= 0) d2 = PART;
  if (d2 > PART)
    {
      printf("Maximum allowed value for d2 is %d\n",PART);
      exit(-1);
    }
    
    d1 = t/d2;
  if (d1 > PART)
    {
      printf("Maximum allowed value for d1 is %d\n",PART);
      exit(-1);
    }
  d1 = d1 * d2;
  if (d1 != t) 
    {
      printf ("d2 does not divide t; t = %d, d2 = %d\n",t,d2);
      exit(-1);
    }

  if(MPI_Init( &argc, &argv) != MPI_SUCCESS){
     printf("Initialization error\n");
     exit(1);
   }

  if((e_code = MPI_Comm_rank( MPI_COMM_WORLD, &mype)) != MPI_SUCCESS)
     err(e_code,"Rank");
  if((e_code = MPI_Comm_size(MPI_COMM_WORLD, &npes)) != MPI_SUCCESS)
     err(e_code,"Size");


  a = (double *)malloc(sizeof(double) * NMATRICES * NONZERO * MATRIXSIZE/npes);
  if (a == (double *)0) {
    printf("Can't malloc space for a matrix\n");
    MPI_Finalize(); 
    exit(-1);
    }


  /* Begin timing here */
  setcputime = cputime();
  setcputime = -cputime();
  
  setrealtime = -wall();

  /* Initialize timing variables */

  if (s8s(n,k,l,t,d,a,ia,mype,npes) == -1){
	err(e_code, "s8s");
     }

  /* End timed setup */
  

  setcputime += cputime(); 
  setrealtime += wall();

  

/*
  printf("Time to set up\n");
  printf("CPU      = %12.4f\n",setcputime);
  printf("WallCloc = %12.4f\n",setrealtime);
*/



  /* Begin timing here */
  runcputime = -cputime();
  runrealtime = -wall();

  /* Do the real work */
  p8s(a,ia,d,b,n,k,l,t,d2,&z, mype, npes); 


  /* End timed run */
  runcputime += cputime();
  runrealtime += wall();

    /* Check the results */
    c8s(n,k,l,t,a,ia,b,d,z,&zz,&ier,mype,npes); 
  if (mype == 0){
    /* Output results */
    r8s(n,k,l,t,d2,z,zz,b,setcputime,setrealtime,runcputime,runrealtime,ier);
    }
  MPI_Finalize(); 
}