      PROGRAM TEST
      implicit real*8(a-h,o-z)
      include "mpif.h"
c
c     /***************************************************************/
c     /* solve a linear system A*x=b by using the SCALAPACK          */
c     /* routine PDGESV.                                             */
c     /* uses 2d block-cyclic distribution of data. the processor    */
c     /* grid is (nxprocs,nyprocs) and the data size is (m,m). the   */
c     /* "block size" (mbb,nbb) is a free parameter.                 */
c     /*                                                             */
c     /* compiling :                                                 */
c     /* make -f makefile_linux clean                                */
c     /* make -f makefile_linux MAIN2                                */
c     /*                                                             */
c     /* running (on 'x' processors)                                 */
c     /* mpirun -np x ./MAIN2                                        */
c     /*                                                             */
c     /*-------------------------------------------------------------*/
c     /* based upon the PDSCAEX example of SCALAPACK                 */
c     /* markus uhlmann                                              */
c     /***************************************************************/
c
c     /* (global) problem size and number of processors */
      parameter(m=2000,n=m,nxprocs=2,nyprocs=2)
      parameter(nrhs=1)
c
c     /* desired block size of block-cyclic description */
      parameter(mbb=5,nbb=5)
c
c     /* derived dimensions: it is decomposed along both indices */
      parameter(nproc=nxprocs*nyprocs)
      parameter(mpdim=m/nxprocs+mbb,npdim=n/nyprocs+nbb)
c
c     /* (local) arrays */
      dimension aloc(mpdim*npdim),adebug(mpdim*npdim)
      dimension bloc(mpdim*nrhs)
      integer ipiv(mpdim+mbb+1) !swap-info from pivoting
      dimension work1(mpdim+mbb+1) !temp-buffer for output
c
c     /* internal scalapack dimensions */
      parameter(IDESCLEN_=9)
c
c     /* array descriptors */
      integer idesca(IDESCLEN_),idescb(IDESCLEN_)
c
c     /* mpi related stuff */
      integer istat(MPI_STATUS_SIZE)
c
c     /* other stuff */
      character*80 filename
c-----|---------------------------------------------------------------|
c     /* initialize MPI */
      call MPI_INIT(ierr)
      t0=MPI_WTIME()
      call MPI_COMM_RANK(MPI_COMM_WORLD,myid,ierr)
      call MPI_COMM_SIZE(MPI_COMM_WORLD,numprocs,ierr)
c
      if(numprocs.ne.nproc)then !thats not what we compiled for
         write(*,*)'difference in no of procs: ',numprocs,nproc
         goto 111
      endif
c
c     /* initialise BLACS process contexts: column-major process grid */
      call SL_INIT(ictxt,nxprocs,nyprocs)
      call BLACS_GRIDINFO(ictxt,nxprocs_0,nyprocs_0,my_row,my_col)
      if(myid.eq.0)then
         write(*,*)'Linear System Solution; Size: ',m
         write(*,*)'Processor Grid Info:',nxprocs,nyprocs
      endif
      write(*,*)'my_id:',myid,' my_col:',my_col,' my_row:',my_row
c
c     /* create the SCALAPACK array descriptors for 2d-block-cyclic */
      mblock=max(1,min(maxblock_scalapack(m,nxprocs),mbb))
      nblock=max(1,min(maxblock_scalapack(n,nyprocs),nbb))
      write(*,*)my_row,my_col,' initializing descriptors: ',
     $     mblock,nblock
      call DESCINIT(idesca,m,n   ,mblock,nblock,0,0,ictxt,mpdim,info)
      call DESCINIT(idescb,m,nrhs,mblock,nblock,0,0,ictxt,mpdim,info)
c
c     /* check local array dimensions */
      mp=NUMROC(m,mblock,my_row,0,nxprocs)
      np=NUMROC(n,nblock,my_col,0,nyprocs)
      if(mp.gt.mpdim.or.np.gt.npdim)then
         write(*,*)my_row,my_col,' insufficient local array dim: ',
     $        mp,mpdim,np,npdim
         call BLACS_ABORT(ictxt,0)
         goto 111
      endif
c
c     /* assign some data locally */
      call assign_matrix(aloc,mpdim,idesca)
      call assign_vector(bloc,idescb)
c
#ifdef DEBUG
c     /* write matrix & vector to screen */
      write(*,*)'the matrix & vector: A,b'
      iunit=6
      call PDLAPRNT(m,n,aloc,1,1,idesca,0,0,'A',iunit,work1)
      call PDLAPRNT(m,nrhs,bloc,1,1,idescb,0,0,'B',iunit,work1)
#endif /* DEBUG */
c
c     /* synchronize for timing */
      if(myid.eq.0)write(*,*)'start timing - next solving...'
      call MPI_BARRIER(MPI_COMM_WORLD,ierr)
      t1=MPI_WTIME()
c
c     /* call the SCALAPACK solver: */
      call PDGESV(m,nrhs,aloc,1,1,idesca,ipiv,bloc,1,1,idescb,info)
c     
      t2=MPI_WTIME()
      if(info.ne.0)then
         write(*,*)my_row,my_col,' PDGESV: info=',info
      endif
c     
c     /* write the result to unit 'iunit' */
      iunit=67
      if(my_row.eq.0.and.my_col.eq.0)
     $     write(*,*)'SCALAPACK - the result is in unit ',iunit
      call PDLAPRNT(m,nrhs,bloc,1,1,idescb,0,0,'B',iunit,work1)
c     
      if(numprocs.eq.1)then
c     /* solution using LAPACK for single processor - for comparison */
         write(*,*)'solving with LAPACK:'
         call assign_matrix1(aloc,mpdim,m,n)
         call assign_vector1(bloc,m)
         t3=MPI_WTIME()
         call DGESV(m,nrhs,aloc,mpdim,ipiv,bloc,mpdim,info)
         t4=MPI_WTIME()
         if(info.ne.0)then
            write(*,*)my_row,my_col,' DGESV: info=',info
         endif
         iunit=66
         write(*,*)'result from LAPACK is in unit ',iunit
         do i=1,m
            write(iunit,'(1x,i5,1x,e30.18)')i,bloc(i)
         enddo
      endif
      if(myid.eq.0)then
         write(*,*)'SCALAPACK solution took: ',t2-t1,' sec ',
     $        '(unit 55)'
         write(55,200)m,n,mblock,nblock,
     $        nxprocs,nyprocs,nxprocs*nyprocs,t2-t1
 200     format(1x,7(i5,1x),10(e15.8))
         if(numprocs.eq.1)then
            write(*,*)'LAPACK solution took: ',t4-t3,' sec '
         endif
      endif
c
c     /* finalize BLACS */
      call BLACS_GRIDEXIT(ictxt)
      call BLACS_EXIT(1)        !1 -> you can still use MPI afterwards
c
c     /* finalize MPI */
 111  write(*,*)myid,' about to finalize MPI'
      call MPI_FINALIZE(ierr)
      write(*,*)'bye. ',myid      
c     
      stop
      end 
c-----|---------------------------------------------------------------|
      integer function maxblock_scalapack(n,numprocs)
      implicit none
c     /* determine the maximum possible block size for scalapack-    */
c     /* style block-data distribution.                              */
      integer n,numprocs
c
      maxblock_scalapack=n/numprocs
      if(mod(n,numprocs).ne.0)maxblock_scalapack=maxblock_scalapack+1
      return
      end function maxblock_scalapack
c-----|---------------------------------------------------------------|
      subroutine assign_matrix(aloc,llda,idesca)
      implicit real*8(a-h,o-z)
c
c     /* assigns a distributed matrix.                               */
c     /* 2d-block-cyclic, scalapack-style.                           */
c     /* NOTE: NOT efficient                                         */
c
      dimension aloc(llda,*),idesca(*)
c
      parameter(ICTXT_=2,M_=3,N_=4)
c----+----------------------------------------------------------------|
c     /* get grid parameters */
      ictxt=idesca(ICTXT_)
      m=idesca(M_)
      n=idesca(N_)
      call BLACS_GRIDINFO(ictxt,nxprocs,nyprocs,my_row,my_col)
c
c     /* assign local part of the data */
      do jglob=1,n
         do iglob=1,m
            call local_from_global2d(idesca,iglob,jglob,
     $           nxprocs,nyprocs,
     $           mrow,mcol,iloc,jloc)
            if(my_row.eq.mrow.and.my_col.eq.mcol)then
               aloc(iloc,jloc)=fmat(iglob,jglob)
            endif
         enddo
      enddo

c
      return
      end
c----+----------------------------------------------------------------|
      subroutine assign_matrix1(aloc,llda,m,n)
      implicit real*8(a-h,o-z)
c     /* assigns a matrix */
c
      dimension aloc(llda,*)
c----+----------------------------------------------------------------|
      do j=1,n
         do i=1,m
            aloc(i,j)=fmat(i,j)
         enddo
      enddo
c
      return
      end
c----+----------------------------------------------------------------|
      subroutine assign_vector(bloc,idescb)
      implicit real*8(a-h,o-z)
c
c     /* assigns a distributed vector                                */
c     /* 1d-block-cyclic, i.e. the data is distributed to the first  */
c     /* column (my_col=0) of a 2d cartesian grid, block-cyclic over */
c     /* the 0:nxprocs-1 processors.                                 */
c
      dimension bloc(*),idescb(*)
c
      parameter(ICTXT_=2,MB_=5,M_=3,N_=4)
c----+----------------------------------------------------------------|
c     /* get grid parameters */
      ictxt=idescb(ICTXT_)
      mb=idescb(MB_)
      m=idescb(M_)
      n=idescb(N_)
      call BLACS_GRIDINFO(ictxt,nxprocs,nyprocs,my_row,my_col)
c
c     /* assign local part of the data */
      do iglob=1,m
         call local_from_global(iglob,mb,nxprocs,
     $        mrow,iloc)
         if(my_row.eq.mrow.and.my_col.eq.0)then
            bloc(iloc)=fvec(iglob)
         endif
      enddo
c
      return
      end
c----+----------------------------------------------------------------|
      subroutine assign_vector1(bloc,m)
      implicit real*8(a-h,o-z)
c     /* assigns a vector*/
c
      dimension bloc(*)
c----+----------------------------------------------------------------|
      do i=1,m
         bloc(i)=fvec(i)
      enddo
c
      return
      end
c----+----------------------------------------------------------------|
      real*8 function fmat(i,j)
      implicit real*8(a-h,o-z)
c
      if(i.eq.j)then
         fmat=1.d4
      else
         fmat=dfloat(i+j)
      endif
c
      return
      end function fmat
c----+----------------------------------------------------------------|
      real*8 function fvec(i)
      implicit real*8(a-h,o-z)
c
      fvec=dfloat(207-i)
c
      return
      end function fvec
c----+----------------------------------------------------------------|
      subroutine local_from_global2d(idesc,iglob,jglob,
     $     nxprocs,nyprocs,
     $     mrow,mcol,iloc,jloc)
      implicit real*8(a-h,o-z)
c
c     /* computes processor coordinates and the local indices of a   */
c     /* 2d array which is block-cyclic distributed                  */
c     /* NOTE: supposes we deal out data starting with the 0th pro-  */
c     /*       cessor.                                               */
c     /*                                                             */
c     /* input                                                       */
c     /* idesc          array descriptor of 2d-block-cyclic array    */
c     /* iglob,jglob    global coordinates of an element             */
c     /* nxprocs,nyprocs total no. of processor in row/col           */
c     /*                                                             */
c     /* output                                                      */
c     /* mrow,mcol      coordinates of processor which holds it      */
c     /* iloc,jloc      local indices of data                        */
c
      dimension idesc(*)
c
      parameter(ICTXT_=2,M_=3,N_=4,MB_=5,NB_=6,RSRC_=7,CSRC_=8,LLD_=9)
c----+----------------------------------------------------------------|
      m=idesc(M_)
      n=idesc(N_)
      mb=idesc(MB_)
      nb=idesc(NB_)
c
      call local_from_global(iglob,mb,nxprocs,
     $     mrow,iloc)
      call local_from_global(jglob,nb,nyprocs,
     $     mcol,jloc)
c
      return
      end subroutine local_from_global2d
c----+----------------------------------------------------------------|
      subroutine local_from_global(iglob,mb,nprocs,
     $     myid,iloc)
      implicit real*8(a-h,o-z)
c
c     /* computes processor coordinates and the local indices of an  */
c     /* array which is block-cyclic distributed                     */
c     /* NOTE: supposes we deal out data starting with the 0th pro-  */
c     /*       cessor.                                               */
c     /*                                                             */
c     /* input                                                       */
c     /* iglob          global coordinates of an element             */
c     /* mb             block size                                   */
c     /* nprocs          total no. of processor in row/col           */
c     /*                                                             */
c     /* output                                                      */
c     /* myid            coordinate  of processor which holds it     */
c     /* iloc           local indicx  of data                        */
c
c----+----------------------------------------------------------------|
      myid=mod((iglob-1)/mb,nprocs)
      l=(iglob-1)/(mb*nprocs)   !note: this intermediate step is important
                                !      because of integer arithmetic
                                !i.e. l=floor(...)
                                !then iloc=l*mb+...
      iloc=l*mb+mod(iglob-1,mb)+1
c
      return
      end subroutine local_from_global
c----+----------------------------------------------------------------|

