program example1
      implicit none
*
*     simple example to show how to generate a ScaLAPACK matrix
*     contribution from Ed d'Azevedo, ORNL, 2005
*
      integer BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
     $        LLD_, MB_, M_, NB_, N_, RSRC_
      parameter ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
     $            CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
     $            RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
*
      integer descA(DLEN_)
*
      integer lwork
      parameter(lwork=10*1000)
      double precision work(lwork)
*
      integer Asize
      parameter(Asize=32*1000*1000)
*
      integer nout
      parameter(nout=16)
*
      double precision aij
      double precision A(Asize)
*
      integer m,n,mb,nb
      integer info, ierr(1)
      integer iam,nprocs
      integer icontext, myprow,mypcol,nprow,npcol
      integer rsrc,csrc,lld,Locp,Locq,Aneed
*
      integer ia,ja,irprnt,icprnt
      logical isroot, isok
*
      integer prow,pcol,lrindx,lcindx,ipos
      logical ismypcol, ismyprow, do_print
*
      double precision t1,t2
      double precision MPI_Wtime
      external MPI_Wtime
*
      integer numroc, indxg2p
      external numroc, indxg2p
      external infog2l, descinit
*
      mb = 50
      nb = 50
      m = 400
      n = 400
*
*     -----------------------
*     setup blacs environment
*     -----------------------
*
      call blacs_pinfo(iam,nprocs)
*
      do nprow=int(sqrt(real(nprocs)))+1,1,-1
         npcol = nprocs/nprow
         if (nprow*npcol.eq.nprocs) goto 11
      enddo
 11   continue
*
      call blacs_get(-1,0,icontext)
      call blacs_gridinit(icontext,'Col-major',nprow,npcol)
*
      call blacs_gridinfo(icontext,nprow,npcol,myprow,mypcol)
      isroot = (myprow.eq.0).and.(mypcol.eq.0)
*
      if (isroot) then
         write(*,*) 'nprow,npcol ', nprow,npcol
         write(*,*) 'm,n ', m,n
         write(*,*) 'mb,nb ', mb,nb
      endif
*
*     ---------------------------------------------------------
*     compute local extent and allocate storage for local piece
*     ---------------------------------------------------------
*
      csrc = 0
      rsrc = 0
*
      Locq = numroc(n,nb,mypcol,csrc,npcol)
      Locq = max(1,Locq)
      Locp = numroc(m,mb,myprow,rsrc,nprow)
*
      lld = max(Locp,1)
      Aneed = lld*Locq
      isok = (Aneed.le.Asize)
      if (.not.isok) then
         if (isroot) then
            write(*,*) 'increase Asize to ',Aneed + 1
         endif
         goto 999
      endif
*
*     ----------------
*     setup descriptor
*     ----------------
*
      call descinit(descA,m,n,mb,nb,rsrc,csrc,icontext,lld,info)
      ierr(1) = info
      call igsum2d(icontext,'All',' ',1,1,ierr,1,-1,-1)
      isok = (info.eq.0)
      if (.not.isok) then
         if (isroot) then
            write(*,*) 'descinit returns info = ',info
         endif
         goto 999
      endif
*
*     ------------------------------
*     simple method to setup matrix
*
*     (ia,ja) are global indices
*     ------------------------------
*
      call blacs_barrier(icontext,'All')
      t1 = MPI_Wtime()
*
      do ja = 1,descA(N_)
         pcol = indxg2p(ja,descA(NB_),mypcol,descA(CSRC_),npcol)
         ismypcol = (pcol.eq.mypcol)
         if (ismypcol) then
            do ia=1,descA(M_)
               prow = indxg2p(ia,descA(MB_),myprow,descA(RSRC_),nprow)
               ismyprow = (prow.eq.myprow)
               if (ismyprow) then
                  call infog2l(ia,ja,descA,nprow,npcol,myprow,mypcol,
     &                 lrindx,lcindx, rsrc,csrc)

                  ipos = lrindx + (lcindx-1)*descA(LLD_)
* 
                  call generate_Aij(m,n,ia,ja,aij)
                  A(ipos) = aij
               endif
            enddo
         endif
      enddo
*           
      call blacs_barrier( icontext, 'All')
      t2 = MPI_Wtime()
      if (isroot) then
         write(*,*) 'time to build matrix is ',t2-t1,' sec'
      endif
*
*        --------------------------------------------
*        if matrix is not too big, print out content
*        for debugging
*        --------------------------------------------
*
      do_print = ((m*n.le.200*1000).and.(lwork.ge.mb))
      if (do_print) then
*
         ia = 1
         ja = 1
         irprnt = 0
         icprnt = 0
         call pdlaprnt(m,n,A,ia,ja,descA,irprnt,icprnt,'A',nout,work)
*
      endif
*
999   continue
*
*     ---------------
*     prepare to exit
*     ---------------
*
      call blacs_barrier(icontext, 'All')
      call blacs_gridexit( icontext )
      call blacs_exit(0)
      stop
      end
*
      subroutine generate_Aij( m,n, ia,ja, aij )
      implicit none
      integer m,n, ia,ja
      double precision aij
*
      aij = dble(ia) + dble(ja-1)*dble(m)
      return
      end
*