
     program gemm_example
!
!   demonstrates the usage of pdgemv, a parallel matrix-vector
!   multiply routine; requires BLACS, MPI, PBLAS libraries
!
      implicit none
!      include "mpif.h"
      include "param.h"

! parameters ...
      integer :: nprows=3, npcols=80
      integer :: m=74612736, n=80, nrhs=80
      integer :: mb, nb, nbrhs

      integer MAXLDA, MAXLDX, MAXSDX, MAXSDA
! local arrays ...
      real*8, dimension(:,:),allocatable ::A, x, b
      integer adesc(DESC_DIM), bdesc(DESC_DIM), xdesc(DESC_DIM)
! local scalars ...
      integer i, j, myrow, mycol, iam, nprocs, dim1, dim2
      integer icontext, aaa, bbb, param
! external function ...
      integer numroc
      external numroc

!! ______________________________________________
      mb=m/nprows; nb=n/npcols; nbrhs=nrhs/npcols

!Before performing any parallel processing, some BLACS subroutines must
!be called in order to initiate communication; define the number of
!process rows and columns (nprows and npcols); query for th
!e current rank ( iam), the number of processors (nprocs), as well as the
!current rank as represented in the 2D topology (myrow and mycol).

!
!  set up communication via BLACS
!
      call blacs_get(0,0,icontext)
      call blacs_gridinit(icontext,'c',nprows,npcols)
      call blacs_pinfo(iam,nprocs)

!  initialize using all of the target number of rows and columns

      call blacs_gridinfo(icontext,nprows,npcols,myrow,mycol)

!Next, a call to the subroutine desc_setup is made to register pertinent
!information regarding the array A, in an array called adesc. Similarly,
!calls to desc_setup are also made for the matrices X
! and B. genmat is then called to define A and X.


!  define the "desc" parameter for arrays A, x, and b

        MAXLDA= max(1,numroc(m, mb, myrow, 0, nprows))
        MAXSDA = max(1,numroc(n, nb, mycol, 0, npcols))
        MAXLDX = max(1,numroc(n, nb, myrow, 0, nprows))
        MAXSDX = max(1,numroc(nrhs, nbrhs, mycol, 0, npcols))
!      if(myrow==0)then
        allocate (A(MAXLDA,MAXSDA))
        allocate (X(MAXLDX,MAXSDX))
        allocate (B(MAXLDA,MAXSDX))
        aaa=1
        bbb=2
        
        do i=1,maxlda
        do j=1,maxsda
        a(i,j)=1.d0
        enddo
        enddo
        do i=1,maxldx
        do j=1,maxsdx
        x(i,j)=1.d0
        enddo
        enddo
!       endif
!       print *,a
      call desc_setup(adesc,m,   n,mb,   nb,icontext,MAXLDA)
      call desc_setup(xdesc,n,nrhs,nb,nbrhs,icontext,MAXLDX)
      call desc_setup(bdesc,m,nrhs,mb,nbrhs,icontext,MAXLDA)

! generate matrices A and x
!      call genmat(adesc,A,xdesc,x,myrow,mycol,MAXLDA,MAXLDX)

!With A and X defined, the PBLAS routine PDGEMM is called to compute the
!matrix-vector product, B. To wrap things up, calls to the BLACS
!routines are made.
!        write(*,*)'A', A
!        write(*,*)'X', X
!      ia=mb*myrow+1
!      ja=nb*mycol+1
!      ix=nb*myrow+1
!      jx=nbrhs*mycol+1
!      ib=mb*myrow+1
!      jb=nbrhs*mycol+1

       print *, 'before pdgemmmmmmmmmmmmmmmmmmmm'
      call pdgemm('N','N',m,nrhs,n,   &
           one,a,ia,ja,adesc,x,ix,jx,xdesc, &
           zero,b,ib,jb,bdesc)

       print *, 'after pdgemmmmmmmmmmmmmmmmmmmm'
! write solution b to screen
      if(mycol.lt.3) then
        dim1 = max(1,numroc(m, mb, myrow, 0, nprows))
        dim2 = max(1,numroc(nrhs, nbrhs, mycol, 0, npcols))
!        write(*,"('Process Grid Coor.:',4i3,':',(5f8.1))")myrow,mycol, &
!           dim1,dim2,((b(i,j),i=1,dim1),j=1,dim2)
      endif

      call blacs_gridexit(icontext)
      call blacs_exit()

      stop
      end

!Subroutine genmat



!Subroutine desc_setup
!In ScaLAPACK, there are supplemental routines provided as tools. Among
!these is a routine called DESCINIT to initialize descriptors that
!describe the matrix, such as matrix size (m x n, among othe
!r items. With IBM's PESSL package, a similar routine, called
!ARRAY_CREATE, exists. These two equivalent subroutines are sufficiently
!general to handle all situations. For this simple example, for t
!he sake of simplicity and portability between the two packages, a much
!simplified version is used here.


      subroutine desc_setup(Adesc,m,n,mb,nb,icontext,LDA)
      implicit none
      include "param.h"
      integer m, n, mb, nb, icontext, LDA, Adesc(DESC_DIM)

      Adesc(DTYPE_) = TWO_D_TYPE
      Adesc(M_) = m
      Adesc(N_) = n
      Adesc(MB_) = mb
      Adesc(NB_) = nb
      Adesc(RSRC_) = rsrc
      Adesc(CSRC_) = csrc
      Adesc(CTXT_) = icontext
      Adesc(LLD_) = LDA

      return
      end




