
c********************************************************************
c
c                        Example 15 vv_mult_blk_cyclic.f
c
c  Objective           : Vector_Vector Multiplication (Using Cyclic data 
c                        Partitioning ) 
c
c  Input               : Read files vdata1.inp for Vector_A and vdata2.inp 
c                        for Vector_B 
c
c  Output              : Process with Rank 0 prints the Vector-Vector 
c                        Multiplication results 
c
c  Necessary Condition : Number of Processes should be
c                        less than or equal to 8
c
c*******************************************************************



	 program main

	 include 'mpif.h'
	 parameter (Epsilon = 1.0E-10, DATA_SIZE = 16)

  	 integer	Numprocs, MyRank
  	 integer 	VectorSize, ScatterSize 
  	 integer      VectorSize_A, VectorSize_B 
  	 integer	index, iproc
  	 integer	Root
  	 real         Mybuffer_A(DATA_SIZE), Mybuffer_B(DATA_SIZE)  
	 real	 	MyFinalVector, FinalAnswer
  	 real 		Vector_A(DATA_SIZE), Vector_B(DATA_SIZE)
  	 real 		Buffer_A(DATA_SIZE), Buffer_B(DATA_SIZE)


C   ........MPI Initialisation .......

   	 call MPI_INIT(ierror) 
   	 call MPI_COMM_RANK(MPI_COMM_WORLD, MyRank, ierror)
   	 call MPI_COMM_SIZE(MPI_COMM_WORLD, Numprocs, ierror)

	 Root = 0


   	 if(MyRank .eq. Root) then 
 	    open(unit=12, file = './data/vdata1.inp')
	    read(12,*) VectorSize_A
	    read(12,*) (Vector_A(i), i=1,VectorSize_A)

 	    open(unit=12, file = './data/vdata2.inp')
	    read(12,*) VectorSize_B
	    read(12,*) (Vector_B(i), i=1,VectorSize_B)
	 endif

C     .......Allocate memory and read data for vector A .....

      call MPI_BCAST(VectorSize_A, 1, MPI_INTEGER, Root, 
     $		 MPI_COMM_WORLD, ierror)

      call MPI_BCAST(VectorSize_B, 1, MPI_INTEGER, Root, 
     $		 MPI_COMM_WORLD, ierror)

       if(VectorSize_A .ne. VectorSize_B) then
	  if(MyRank .eq. Root) then
          print *,"Two Vectors can not be multiplied - exit"
	  endif
	  goto 100
       else
        VectorSize = VectorSize_A
	endif
	 
   	 if(VectorSize .lt. Numprocs) then
	    if(MyRank .eq. Root) then
              print*,"VectorSize should be more than No of Processors"
	    endif
	    goto 100
	 endif

   	 if(mod(VectorSize, Numprocs) .ne. 0) then
     	   if(MyRank .eq. Root) then
	      	print*,"VectorSize Can not be Striped Evenly ..... "
	    endif
	    goto 100
	 endif

C    ....... Rearrangement of data for Cyclic distribution ....... 
	if(MyRank .eq. Root) then
	 index = 1
      	 Block_Size = VectorSize / Numprocs
 	 do iproc = 1, Numprocs
	      do global_index  = iproc, VectorSize, Numprocs
               Buffer_A(index) = Vector_A(global_index)
               Buffer_B(index) = Vector_B(global_index)
		         index = index + 1
	       enddo					
	 enddo					
	endif

C    ...... Scatter vector A and B ......

       ScatterSize = VectorSize / Numprocs
       call MPI_SCATTER(Buffer_A,ScatterSize,MPI_REAL,Mybuffer_A, 
     $		 ScatterSize,MPI_REAL,Root,MPI_COMM_WORLD, ierror)
			 
       call MPI_SCATTER( Buffer_B,ScatterSize,MPI_REAL,Mybuffer_B, 
     $		 ScatterSize,MPI_REAL,Root,MPI_COMM_WORLD,ierror)

C	  ........ Calculate partial sum ........ 

       MyFinalVector = 0.0
       do index = 1, ScatterSize
        MyFinalVector = MyFinalVector + 
     $ 	                Mybuffer_A(index)*Mybuffer_B(index)
       enddo

C .......... Collective computation : Final answer on process 0 ....
      call MPI_REDUCE(MyFinalVector, FinalAnswer, 1, MPI_REAL,
     $		 MPI_SUM, Root, MPI_COMM_WORLD, ierror)

	if(MyRank .eq. 0) then
 	  print*,'FinalAnswer', FinalAnswer
	 endif

 100	 call MPI_FINALIZE(ierror)
         stop
	 end



