      PROGRAM SOLVE_LU

C     ===============================================================
C     ===============================================================
C     THIS PROGRAM SOLVE THE SYSTEM OF LINEAR SIMULTANEOUS EQNS AX = B
C     WHERE THE MATRIX A AND RHS VECTOR B ARE ASSIGNED RANDOM ENTRIES IN [0,1)
C     THE DIMENSION OF THE PROBLEM (I.E. NO. OF EQNS WE SOLVE) AND THE 
C     DIMENSIONS OF THE BLACS PROCESSOR GRID CAN BE SPECIFIED ON THE COMMAND LINE
C     AT RUN TIME

 

C      COMMAND LINE OPTIONS SUPPLIED AFTER mpirun -n ** para.exe 
C      SINGLE INTEGER -> N
C      TWO INTEGERS -> ILLEGAL, PROGRAM TERMINATES
C      THREE INTEGERS -> N, NPROW, NPCOL
C                        WHERE NPROW, NPCOL DENOTE NO. OF ROWS
C                        AND COLS IN BLACS PROCESS GRID    
C      FOUR INTEGERS -> N, NPROW, NPCOL, MB
C                       WHERE MB*MB BLACS BLOCK SIZE IS USED


C      *YOU* MUST ENSURE THE NUMBER OF PROCESSES CALLED BY MPIRUN 
C       = NPROW*NPCOL
 

C      RUN WITH e.g.:
C      mpirun -np 2 ./para.exe  3600  2 1 32  [SOLVE 3600 LIN EQNS USING 
C                               2X1 PROCESS GRID WITH 32X32 BLOCK SIZE]
 
C      ================================================================
C      ================================================================
      USE MPI
      IMPLICIT NONE
      INTEGER :: N
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: LOCAL_A
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: LOCAL_B 
      INTEGER :: ISTATUS
C     FOR LAPACK PDGESV CALL
      INTEGER  :: INFO,  NRHS, IA, JA, IB, JB
      INTEGER, ALLOCATABLE, DIMENSION (:) :: IPIV
c     FOR READING COMMAND LINE ARGUMENTS
      INTEGER :: IARGC, N_COMMAND_ARG
      CHARACTER :: ARGV*10 
C     WE USE FOLLOWING COMMAND LINE ARGUMENTS 
C     ARG 1 : N (DIMENSION OF PROBLEM)
C     ARG 2 : NPROW (NO. OF ROWS OF PROCESSES IN A RECTANGULAR ARRAY)
C     ARG 3 : NPCOL (NO. OF COLUMNS OF PROCESSES IN A RECTANGULAR ARRAY)
C     ARG 4 : BLACS BLOCK SIZE MB (BLOCKS ARE OF SIZE MB * MB) 
C   
c     FOR PARALLEL PROCESS ARRAY
      INTEGER  :: NPROW, NPCOL, ICTXT,MYROW, MYCOL, MB, NB, MLOC, NLOC
      INTEGER :: IDESCA(9), IDESCB(9)
      INTEGER :: IERR
      INTEGER :: NUMROC
      
 
c     for random number seed
      INTEGER :: ISEEDSIZE
      INTEGER, ALLOCATABLE, DIMENSION ( :) :: SEED
 
C      ----------------------------------------
C      -------  EXECUTABLE STATEMENTS   -------


C      ===============================================
C      READ IN COMMAND LINE ARGUMENTS IF PRESENT

      N_COMMAND_ARG = iargc()
      IF (N_COMMAND_ARG == 2) THEN
          WRITE(*,*) 'ILLEGAL NO. OF COMMAND LINE PARAMETERS'
          STOP
      ENDIF
      IF (N_COMMAND_ARG .GE. 1)THEN
          CALL GETARG(1,argv)
C          WRITE(*,*)'ARGV = ',ARGV
          READ (ARGV,'(I10)') N
      ELSE
          N = 100
      ENDIF


      IF (N_COMMAND_ARG .GE. 3)THEN
          CALL GETARG(2,argv)
          READ (ARGV,'(I10)') NPROW
          CALL GETARG(3,argv)
          READ (ARGV,'(I10)') NPCOL

      ELSE
          NPROW = 2
          NPCOL = 2
      ENDIF 

      IF (N_COMMAND_ARG .GE. 4)THEN
          CALL GETARG(4,argv)
          READ (ARGV,'(I10)') MB
      ELSE
          MB = 8
      ENDIF
      NB = MB
 
C     ==============================================
C     INITIALISE THE BLACS PROCESS GRID, FIND DIMENSIONS OF LOCAL
C     MATRICES / VECTORS AND ALLOCATE SPACE

      CALL SL_INIT(ICTXT, NPROW, NPCOL)
      CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )

      MLOC = NUMROC(N, MB, MYROW, 0, NPROW)
      NLOC = NUMROC(N, NB, MYCOL, 0, NPCOL)
 
      IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 )WRITE(*,*)
     @       'WE ARE SOLVING A SYSTEM OF ', N, ' LINEAR EQUATIONS'

      WRITE(*,*) 'PROC: ',MYROW, MYCOL,'HAS  MLOC, NLOC =', MLOC,NLOC

c      ==============================================
C     ALLOCATE SPACE FOR MATRIX A AND VECTORS B AND X

      WRITE(*,*) 'PROC: ',MYROW, MYCOL,' ALLOCATING SPACE ...'

      ALLOCATE ( LOCAL_A(MLOC,NLOC), STAT = ISTATUS )
      IF(ISTATUS .NE. 0) THEN
          WRITE(*,*)'UNABLE TO ALLOCATE LOCAL_A, PROCESS: ',MYROW,MYCOL
          STOP
      ENDIF

      ALLOCATE ( LOCAL_B(MLOC), STAT = ISTATUS )
      IF (ISTATUS /= 0) THEN
          WRITE(*,*)
     @ ' FAILED TO ALLOCATE SPACE FOR LOCAL_B, PROCESS: ',MYROW,MYCOL
          STOP
      ENDIF

c     BLACS DESCRIPTOR FOR A AND ITS COPY
      CALL DESCINIT (IDESCA, N, N, MB, NB, 0, 0,
     @               ICTXT, MLOC, IERR)
 

c     BLACS DESCRIPTOR FOR B AND SOLN VECTOR X
      CALL DESCINIT (IDESCB, N, 1, MB, 1, 0, 0, ICTXT, MLOC, IERR)  
 


c      ==============================================
C      FILL ENTRIES OF MATRIX A AND R.H.S. VECTOR B WITH RANDOM ENTRIES
      
      WRITE(*,*)'PROC: ',MYROW, MYCOL,
     @        ' CONSTRUCTING MATRIX A AND RHS VECTOR B ...'
      
      CALL RANDOM_SEED

      CALL RANDOM_SEED ( SIZE = ISEEDSIZE ) ! GET SIZE OF SEED ARRAY

      ALLOCATE ( SEED(1:ISEEDSIZE) )
      CALL RANDOM_SEED ( GET = SEED )
 
      SEED(1) = SEED(1) + NPCOL*MYROW + MYCOL ! ENSURES DIFFERENT SEED
                                              ! FOR EACH PROCESS
      CALL RANDOM_SEED ( PUT = SEED )
  
      CALL RANDOM_NUMBER(LOCAL_B)
       
      CALL RANDOM_NUMBER(LOCAL_A)

c      ==============================================
C      CALL LAPACK LU SOLVER ROUTINE

      WRITE(*,*)'PROC: ',MYROW, MYCOL,
     @    'NOW SOLVING SYSTEM AX = B USING SCALAPACK PDGESV ..'
      ALLOCATE ( IPIV(MLOC + MB), STAT=ISTATUS )
      IF(ISTATUS /= 0) THEN
          WRITE(*,*)'UNABLE TO ALLOCATE IPIV, PROCESS: ',MYROW,MYCOL
          STOP
      ENDIF


      IA = 1
      JA = 1
      IB = 1
      JB = 1
      NRHS = 1
      INFO = 0

      CALL PDGESV(N, NRHS, LOCAL_A, IA, JA, IDESCA, IPIV, 
     @            LOCAL_B, IB, JB, IDESCB, INFO )
 
      IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
          WRITE(*,*)
          WRITE(*,*) 'INFO code returned by PDGESV = ', INFO
          WRITE(*,*)
      END IF
 

c      ==============================================
C     DEALLOCATE MEMORY
      DEALLOCATE(LOCAL_A, STAT=ISTATUS)
      IF(ISTATUS /= 0) THEN
          WRITE(*,*)'UNABLE TO DEALLOCATE ' 
          STOP
      ENDIF   


      DEALLOCATE(LOCAL_B, STAT=ISTATUS)
      IF(ISTATUS /= 0) THEN
          WRITE(*,*)'UNABLE TO DEALLOCATE ' 
          STOP
      ENDIF   


      DEALLOCATE(IPIV, STAT=ISTATUS)
      IF(ISTATUS /= 0) THEN
          WRITE(*,*)'UNABLE TO DEALLOCATE ' 
          STOP
      ENDIF   
    
c     ===================================================
c     RELEASE BLACS CONTEXT

      CALL BLACS_GRIDEXIT(ictxt)
      CALL BLACS_EXIT(0)


      END PROGRAM SOLVE_LU

