C
C
C  Program to demo parallel performance of the following LAIPE subroutine
C  (1) meSolution_CSG_S that solves a system of linear equations by a
C      multiple-entry solver. The left side matrix is symmetric.
C
C  If a multiple-processor computer is available, this subroutine will
C  be implemented in a loop. The loop bound begins with one to the number
C  of system processors, and each loop sets the loop counter as the
C  number of employed processors. Keep computer STANDING ALONE while
C  running this program, and elapsed time with respect to different
C  number of precossors will be reported.
C
C  Running this program, user has to provide the matrix order, i.e., 40000,
C  the lower bandwidth, i.e., 15, and the number of vectors in the right
C  side of equations, i.e., 2. Larger order may have a better performance.
C
C  define global variables
C
      PARAMETER (LIMIT=7000000)
      REAL*4 A(LIMIT)
C
C  variables for memory distribution
C  where IndexA : index to the left side matrix [A]
C        IndexX : index to the right side vector(s)
C        IndexSolution : index to the exact solution that is used to
C                        check accuracy
C        IndexWorking : index to a working space
C
      INTEGER*4 IndexA,IndexX,IndexWorking,IndexSolution
C
C  variable for collecting elapsed time
C
      REAL*4 Second
C
C  variables for calculating accuracy
C
      REAL*4 Lower,Upper,R4TEMP
C
C  variables for problem size
C  where N : matrix order
C        M : lower bandwidth
C        Nset : number of vectors in the right side of equations
C
      INTEGER*4 N,M,Nset
C
C  variables for number of processors
C  where CPUs : number of system processors
C        Processors: number of employed processors
C
      INTEGER*4 Processors,CPUs
C
C  execution flag
C
      LOGICAL*4 NoGood
C
C  numerical zero
C
      REAL*4 ZERO
      DATA ZERO/0.0001/
C
C  enter order of the system
C
      WRITE(*,'('' Enter order (i.e., 40000 or something else): '',$)')
      READ(*,*) N
C
C  enter half bandwidth
C
      WRITE(*,
     &'('' Enter half bandwidth (i.e., 15 or something else): '',$)')
      READ(*,*) M
C
C  enter number of vectors
C
      WRITE(*,
     &'('' Enter number of vectors (i.e., 2 or something else): '',$)')
      READ(*,*) Nset
C
C  memory distribution
C
      IndexX=1
      IndexSolution=IndexX+N*Nset
      IndexA=IndexSolution+N*Nset
      IndexWorking=IndexA+(N-1)*M+N
C
C  check memory space
C
      IF(IndexWorking+2*M*N.GT.LIMIT) THEN
         WRITE(*,'('' Memory overflow'')')
         STOP
      END IF
C
C  get number of system processors
C
      CALL GetCPUs(CPUs)
C
C  set processors to execute the following statements
C
      DO Processors=1,CPUS
         CALL SetEmployedProcessors(Processors)
C
C  generate left side matrix and right side vector
C
         Write(*,*)
         Write(*,'('' Generate data in a pseudo-random procedure...'',
     &             $)')
         CALL GenerateData(A(IndexX),A(IndexA),
     &                     A(IndexSolution),N,M,Nset)
C
C  output number of employed processors
C
         WRITE(*,
     &   '(/,'' Number of employed processors: '',I3)') Processors
C
C  start collecting elapsed time spent in the solution
C
         CALL CollectElapsedTime
C      
C  solve in parallel
C
         CALL meSolution_CSG_S
     &       (
     &        A(IndexA),
     &        N,
     &        M,
     &        A(IndexX),
     &        Nset,
     &        A(IndexWorking),
     &        Zero,
     &        NoGood
     &       )
C
C  output elapsed time
C
         CALL GetElapsedTime(Second)
         Write(*,'('' Solve>>> Elapsed Time (seconds): '',F8.2)')
     &                                                   Second
C
C  stop if the system is not suitable for routine Decompose_DSG
C
         IF(NoGood) THEN
            WRITE(*,'('' The system is not suit for routine'',
     &                '' Decompose_CSG.'')')
            STOP
         END IF
C
C  calculate lower and upper bounds of relative errror
C
         KK=IndexX
         JJ=IndexSolution
         DO II=1,Nset
            Lower=(A(JJ)-A(KK))/A(JJ)
            Upper=Lower
            DO I=2,N
               JJ=JJ+1
               KK=KK+1
               R4TEMP=(A(JJ)-A(KK))/A(JJ)
               IF(R4TEMP.LT.Lower) THEN
                  Lower=R4TEMP
               ELSE IF(R4TEMP.GT.Upper) THEN
                  Upper=R4TEMP
               END IF
            END DO
            JJ=JJ+1
            KK=KK+1
C
C  output lower bound and upper bound of relative error
C
            WRITE(*,*) 'Lower and upper bound of the relative error:'
            WRITE(*,*) Lower,Upper
         END DO
      END DO
      STOP
      END
      SUBROUTINE GenerateData(X,A,Solution,N,M,Nset)
C
C
C  routine to generate the left side matrix, and the right side vector
C  (A)FORTRAN CALL: CALL GenerateData(X,A,Solution,N,M,Nset)
C     1.X: <R4> right side vector, dimension(N,Nset)
C     2.A: <R4> left side matrix, dimension(*)
C     3.Solution: <R4> return the solution for checking accuracy,
C                      dimension(N,Nset)
C     3.N: <I4> order
C     4.M: <I4> half bandwidth
C     5.Nset: <I4> number of vectors
C
      INTEGER*4 N,M,SEED,Nset
      REAL*4 A(M,1),X(N,1),Solution(N,1)
C
C  private variables
C
      REAL*4 R4TEMP$
      AUTOMATIC R4TEMP$
C
C  left side matrix
C
      Seed=456789
      DO I=1,N
         DO J=I,MIN0(I+M,N)
            A(J,I)=RAN(SEED)*100.0
         END DO
      END DO
      DO I=1,N
         DO J=I+1,MIN0(I+M,N)
            IF(A(J,I).GT.A(I,I)) THEN
               R4TEMP$=A(J,I)
               A(J,I)=A(I,I)
               A(I,I)=R4TEMP$
            END IF
         END DO
      END DO
      DO I=1,N
         A(I,I)=A(I,I)*50.0
      END DO
C
C  the solution
C
      DO II=1,Nset
         DO I=1,N
            Solution(I,II)=Ran(Seed)*500.0
         END DO
      END DO
C
C  right side vector
C
      DO II=1,Nset
         DO I=1,N
            X(I,II)=A(I,I)*Solution(I,II)
         END DO
         DO I=1,N
            DO J=I+1,MIN0(N,I+M)
               X(I,II)=X(I,II)+A(J,I)*Solution(J,II)
               X(J,II)=X(J,II)+A(J,I)*Solution(I,II)
            END DO
         END DO
      END DO
C
      RETURN
      END
      REAL*4 FUNCTION RAN(SEED)
C
C
C  function to generate a pseudo-random number
C  (A)FORTRAN ENTRY: RAN(SEED)
C     1.RAN: <R4> return a random number
C     2.SEED: <I4> a seed (updated)
C
      INTEGER*4 SEED,A
      REAL*8 RTEMPA,RTEMPB
      DATA M/2147483647/,A/16807/
C
C  generate a random number
C
      RTEMPA=DFLOAT(SEED)*A
      ITEMP=INT(RTEMPA/DFLOAT(M))
      RTEMPB=DFLOAT(M)*DFLOAT(ITEMP)
      SEED=INT(RTEMPA-RTEMPB)
      RAN=DFLOAT(SEED)/DFLOAT(M)
C
      RETURN
      END
