src/hamiltonian_expectation_value/sparse_matrix_parallel.c

   1
   2                    //  PHYS 405                                                //
   3
   4                    //  expectation expectation_value of a sparse Hamiltonian matrix    //
   5                    //  H matrix filled with random numbers at random locations         //
   6                    //  parallel version                                                //
   7
   8                    //  assignment 2 - part b               //
   9
  10                                          // Michel Vallieres, Fall 2007 //
  11
  12 #include <stdio.h>
  13 #include <math.h>
  14 #include <stdlib.h>
  15 #include <mpi.h>          /* MPI header file */
  16
  17
  18
  19 int main( int argc, char * argv[] )
  20 {
  21   double Hij, symm, local_value, expectation_value, norm, *v;
  22   int    N, Nelements, i, j, element, Nlocal, Nshift;
  23   FILE   *fp;
  24   int    rank, size;
  25
  26                                           /* init into MPI */
  27   MPI_Init(&argc, &argv);
  28                                           /* my rank - my id */
  29   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  30                                           /* how many processes in the virtual machine */
  31   MPI_Comm_size(MPI_COMM_WORLD, &size);
  32
  33                                           /* data file size */
  34   fp = fopen( "/tmp/michel_sparse_Hamiltonian_size", "r" );
  35   fscanf( fp, "%d %d\n", &N, &Nelements);
  36   fclose( fp );
  37
  38                                           /* my data slice */
  39   Nlocal = Nelements / size;
  40   Nshift = Nlocal * rank;
  41   if ( rank == size-1 )
  42     Nlocal = Nelements - Nshift;
  43   fprintf( stderr, " Rank, size, H size, non-zero, Nlocal, Nshift: %d %d %d %d %d %d\n",
  44            rank, size, N, Nelements, Nlocal, Nshift );
  45
  46                                          /* unit vector v*v = 1 */
  47   v = (double *)malloc( N * sizeof(double) );
  48   norm = 1.0 / sqrt( (float)N );
  49   for ( i=0 ; i<N ; i++ )
  50       v[i] = norm;
  51
  52                                           /* data file */
  53   fp = fopen( "/tmp/michel_sparse_Hamiltonian", "r" );
  54
  55                                            /* skip values not pertaining to rank */
  56   if ( Nshift != 0 )
  57     {
  58        for ( element=0 ; element<Nshift ; element++ )
  59                fscanf( fp, "%d %d %lf", &i, &j, &Hij );
  60     }
  61
  62                                           /* loop over non-zero Hij */
  63   local_value = 0.0;
  64   for ( element=0 ; element<Nlocal ; element++ )
  65     {
  66                                           /* read  in Hij */
  67       fscanf( fp, "%d %d %lf", &i, &j, &Hij );
  68                                           /* symmetric matrix -- 2 off-diag */
  69       symm = 1.0;
  70       if ( i != j )
  71         symm = 2.0;
  72                                           /* matrix multiplication */
  73       local_value = local_value + symm * v[i] * Hij * v[j];
  74
  75     }
  76
  77   fclose( fp );
  78   free( v );
  79
  80
  81   fprintf( stderr, " Local contribution in rank: %f %d\n",
  82                                                 local_value, rank );
  83
  84                                  // sum up the contributions
  85   MPI_Reduce ( &local_value, &expectation_value, 1,
  86                MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
  87
  88                                  // final answer
  89   if ( rank == 0 )
  90     fprintf( stderr, "\n Expectation value: %f\n\n", expectation_value );
  91
  92   MPI_Finalize();
  93   exit(0);
  94
  95 }
  96