!  sctest1.f  02 February 1998
!
!  A simple ScaLAPACK example.
!
      program EXAMPLE1
!
!***********************************************************************
!
      real A ( 5, 5 )
      real B ( 5 )
      integer DESCA ( 9 )
      integer DESCB ( 9 )
      integer i
      integer IAM
      integer ICTXT
      integer ihi
      integer ilo
      integer INFO
      integer IPIV ( 10 )    ! Dimensioning of IPIV is a question.
      integer istore
      integer j
      integer jhi
      integer jlo
      integer jstore
      integer M
      integer MYCOL
      integer MYROW
      integer N
      integer NBRHS
      integer NPCOL
      integer NPROCS
      integer NPROW
!
!***********************************************************************
!  1. SET UP THE PROCESSES, AND DEFINE THE GRID.
!***********************************************************************
!
!  Define the number of rows and columns in the processor grid.
!
      NPROW = 2
      NPCOL = 2
!
!  BLACS_PINFO tells this process its number, and the total number of
!  processes.
!  
      call BLACS_PINFO ( IAM, NPROCS )
!
!  If no processes have been set up, then we must set up all NPROW * NPCOL
!  processes via a call to BLACS_SETUP.  
!
!  We only want one process to do the setup.  There is always a process
!  number 0, so that's the one we use to take care of this chore.
!
      if ( NPROCS .lt. 1 ) then

        if ( IAM .eq. 0 ) then
          NPROCS = NPROW * NPCOL
        end if

        call BLACS_SETUP ( IAM, NPROCS )

      end if
!  
!  Ask BLACS_GET to return the default system context in ICTXT.
!
      call BLACS_GET ( -1, 0, ICTXT )
!
!  BLACS_GRIDINIT sets up the process grid of NPROW by NPCOL
!  processes, in row-major order.  
!
!  Note that we pass the default system context in ICTXT, and
!  BLACS_GRIDINIT replaces it with the BLACS context for this grid.
!
      call BLACS_GRIDINIT ( ICTXT, 'Row-major', NPROW, NPCOL )
!
!  BLACS_GRIDINFO tells each process the "position" it has been
!  assigned in the process grid.
!
!  ICTXT is input, the context handle.
!  NPROW and NPCOL are output, the number of processor rows and columns.
!  MYROW and MYCOL are output, this process's row and column.
!
      call BLACS_GRIDINFO ( ICTXT, NPROW, NPCOL, MYROW, MYCOL )

      if ( MYROW .eq. -1 ) then
        call BLACS_EXIT ( 0 )
        stop
      end if
!
!***********************************************************************
!  2. DEFINE THE MATRIX AND RIGHT HAND SIDE.
!***********************************************************************
!  
!  Set the array descriptors for A and B.
!
      M = 10
      N = 10
      NOUT = 6
      call DESCINIT ( DESCA, M, N, 5, 5, 0, 0, ICTXT, 5, INFO )

      call DESCINIT ( DESCB, N, 1, 5, 1, 0, 0, ICTXT, 5, INFO )
!
!  Store the portion of A and B that belongs to this processor.
!
      do i = 1, 5
        do j = 1, 5
          A(I,J) = min ( MYROW*5 + i, MYCOL*5 + j )
        end do
      end do

      if ( MYCOL .eq. 0 ) then
        do i = 1, 5
          if ( MYROW .eq. 1 .and. i .eq. 1 ) then
            B(I) = 1.0
          else
            B(I) = 0.0
          end if
        end do
      end if
!
!***********************************************************************
!  3. SOLVE THE SYSTEM.
!***********************************************************************
!  
      if ( MYROW .eq. 0 .and. MYCOL .eq. 0 ) then
        write ( *, * ) ' '
        write ( *, * ) '  SCTEST1'
        write ( *, * ) '  ScaLAPACK demonstration.'
        write ( *, * ) ' '
        write ( *, * ) '  Solve A*X=B using PSGESV, for dense matrices.'
        write ( *, * ) '  The matrix A is defined as A(I,J) = min(I,J).'
      end if
*
      IF( IAM.EQ.0 ) THEN
         WRITE( NOUT, FMT = * )
         WRITE( NOUT, FMT = * ) 'Matrix A:'
         WRITE( NOUT, FMT = * )
      END IF
      CALL PSLAPRNT( N, N, A, 1, 1, DESCA, 0, 0,
     $               'A', NOUT, IPIV )
      IF( IAM.EQ.0 ) THEN
         WRITE( NOUT, FMT = * )
         WRITE( NOUT, FMT = * ) 'Matrix B:'
         WRITE( NOUT, FMT = * )
      END IF
      CALL PSLAPRNT( N, 1, B, 1, 1, DESCB, 0, 0,
     $               'B', NOUT, IPIV )
*
!
!  Call the ScaLAPACK routine PSGESV to solve the linear system A*X=B.
!
      call PSGESV ( N, 1, A, 1, 1, DESCA, IPIV, B, 1, 1, DESCB, INFO )
*
      IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
         WRITE( NOUT, FMT = * )
         WRITE( NOUT, FMT = * ) 'INFO code returned by PDGESV = ', INFO
         WRITE( NOUT, FMT = * )
         WRITE( NOUT, FMT = * ) 'Matrix X = A^{-1} * B'
         WRITE( NOUT, FMT = * )
      END IF
      CALL PSLAPRNT( N, 1, B, 1, 1, DESCB, 0, 0, 'X', NOUT,
     $               IPIV )
!
!***********************************************************************
!  4. SHUT DOWN THE PROCESSES.
!***********************************************************************
!  
!  Free the BLACS context.
!
      call BLACS_GRIDEXIT ( ICTXT )
!
!  Break the connection of this process to the BLACS.
!
      call BLACS_EXIT ( 0 )

      stop
      end