
! *** generated by SAPFOR with version 1757 and build date: Mar 26 2021 10:17:52

! *** generated by SAPFOR with version 1651 and build date: Oct  5 2020 10:15:03 
!-------------------------------------------------------------------------! 
!                                                                         ! 
!        N  A  S     P A R A L L E L     B E N C H M A R K S  3.3         ! 
!                                                                         ! 
!                                   C G                                   ! 
!                                                                         ! 
!-------------------------------------------------------------------------! 
!                                                                         ! 
!    This benchmark is part of the NAS Parallel Benchmark 3.3 suite.      ! 
!    It is described in NAS Technical Reports 95-020 and 02-007           ! 
!                                                                         ! 
!    Permission to use, copy, distribute and modify this software         ! 
!    for any purpose with or without fee is hereby granted.  We           ! 
!    request, however, that all derived work reference the NAS            ! 
!    Parallel Benchmarks 3.3. This software is provided "as is"           ! 
!    without express or implied warranty.                                 ! 
!                                                                         ! 
!    Information on NPB 3.3, including the technical report, the          ! 
!    original specifications, source code, results and information        ! 
!    on how to submit new results, is available at:                       ! 
!                                                                         ! 
!           http://www.nas.nasa.gov/Software/NPB/                         ! 
!                                                                         ! 
!    Send comments or suggestions to  npb@nas.nasa.gov                    ! 
!                                                                         ! 
!          NAS Parallel Benchmarks Group                                  ! 
!          NASA Ames Research Center                                      ! 
!          Mail Stop: T27A-1                                              ! 
!          Moffett Field, CA   94035-1000                                 ! 
!                                                                         ! 
!          E-mail:  npb@nas.nasa.gov                                      ! 
!          Fax:     (650) 604-3957                                        ! 
!                                                                         ! 
!-------------------------------------------------------------------------! 
!--------------------------------------------------------------------- 
!      
! Authors: M. Yarrow 
!          C. Kuszmaul 
!          R. F. Van der Wijngaart 
!          H. Jin 
!      
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         PROGRAM CG

!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         IMPLICIT NONE

      INCLUDE 'mpinpb.h' 
      INCLUDE 'timing.h' 
         INTEGER ::  STATUS(MPI_STATUS_SIZE),REQUEST,IERR

      INCLUDE 'npbparams.h' 
         INTEGER ::  NUM_PROCS
         PARAMETER (NUM_PROCS = NUM_PROC_COLS * NUM_PROC_ROWS)

!--------------------------------------------------------------------- 
!  Class specific parameters:  
!  It appears here for reference only. 
!  These are their values, however, this info is imported in the npbparams.h 
!  include file, which is written by the sys/setparams.c program. 
!--------------------------------------------------------------------- 
!---------- 
!  Class S: 
!---------- 
!C       parameter( na=1400,  
!C      >           nonzer=7,  
!C      >           shift=10.,  
!C      >           niter=15, 
!C      >           rcond=1.0d-1 ) 
!---------- 
!  Class W: 
!---------- 
!C       parameter( na=7000, 
!C      >           nonzer=8,  
!C      >           shift=12.,  
!C      >           niter=15, 
!C      >           rcond=1.0d-1 ) 
!---------- 
!  Class A: 
!---------- 
!C       parameter( na=14000, 
!C      >           nonzer=11,  
!C      >           shift=20.,  
!C      >           niter=15, 
!C      >           rcond=1.0d-1 ) 
!---------- 
!  Class B: 
!---------- 
!C       parameter( na=75000,  
!C      >           nonzer=13,  
!C      >           shift=60.,  
!C      >           niter=75, 
!C      >           rcond=1.0d-1 ) 
!---------- 
!  Class C: 
!---------- 
!C       parameter( na=150000,  
!C      >           nonzer=15,  
!C      >           shift=110.,  
!C      >           niter=75, 
!C      >           rcond=1.0d-1 ) 
!---------- 
!  Class D: 
!---------- 
!C       parameter( na=1500000,  
!C      >           nonzer=21,  
!C      >           shift=500.,  
!C      >           niter=100, 
!C      >           rcond=1.0d-1 ) 
!---------- 
!  Class E: 
!---------- 
!C       parameter( na=9000000,  
!C      >           nonzer=26,  
!C      >           shift=1500.,  
!C      >           niter=100, 
!C      >           rcond=1.0d-1 ) 
         INTEGER ::  NZ
         PARAMETER (NZ = NA * (NONZER + 1) / NUM_PROCS * (NONZER + 1) + 
     &NONZER + NA * (NONZER + 2 + NUM_PROCS / 256) / NUM_PROC_COLS)
         COMMON /PARTIT_SIZE/NAA,NZZ,NPCOLS,NPROWS,PROC_COL,PROC_ROW,FIR
     &STROW,LASTROW,FIRSTCOL,LASTCOL,EXCH_PROC,EXCH_RECV_LENGTH,SEND_STA
     &RT,SEND_LEN
         INTEGER ::  NAA,NZZ,NPCOLS,NPROWS,PROC_COL,PROC_ROW,FIRSTROW,LA
     &STROW,FIRSTCOL,LASTCOL,EXCH_PROC,EXCH_RECV_LENGTH,SEND_START,SEND_
     &LEN
         COMMON /MAIN_INT_MEM/COLIDX,ROWSTR,IV,AROW,ACOL
         INTEGER ::  COLIDX(NZ),ROWSTR(NA + 1),IV(2 * NA + 1),AROW(NZ),A
     &COL(NZ)
         COMMON /MAIN_FLT_MEM/V,AELT,A,X,Z,P,Q,R,W
         DOUBLE PRECISION ::  V(NA + 1),AELT(NZ),A(NZ),X(NA / NUM_PROC_R
     &OWS + 2),Z(NA / NUM_PROC_ROWS + 2),P(NA / NUM_PROC_ROWS + 2),Q(NA 
     &/ NUM_PROC_ROWS + 2),R(NA / NUM_PROC_ROWS + 2),W(NA / NUM_PROC_ROW
     &S + 2)
         COMMON /URANDO/AMULT,TRAN
         DOUBLE PRECISION ::  AMULT,TRAN
         INTEGER ::  L2NPCOLS
         INTEGER ::  REDUCE_EXCH_PROC(NUM_PROC_COLS)
         INTEGER ::  REDUCE_SEND_STARTS(NUM_PROC_COLS)
         INTEGER ::  REDUCE_SEND_LENGTHS(NUM_PROC_COLS)
         INTEGER ::  REDUCE_RECV_STARTS(NUM_PROC_COLS)
         INTEGER ::  REDUCE_RECV_LENGTHS(NUM_PROC_COLS)
         INTEGER ::  I,J,K,IT
         DOUBLE PRECISION ::  ZETA,RANDLC
         EXTERNAL RANDLC
         DOUBLE PRECISION ::  RNORM
         DOUBLE PRECISION ::  NORM_TEMP1(2),NORM_TEMP2(2)
         DOUBLE PRECISION ::  T,TMAX,MFLOPS
         EXTERNAL TIMER_READ
         DOUBLE PRECISION ::  TIMER_READ
         CHARACTER ::  CLASS
         LOGICAL ::  VERIFIED
         DOUBLE PRECISION ::  ZETA_VERIFY_VALUE,EPSILON,ERR
         DOUBLE PRECISION ::  TSUM(T_LAST + 2),T1(T_LAST + 2),TMING(T_LA
     &ST + 2),TMAXG(T_LAST + 2)
         CHARACTER ::  T_RECS(T_LAST + 2)*8
         DATA T_RECS/'total', 'conjg', 'rcomm', 'ncomm',            ' to
     &tcomp', ' totcomm'/ 
         INTERFACE 
            SUBROUTINE CONJ_GRAD (COLIDX, ROWSTR, X, Z, A, P, Q, R, W, R
     &NORM, L2NPCOLS, REDUCE_EXCH_PROC, REDUCE_SEND_STARTS, REDUCE_SEND_
     &LENGTHS, REDUCE_RECV_STARTS, REDUCE_RECV_LENGTHS)
            IMPLICIT NONE
            COMMON /PARTIT_SIZE/NAA,NZZ,NPCOLS,NPROWS,PROC_COL,PROC_ROW,
     &FIRSTROW,LASTROW,FIRSTCOL,LASTCOL,EXCH_PROC,EXCH_RECV_LENGTH,SEND_
     &START,SEND_LEN
            INTEGER ::  NAA,NZZ,NPCOLS,NPROWS,PROC_COL,PROC_ROW,FIRSTROW
     &,LASTROW,FIRSTCOL,LASTCOL,EXCH_PROC,EXCH_RECV_LENGTH,SEND_START,SE
     &ND_LEN
            DOUBLE PRECISION ::  X(:),Z(:),P(:),Q(:),R(:),W(:)
            DOUBLE PRECISION ::  A(NZZ)
            INTEGER ::  COLIDX(NZZ),ROWSTR(NAA + 1)
            INTEGER ::  L2NPCOLS
            INTEGER ::  REDUCE_EXCH_PROC(L2NPCOLS)
            INTEGER ::  REDUCE_SEND_STARTS(L2NPCOLS)
            INTEGER ::  REDUCE_SEND_LENGTHS(L2NPCOLS)
            INTEGER ::  REDUCE_RECV_STARTS(L2NPCOLS)
            INTEGER ::  REDUCE_RECV_LENGTHS(L2NPCOLS)
            DOUBLE PRECISION ::  D,SUM,RHO,RHO0,ALPHA,BETA,RNORM
            END SUBROUTINE  
         END INTERFACE

!--------------------------------------------------------------------- 
!  Set up mpi initialization and number of proc testing 
!--------------------------------------------------------------------- 
         CALL INITIALIZE_MPI()
         IF (NA .EQ. 1400 .AND. NONZER .EQ. 7 .AND. NITER .EQ. 15 .AND. 
     &SHIFT .EQ. 10.D0) THEN
            CLASS = 'S'
            ZETA_VERIFY_VALUE = 8.5971775078648D0
         ELSE IF (NA .EQ. 7000 .AND. NONZER .EQ. 8 .AND. NITER .EQ. 15 .
     &AND. SHIFT .EQ. 12.D0) THEN  
            CLASS = 'W'
            ZETA_VERIFY_VALUE = 10.362595087124D0
         ELSE IF (NA .EQ. 14000 .AND. NONZER .EQ. 11 .AND. NITER .EQ. 15
     & .AND. SHIFT .EQ. 20.D0) THEN  
            CLASS = 'A'
            ZETA_VERIFY_VALUE = 17.130235054029D0
         ELSE IF (NA .EQ. 75000 .AND. NONZER .EQ. 13 .AND. NITER .EQ. 75
     & .AND. SHIFT .EQ. 60.D0) THEN  
            CLASS = 'B'
            ZETA_VERIFY_VALUE = 22.712745482631D0
         ELSE IF (NA .EQ. 150000 .AND. NONZER .EQ. 15 .AND. NITER .EQ. 7
     &5 .AND. SHIFT .EQ. 110.D0) THEN  
            CLASS = 'C'
            ZETA_VERIFY_VALUE = 28.973605592845D0
         ELSE IF (NA .EQ. 1500000 .AND. NONZER .EQ. 21 .AND. NITER .EQ. 
     &100 .AND. SHIFT .EQ. 500.D0) THEN  
            CLASS = 'D'
            ZETA_VERIFY_VALUE = 52.514532105794D0
         ELSE IF (NA .EQ. 9000000 .AND. NONZER .EQ. 26 .AND. NITER .EQ. 
     &100 .AND. SHIFT .EQ. 1.5D3) THEN  
            CLASS = 'E'
            ZETA_VERIFY_VALUE = 77.522164599383D0
         ELSE  
            CLASS = 'U'
         ENDIF  
         IF (ME .EQ. ROOT) THEN
            WRITE (UNIT = *,FMT = 1000) 
            WRITE (UNIT = *,FMT = 1001) NA
            WRITE (UNIT = *,FMT = 1002) NITER
            WRITE (UNIT = *,FMT = 1003) NPROCS
            WRITE (UNIT = *,FMT = 1004) NONZER
            WRITE (UNIT = *,FMT = 1005) SHIFT
1000              FORMAT(//,' NAS Parallel Benchmarks 3.3 -- CG Benchmar
     &k', /)
1001              FORMAT(' Size: ', I10 ) 
1002              FORMAT(' Iterations: ', I5 ) 
1003              FORMAT(' Number of active processes: ', I5 ) 
1004              FORMAT(' Number of nonzeroes per row: ', I8) 
1005              FORMAT(' Eigenvalue shift: ', E8.3) 
         ENDIF  
         IF (.NOT.(CONVERTDOUBLE)) THEN
            DP_TYPE = MPI_DOUBLE_PRECISION
         ELSE  
            DP_TYPE = MPI_REAL
         ENDIF  
         NAA = NA
         NZZ = NZ

!--------------------------------------------------------------------- 
!  Set up processor info, such as whether sq num of procs, etc 
!--------------------------------------------------------------------- 
         CALL SETUP_PROC_INFO(NUM_PROCS,NUM_PROC_ROWS,NUM_PROC_COLS)

!--------------------------------------------------------------------- 
!  Set up partition's submatrix info: firstcol, lastcol, firstrow, lastrow 
!--------------------------------------------------------------------- 
         CALL SETUP_SUBMATRIX_INFO(L2NPCOLS,REDUCE_EXCH_PROC,REDUCE_SEND
     &_STARTS,REDUCE_SEND_LENGTHS,REDUCE_RECV_STARTS,REDUCE_RECV_LENGTHS
     &)
         DO  I = 1,T_LAST
            CALL TIMER_CLEAR(I)
         ENDDO  

!--------------------------------------------------------------------- 
!  Inialize random number generator 
!--------------------------------------------------------------------- 
         TRAN = 314159265.0D0
         AMULT = 1220703125.0D0
         ZETA = RANDLC (TRAN,AMULT)

!--------------------------------------------------------------------- 
!  Set up partition's sparse random matrix for given class size 
!--------------------------------------------------------------------- 
         CALL MAKEA(NAA,NZZ,A,COLIDX,ROWSTR,NONZER,FIRSTROW,LASTROW,FIRS
     &TCOL,LASTCOL,RCOND,AROW,ACOL,AELT,V,IV,SHIFT)

!--------------------------------------------------------------------- 
!  Note: as a result of the above call to makea: 
!        values of j used in indexing rowstr go from 1 --> lastrow-firstrow+1 
!        values of colidx which are col indexes go from firstcol --> lastcol 
!        So: 
!        Shift the col index vals from actual (firstcol --> lastcol )  
!        to local, i.e., (1 --> lastcol-firstcol+1) 
!--------------------------------------------------------------------- 
!DVM$    GET_ACTUAL (COLIDX,ROWSTR)
         DO  J = 1,LASTROW - FIRSTROW + 1
            DO  K = ROWSTR(J),ROWSTR(J + 1) - 1
               COLIDX(K) = COLIDX(K) - FIRSTCOL + 1
            ENDDO  
         ENDDO  
!DVM$    ACTUAL (COLIDX)

!--------------------------------------------------------------------- 
!  set starting vector to (1, 1, .... 1) 
!--------------------------------------------------------------------- 
!DVM$    REGION 
!DVM$    PARALLEL (I), PRIVATE (I),TIE (X(I))
         DO  I = 1,NA / NUM_PROC_ROWS + 1
            X(I) = 1.0D0
         ENDDO  
!DVM$    END REGION
         ZETA = 0.0D0

!--------------------------------------------------------------------- 
!----> 
!  Do one iteration untimed to init all code and data page tables 
!---->                    (then reinit, start timing, to niter its) 
!--------------------------------------------------------------------- 
         DO  IT = 1,1

!--------------------------------------------------------------------- 
!  The call to the conjugate gradient routine: 
!--------------------------------------------------------------------- 
            CALL CONJ_GRAD(COLIDX,ROWSTR,X,Z,A,P,Q,R,W,RNORM,L2NPCOLS,RE
     &DUCE_EXCH_PROC,REDUCE_SEND_STARTS,REDUCE_SEND_LENGTHS,REDUCE_RECV_
     &STARTS,REDUCE_RECV_LENGTHS)

!--------------------------------------------------------------------- 
!  zeta = shift + 1/(x.z) 
!  So, first: (x.z) 
!  Also, find norm of z 
!  So, first: (z.z) 
!--------------------------------------------------------------------- 
            NORM_TEMP1(1) = 0.0D0
            NORM_TEMP1(2) = 0.0D0
!DVM$       REGION 
!DVM$       PARALLEL (J), PRIVATE (J),TIE (X(J),Z(J)),REDUCTION (SUM (NO
!DVM$&RM_TEMP1))
            DO  J = 1,LASTCOL - FIRSTCOL + 1
               NORM_TEMP1(1) = NORM_TEMP1(1) + X(J) * Z(J)
               NORM_TEMP1(2) = NORM_TEMP1(2) + Z(J) * Z(J)
            ENDDO  
!DVM$       END REGION
            DO  I = 1,L2NPCOLS
               IF (TIMERON)  CALL TIMER_START(T_NCOMM)
               CALL MPI_IRECV(NORM_TEMP2,2,DP_TYPE,REDUCE_EXCH_PROC(I),I
     &,MPI_COMM_WORLD,REQUEST,IERR)
               CALL MPI_SEND(NORM_TEMP1,2,DP_TYPE,REDUCE_EXCH_PROC(I),I,
     &MPI_COMM_WORLD,IERR)
               CALL MPI_WAIT(REQUEST,STATUS,IERR)
               IF (TIMERON)  CALL TIMER_STOP(T_NCOMM)
               NORM_TEMP1(1) = NORM_TEMP1(1) + NORM_TEMP2(1)
               NORM_TEMP1(2) = NORM_TEMP1(2) + NORM_TEMP2(2)
            ENDDO  
!DVM$       GET_ACTUAL (NORM_TEMP1)
            NORM_TEMP1(2) = 1.0D0 / SQRT (NORM_TEMP1(2))
!DVM$       ACTUAL (NORM_TEMP1(2))

!--------------------------------------------------------------------- 
!  Normalize z to obtain x 
!--------------------------------------------------------------------- 
!DVM$       REGION 
!DVM$       PARALLEL (J), PRIVATE (J),TIE (X(J),Z(J))
            DO  J = 1,LASTCOL - FIRSTCOL + 1
               X(J) = NORM_TEMP1(2) * Z(J)
            ENDDO  
!DVM$       END REGION

! end of do one iteration untimed 
         ENDDO  

!--------------------------------------------------------------------- 
!  set starting vector to (1, 1, .... 1) 
!--------------------------------------------------------------------- 
!      
!  NOTE: a questionable limit on size:  should this be na/num_proc_cols+1 ? 
!      
!DVM$    REGION 
!DVM$    PARALLEL (I), PRIVATE (I),TIE (X(I))
         DO  I = 1,NA / NUM_PROC_ROWS + 1
            X(I) = 1.0D0
         ENDDO  
!DVM$    END REGION
         ZETA = 0.0D0

!--------------------------------------------------------------------- 
!  Synchronize and start timing 
!--------------------------------------------------------------------- 
         DO  I = 1,T_LAST
            CALL TIMER_CLEAR(I)
         ENDDO  
         CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
         CALL TIMER_CLEAR(1)
         CALL TIMER_START(1)

!--------------------------------------------------------------------- 
!----> 
!  Main Iteration for inverse power method 
!----> 
!--------------------------------------------------------------------- 
         DO  IT = 1,NITER

!--------------------------------------------------------------------- 
!  The call to the conjugate gradient routine: 
!--------------------------------------------------------------------- 
            CALL CONJ_GRAD(COLIDX,ROWSTR,X,Z,A,P,Q,R,W,RNORM,L2NPCOLS,RE
     &DUCE_EXCH_PROC,REDUCE_SEND_STARTS,REDUCE_SEND_LENGTHS,REDUCE_RECV_
     &STARTS,REDUCE_RECV_LENGTHS)

!--------------------------------------------------------------------- 
!  zeta = shift + 1/(x.z) 
!  So, first: (x.z) 
!  Also, find norm of z 
!  So, first: (z.z) 
!--------------------------------------------------------------------- 
            NORM_TEMP1(1) = 0.0D0
            NORM_TEMP1(2) = 0.0D0
!DVM$       REGION 
!DVM$       PARALLEL (J), PRIVATE (J),TIE (X(J),Z(J)),REDUCTION (SUM (NO
!DVM$&RM_TEMP1))
            DO  J = 1,LASTCOL - FIRSTCOL + 1
               NORM_TEMP1(1) = NORM_TEMP1(1) + X(J) * Z(J)
               NORM_TEMP1(2) = NORM_TEMP1(2) + Z(J) * Z(J)
            ENDDO  
!DVM$       END REGION
            DO  I = 1,L2NPCOLS
               IF (TIMERON)  CALL TIMER_START(T_NCOMM)
               CALL MPI_IRECV(NORM_TEMP2,2,DP_TYPE,REDUCE_EXCH_PROC(I),I
     &,MPI_COMM_WORLD,REQUEST,IERR)
               CALL MPI_SEND(NORM_TEMP1,2,DP_TYPE,REDUCE_EXCH_PROC(I),I,
     &MPI_COMM_WORLD,IERR)
               CALL MPI_WAIT(REQUEST,STATUS,IERR)
               IF (TIMERON)  CALL TIMER_STOP(T_NCOMM)
               NORM_TEMP1(1) = NORM_TEMP1(1) + NORM_TEMP2(1)
               NORM_TEMP1(2) = NORM_TEMP1(2) + NORM_TEMP2(2)
            ENDDO  
!DVM$       GET_ACTUAL (NORM_TEMP1)
            NORM_TEMP1(2) = 1.0D0 / SQRT (NORM_TEMP1(2))
!DVM$       ACTUAL (NORM_TEMP1(2))
            IF (ME .EQ. ROOT) THEN
               ZETA = SHIFT + 1.0D0 / NORM_TEMP1(1)
               IF (IT .EQ. 1)  WRITE (UNIT = *,FMT = 9000) 
               WRITE (UNIT = *,FMT = 9001) IT,RNORM,ZETA
            ENDIF  
9000              FORMAT( /,'   iteration           ||r||               
     &  zeta' ) 
9001              FORMAT( 4X, I5, 7X, E20.14, F20.13 ) 

!--------------------------------------------------------------------- 
!  Normalize z to obtain x 
!--------------------------------------------------------------------- 
!DVM$       REGION 
!DVM$       PARALLEL (J), PRIVATE (J),TIE (X(J),Z(J))
            DO  J = 1,LASTCOL - FIRSTCOL + 1
               X(J) = NORM_TEMP1(2) * Z(J)
            ENDDO  
!DVM$       END REGION

! end of main iter inv pow meth 
         ENDDO  
         CALL TIMER_STOP(1)

!--------------------------------------------------------------------- 
!  End of timed section 
!--------------------------------------------------------------------- 
         T = TIMER_READ (1)
         CALL MPI_REDUCE(T,TMAX,1,DP_TYPE,MPI_MAX,ROOT,MPI_COMM_WORLD,IE
     &RR)
         IF (ME .EQ. ROOT) THEN
            WRITE (UNIT = *,FMT = 100) 
100                  FORMAT(' Benchmark completed ') 
            EPSILON = 1.D-10
            IF (CLASS .NE. 'U') THEN
               ERR = ABS (ZETA - ZETA_VERIFY_VALUE) / ZETA_VERIFY_VALUE
               IF (ERR .LE. EPSILON) THEN
                  VERIFIED = .TRUE.
                  WRITE (UNIT = *,FMT = 200) 
                  WRITE (UNIT = *,FMT = 201) ZETA
                  WRITE (UNIT = *,FMT = 202) ERR
200                                          FORMAT(' VERIFICATION SUCCE
     &SSFUL ') 
201                                          FORMAT(' Zeta is    ', E20.
     &13) 
202                                          FORMAT(' Error is   ', E20.
     &13) 
               ELSE  
                  VERIFIED = .FALSE.
                  WRITE (UNIT = *,FMT = 300) 
                  WRITE (UNIT = *,FMT = 301) ZETA
                  WRITE (UNIT = *,FMT = 302) ZETA_VERIFY_VALUE
300                                          FORMAT(' VERIFICATION FAILE
     &D') 
301                                          FORMAT(' Zeta              
     &  ', E20.13)
302                                          FORMAT(' The correct zeta i
     &s ', E20.13)
               ENDIF  
            ELSE  
               VERIFIED = .FALSE.
               WRITE (UNIT = *,FMT = 400) 
               WRITE (UNIT = *,FMT = 401) 
               WRITE (UNIT = *,FMT = 201) ZETA
400                              FORMAT(' Problem size unknown') 
401                              FORMAT(' NO VERIFICATION PERFORMED') 
            ENDIF  
            IF (TMAX .NE. 0.) THEN
               MFLOPS = FLOAT (2 * 75 * 150000) * (3. + FLOAT (15 * (15 
     &+ 1)) + 25. * (5. + FLOAT (15 * (15 + 1))) + 3.) / TMAX / 1000000.
     &0
            ELSE  
               MFLOPS = 0.0
            ENDIF  
            CALL PRINT_RESULTS('CG',CLASS,NA,0,0,NITER,NNODES_COMPILED,N
     &PROCS,TMAX,MFLOPS,'          floating point',VERIFIED,NPBVERSION,C
     &OMPILETIME,CS1,CS2,CS3,CS4,CS5,CS6,CS7)
         ENDIF  
         IF (.NOT.(TIMERON))  GOTO 999
!DVM$    GET_ACTUAL (T1)
         DO  I = 1,T_LAST
            T1(I) = TIMER_READ (I)
         ENDDO  
!DVM$    ACTUAL (T1)
         T1(T_CONJG) = T1(T_CONJG) - T1(T_RCOMM)
!DVM$    ACTUAL (T1(T_CONJG))
         T1(T_LAST + 2) = T1(T_RCOMM) + T1(T_NCOMM)
!DVM$    ACTUAL (T1(T_LAST + 2))
         T1(T_LAST + 1) = T1(T_TOTAL) - T1(T_LAST + 2)
!DVM$    ACTUAL (T1(T_LAST + 1))
         CALL MPI_REDUCE(T1,TSUM,4 + 2,DP_TYPE,MPI_SUM,0,MPI_COMM_WORLD,
     &IERR)
!DVM$    GET_ACTUAL (T1)
         CALL MPI_REDUCE(T1,TMING,4 + 2,DP_TYPE,MPI_MIN,0,MPI_COMM_WORLD
     &,IERR)
!DVM$    GET_ACTUAL (T1)
         CALL MPI_REDUCE(T1,TMAXG,4 + 2,DP_TYPE,MPI_MAX,0,MPI_COMM_WORLD
     &,IERR)
         IF (ME .EQ. 0) THEN
            WRITE (UNIT = *,FMT = 800) NPROCS
!DVM$       GET_ACTUAL (T_RECS,TMAXG,TMING,TSUM)
            DO  I = 1,T_LAST + 2
               TSUM(I) = TSUM(I) / NPROCS
               WRITE (UNIT = *,FMT = 810) I,T_RECS(I),TMING(I),TMAXG(I),
     &TSUM(I)
            ENDDO  
!DVM$       ACTUAL (TSUM)
         ENDIF  
800      FORMAT(' nprocs =', I6, 11X, 'minimum', 5X, 'maximum',        5
     &X, 'average') 
810      FORMAT(' timer ', I2, '(', A8, ') :', 3(2X,F10.4)) 
999      CONTINUE
         CALL MPI_FINALIZE(IERR)

! end main 
         END


!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         SUBROUTINE INITIALIZE_MPI ()

!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         IMPLICIT NONE

      INCLUDE 'mpinpb.h' 
      INCLUDE 'timing.h' 
         INTEGER ::  IERR,FSTATUS
         CALL MPI_INIT(IERR)
         CALL MPI_COMM_RANK(MPI_COMM_WORLD,ME,IERR)
         CALL MPI_COMM_SIZE(MPI_COMM_WORLD,NPROCS,IERR)
         ROOT = 0
         IF (ME .EQ. ROOT) THEN
            OPEN (UNIT = 2,FILE = 'timer.flag',STATUS = 'old',IOSTAT = F
     &STATUS)
            TIMERON = .FALSE.
            IF (FSTATUS .EQ. 0) THEN
               TIMERON = .TRUE.
               CLOSE (UNIT = 2)
            ENDIF  
         ENDIF  
         CALL MPI_BCAST(TIMERON,1,MPI_LOGICAL,0,MPI_COMM_WORLD,IERR)
         RETURN 
         END


!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         SUBROUTINE SETUP_PROC_INFO (NUM_PROCS, NUM_PROC_ROWS, NUM_PROC_
     &COLS)

!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         IMPLICIT NONE

      INCLUDE 'mpinpb.h' 
         COMMON /PARTIT_SIZE/NAA,NZZ,NPCOLS,NPROWS,PROC_COL,PROC_ROW,FIR
     &STROW,LASTROW,FIRSTCOL,LASTCOL,EXCH_PROC,EXCH_RECV_LENGTH,SEND_STA
     &RT,SEND_LEN
         INTEGER ::  NAA,NZZ,NPCOLS,NPROWS,PROC_COL,PROC_ROW,FIRSTROW,LA
     &STROW,FIRSTCOL,LASTCOL,EXCH_PROC,EXCH_RECV_LENGTH,SEND_START,SEND_
     &LEN
         INTEGER ::  NUM_PROCS,NUM_PROC_COLS,NUM_PROC_ROWS
         INTEGER ::  I,IERR
         INTEGER ::  LOG2NPROCS
         INTENT(IN) NUM_PROC_COLS,NUM_PROC_ROWS,NUM_PROCS

!--------------------------------------------------------------------- 
!  num_procs must be a power of 2, and num_procs=num_proc_cols*num_proc_rows 
!  When num_procs is not square, then num_proc_cols = 2*num_proc_rows 
!--------------------------------------------------------------------- 
!  First, number of procs must be power of two.  
!--------------------------------------------------------------------- 
         IF (NPROCS .NE. NUM_PROCS) THEN
            IF (ME .EQ. ROOT)  WRITE (UNIT = *,FMT = 9000) NPROCS,NUM_PR
     &OCS
9000                  FORMAT(      /,'Error: ',/,'num of procs allocated
     &   (',                  I4, ' )',                 /,'is not equal 
     &to',/,                 'compiled number of procs (',              
     &   I4, ' )',/   ) 
            CALL MPI_FINALIZE(IERR)
            STOP
         ENDIF  
         I = NUM_PROC_COLS
100      CONTINUE
         IF (I .NE. 1 .AND. I / 2 * 2 .NE. I) THEN
            IF (ME .EQ. ROOT) THEN
               WRITE (UNIT = *,FMT = *) 'Error: num_proc_cols is ',NUM_P
     &ROC_COLS,' which is not a power of two'
            ENDIF  
            CALL MPI_FINALIZE(IERR)
            STOP
         ENDIF  
         I = I / 2
         IF (I .NE. 0) THEN
            GOTO 100
         ENDIF  
         I = NUM_PROC_ROWS
200      CONTINUE
         IF (I .NE. 1 .AND. I / 2 * 2 .NE. I) THEN
            IF (ME .EQ. ROOT) THEN
               WRITE (UNIT = *,FMT = *) 'Error: num_proc_rows is ',NUM_P
     &ROC_ROWS,' which is not a power of two'
            ENDIF  
            CALL MPI_FINALIZE(IERR)
            STOP
         ENDIF  
         I = I / 2
         IF (I .NE. 0) THEN
            GOTO 200
         ENDIF  
         LOG2NPROCS = 0
         I = NPROCS
300      CONTINUE
         IF (I .NE. 1 .AND. I / 2 * 2 .NE. I) THEN
            WRITE (UNIT = *,FMT = *) 'Error: nprocs is ',NPROCS,' which 
     &is not a power of two'
            CALL MPI_FINALIZE(IERR)
            STOP
         ENDIF  
         I = I / 2
         IF (I .NE. 0) THEN
            LOG2NPROCS = LOG2NPROCS + 1
            GOTO 300
         ENDIF  

!C       write( *,* ) 'nprocs, log2nprocs: ',nprocs,log2nprocs 
         NPCOLS = NUM_PROC_COLS
         NPROWS = NUM_PROC_ROWS
         RETURN 
         END


!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         SUBROUTINE SETUP_SUBMATRIX_INFO (L2NPCOLS, REDUCE_EXCH_PROC, RE
     &DUCE_SEND_STARTS, REDUCE_SEND_LENGTHS, REDUCE_RECV_STARTS, REDUCE_
     &RECV_LENGTHS)

!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         IMPLICIT NONE

      INCLUDE 'mpinpb.h' 
      INCLUDE 'npbparams.h' 
         INTEGER ::  COL_SIZE,ROW_SIZE
         COMMON /PARTIT_SIZE/NAA,NZZ,NPCOLS,NPROWS,PROC_COL,PROC_ROW,FIR
     &STROW,LASTROW,FIRSTCOL,LASTCOL,EXCH_PROC,EXCH_RECV_LENGTH,SEND_STA
     &RT,SEND_LEN
         INTEGER ::  NAA,NZZ,NPCOLS,NPROWS,PROC_COL,PROC_ROW,FIRSTROW,LA
     &STROW,FIRSTCOL,LASTCOL,EXCH_PROC,EXCH_RECV_LENGTH,SEND_START,SEND_
     &LEN
         INTEGER ::  REDUCE_EXCH_PROC(NUM_PROC_COLS)
         INTEGER ::  REDUCE_SEND_STARTS(NUM_PROC_COLS)
         INTEGER ::  REDUCE_SEND_LENGTHS(NUM_PROC_COLS)
         INTEGER ::  REDUCE_RECV_STARTS(NUM_PROC_COLS)
         INTEGER ::  REDUCE_RECV_LENGTHS(NUM_PROC_COLS)
         INTEGER ::  I,J
         INTEGER ::  DIV_FACTOR
         INTEGER ::  L2NPCOLS
         INTENT(INOUT) L2NPCOLS
         INTENT(OUT) REDUCE_RECV_LENGTHS,REDUCE_RECV_STARTS,REDUCE_SEND_
     &LENGTHS,REDUCE_SEND_STARTS,REDUCE_EXCH_PROC
         PROC_ROW = ME / NPCOLS
         PROC_COL = ME - PROC_ROW * NPCOLS

!--------------------------------------------------------------------- 
!  If naa evenly divisible by npcols, then it is evenly divisible  
!  by nprows  
!--------------------------------------------------------------------- 
         IF (NAA / NPCOLS * NPCOLS .EQ. NAA) THEN
            COL_SIZE = NAA / NPCOLS
            FIRSTCOL = PROC_COL * COL_SIZE + 1
            LASTCOL = FIRSTCOL - 1 + COL_SIZE
            ROW_SIZE = NAA / NPROWS
            FIRSTROW = PROC_ROW * ROW_SIZE + 1
            LASTROW = FIRSTROW - 1 + ROW_SIZE

!--------------------------------------------------------------------- 
!  If naa not evenly divisible by npcols, then first subdivide for nprows 
!  and then, if npcols not equal to nprows (i.e., not a sq number of procs),  
!  get col subdivisions by dividing by 2 each row subdivision. 
!--------------------------------------------------------------------- 
         ELSE  
            IF (PROC_ROW .LT. NAA - NAA / NPROWS * NPROWS) THEN
               ROW_SIZE = NAA / NPROWS + 1
               FIRSTROW = PROC_ROW * ROW_SIZE + 1
               LASTROW = FIRSTROW - 1 + ROW_SIZE
            ELSE  
               ROW_SIZE = NAA / NPROWS
               FIRSTROW = (NAA - NAA / NPROWS * NPROWS) * (ROW_SIZE + 1)
     & + (PROC_ROW - (NAA - NAA / NPROWS * NPROWS)) * ROW_SIZE + 1
               LASTROW = FIRSTROW - 1 + ROW_SIZE
            ENDIF  
            IF (NPCOLS .EQ. NPROWS) THEN
               IF (PROC_COL .LT. NAA - NAA / NPCOLS * NPCOLS) THEN
                  COL_SIZE = NAA / NPCOLS + 1
                  FIRSTCOL = PROC_COL * COL_SIZE + 1
                  LASTCOL = FIRSTCOL - 1 + COL_SIZE
               ELSE  
                  COL_SIZE = NAA / NPCOLS
                  FIRSTCOL = (NAA - NAA / NPCOLS * NPCOLS) * (COL_SIZE +
     & 1) + (PROC_COL - (NAA - NAA / NPCOLS * NPCOLS)) * COL_SIZE + 1
                  LASTCOL = FIRSTCOL - 1 + COL_SIZE
               ENDIF  
            ELSE  
               IF (PROC_COL / 2 .LT. NAA - NAA / (NPCOLS / 2) * (NPCOLS 
     &/ 2)) THEN
                  COL_SIZE = NAA / (NPCOLS / 2) + 1
                  FIRSTCOL = PROC_COL / 2 * COL_SIZE + 1
                  LASTCOL = FIRSTCOL - 1 + COL_SIZE
               ELSE  
                  COL_SIZE = NAA / (NPCOLS / 2)
                  FIRSTCOL = (NAA - NAA / (NPCOLS / 2) * (NPCOLS / 2)) *
     & (COL_SIZE + 1) + (PROC_COL / 2 - (NAA - NAA / (NPCOLS / 2) * (NPC
     &OLS / 2))) * COL_SIZE + 1
                  LASTCOL = FIRSTCOL - 1 + COL_SIZE
               ENDIF  

!C               write( *,* ) col_size,firstcol,lastcol 
               IF (MOD (ME,2) .EQ. 0) THEN
                  LASTCOL = FIRSTCOL - 1 + (COL_SIZE - 1) / 2 + 1
               ELSE  
                  FIRSTCOL = FIRSTCOL + (COL_SIZE - 1) / 2 + 1
                  LASTCOL = FIRSTCOL - 1 + COL_SIZE / 2

!C                   write( *,* ) firstcol,lastcol 
               ENDIF  
            ENDIF  
         ENDIF  
         IF (NPCOLS .EQ. NPROWS) THEN
            SEND_START = 1
            SEND_LEN = LASTROW - FIRSTROW + 1
         ELSE  
            IF (MOD (ME,2) .EQ. 0) THEN
               SEND_START = 1
               SEND_LEN = (1 + LASTROW - FIRSTROW + 1) / 2
            ELSE  
               SEND_START = (1 + LASTROW - FIRSTROW + 1) / 2 + 1
               SEND_LEN = (LASTROW - FIRSTROW + 1) / 2
            ENDIF  
         ENDIF  

!--------------------------------------------------------------------- 
!  Transpose exchange processor 
!--------------------------------------------------------------------- 
         IF (NPCOLS .EQ. NPROWS) THEN
            EXCH_PROC = MOD (ME,NPROWS) * NPROWS + ME / NPROWS
         ELSE  
            EXCH_PROC = 2 * (MOD (ME / 2,NPROWS) * NPROWS + ME / 2 / NPR
     &OWS) + MOD (ME,2)
         ENDIF  
         I = NPCOLS / 2
         L2NPCOLS = 0
         DO WHILE (I .GT. 0)
            L2NPCOLS = L2NPCOLS + 1
            I = I / 2
         ENDDO  

!--------------------------------------------------------------------- 
!  Set up the reduce phase schedules... 
!--------------------------------------------------------------------- 
         DIV_FACTOR = NPCOLS
!DVM$    GET_ACTUAL (REDUCE_EXCH_PROC)
         DO  I = 1,L2NPCOLS
            J = MOD (PROC_COL + DIV_FACTOR / 2,DIV_FACTOR) + PROC_COL / 
     &DIV_FACTOR * DIV_FACTOR
            REDUCE_EXCH_PROC(I) = PROC_ROW * NPCOLS + J
            DIV_FACTOR = DIV_FACTOR / 2
         ENDDO  
!DVM$    ACTUAL (REDUCE_EXCH_PROC)
!DVM$    GET_ACTUAL (REDUCE_RECV_LENGTHS,REDUCE_RECV_STARTS,REDUCE_SEND_
!DVM$&LENGTHS,REDUCE_SEND_STARTS)
         DO  I = L2NPCOLS,1,(-(1))
            IF (NPROWS .EQ. NPCOLS) THEN
               REDUCE_SEND_STARTS(I) = SEND_START
               REDUCE_SEND_LENGTHS(I) = SEND_LEN
               REDUCE_RECV_LENGTHS(I) = LASTROW - FIRSTROW + 1
            ELSE  
               REDUCE_RECV_LENGTHS(I) = SEND_LEN
               IF (I .EQ. L2NPCOLS) THEN
                  REDUCE_SEND_LENGTHS(I) = LASTROW - FIRSTROW + 1 - SEND
     &_LEN
                  IF (ME / 2 * 2 .EQ. ME) THEN
                     REDUCE_SEND_STARTS(I) = SEND_START + SEND_LEN
                  ELSE  
                     REDUCE_SEND_STARTS(I) = 1
                  ENDIF  
               ELSE  
                  REDUCE_SEND_LENGTHS(I) = SEND_LEN
                  REDUCE_SEND_STARTS(I) = SEND_START
               ENDIF  
            ENDIF  
            REDUCE_RECV_STARTS(I) = SEND_START
         ENDDO  
!DVM$    ACTUAL (REDUCE_RECV_LENGTHS,REDUCE_RECV_STARTS,REDUCE_SEND_LENG
!DVM$&THS,REDUCE_SEND_STARTS)
         EXCH_RECV_LENGTH = LASTCOL - FIRSTCOL + 1
         RETURN 
         END


!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         SUBROUTINE CONJ_GRAD (COLIDX, ROWSTR, X, Z, A, P, Q, R, W, RNOR
     &M, L2NPCOLS, REDUCE_EXCH_PROC, REDUCE_SEND_STARTS, REDUCE_SEND_LEN
     &GTHS, REDUCE_RECV_STARTS, REDUCE_RECV_LENGTHS)

!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!  Floaging point arrays here are named as in NPB1 spec discussion of  
!  CG algorithm 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         IMPLICIT NONE

      INCLUDE 'mpinpb.h' 
      INCLUDE 'timing.h' 
         INTEGER ::  STATUS(MPI_STATUS_SIZE),REQUEST
         COMMON /PARTIT_SIZE/NAA,NZZ,NPCOLS,NPROWS,PROC_COL,PROC_ROW,FIR
     &STROW,LASTROW,FIRSTCOL,LASTCOL,EXCH_PROC,EXCH_RECV_LENGTH,SEND_STA
     &RT,SEND_LEN
         INTEGER ::  NAA,NZZ,NPCOLS,NPROWS,PROC_COL,PROC_ROW,FIRSTROW,LA
     &STROW,FIRSTCOL,LASTCOL,EXCH_PROC,EXCH_RECV_LENGTH,SEND_START,SEND_
     &LEN
         DOUBLE PRECISION ::  X(:),Z(:),A(NZZ)
         INTEGER ::  COLIDX(NZZ),ROWSTR(NAA + 1)

! used as work temporary 
         DOUBLE PRECISION ::  P(:),Q(:),R(:),W(:)
         INTEGER ::  L2NPCOLS
         INTEGER ::  REDUCE_EXCH_PROC(L2NPCOLS)
         INTEGER ::  REDUCE_SEND_STARTS(L2NPCOLS)
         INTEGER ::  REDUCE_SEND_LENGTHS(L2NPCOLS)
         INTEGER ::  REDUCE_RECV_STARTS(L2NPCOLS)
         INTEGER ::  REDUCE_RECV_LENGTHS(L2NPCOLS)
         INTEGER ::  I,J,K,IERR
         INTEGER ::  CGIT,CGITMAX
         DOUBLE PRECISION ::  D,SUM,RHO,RHO0,ALPHA,BETA,RNORM
         EXTERNAL TIMER_READ
         DOUBLE PRECISION ::  TIMER_READ
         DATA      CGITMAX / 25 / 
         INTENT(INOUT) W,R,Q,P,Z
         INTENT(IN) REDUCE_RECV_LENGTHS,REDUCE_RECV_STARTS,REDUCE_SEND_L
     &ENGTHS,REDUCE_SEND_STARTS,REDUCE_EXCH_PROC,L2NPCOLS,A,X,ROWSTR,COL
     &IDX
         INTENT(OUT) RNORM
         IF (TIMERON)  CALL TIMER_START(T_CONJG)

!--------------------------------------------------------------------- 
!  Initialize the CG algorithm: 
!--------------------------------------------------------------------- 
!DVM$    REGION 
!DVM$    PARALLEL (J), PRIVATE (J),TIE (P(J),Q(J),R(J),W(J),X(J),Z(J))
         DO  J = 1,NAA / NPROWS + 1
            Q(J) = 0.0D0
            Z(J) = 0.0D0
            R(J) = X(J)
            P(J) = R(J)
            W(J) = 0.0D0
         ENDDO  
!DVM$    END REGION

!--------------------------------------------------------------------- 
!  rho = r.r 
!  Now, obtain the norm of r: First, sum squares of r elements locally... 
!--------------------------------------------------------------------- 
         SUM = 0.0D0
!DVM$    REGION 
!DVM$    PARALLEL (J), PRIVATE (J),TIE (R(J)),REDUCTION (SUM (SUM))
         DO  J = 1,LASTCOL - FIRSTCOL + 1
            SUM = SUM + R(J) * R(J)
         ENDDO  
!DVM$    END REGION

!--------------------------------------------------------------------- 
!  Exchange and sum with procs identified in reduce_exch_proc 
!  (This is equivalent to mpi_allreduce.) 
!  Sum the partial sums of rho, leaving rho on all processors 
!--------------------------------------------------------------------- 
         DO  I = 1,L2NPCOLS
            IF (TIMERON)  CALL TIMER_START(T_RCOMM)
            CALL MPI_IRECV(RHO,1,DP_TYPE,REDUCE_EXCH_PROC(I),I,MPI_COMM_
     &WORLD,REQUEST,IERR)
            CALL MPI_SEND(SUM,1,DP_TYPE,REDUCE_EXCH_PROC(I),I,MPI_COMM_W
     &ORLD,IERR)
            CALL MPI_WAIT(REQUEST,STATUS,IERR)
            IF (TIMERON)  CALL TIMER_STOP(T_RCOMM)
            SUM = SUM + RHO
         ENDDO  
         RHO = SUM

!--------------------------------------------------------------------- 
!----> 
!  The conj grad iteration loop 
!----> 
!--------------------------------------------------------------------- 
         DO  CGIT = 1,CGITMAX

!--------------------------------------------------------------------- 
!  q = A.p 
!  The partition submatrix-vector multiply: use workspace w 
!--------------------------------------------------------------------- 
!DVM$       REGION 
!DVM$       PARALLEL (J), PRIVATE (J,K,SUM),TIE (W(J))
            DO  J = 1,LASTROW - FIRSTROW + 1
               SUM = 0.D0
               DO  K = ROWSTR(J),ROWSTR(J + 1) - 1
                  SUM = SUM + A(K) * P(COLIDX(K))
               ENDDO  
               W(J) = SUM
            ENDDO  
!DVM$       END REGION

!--------------------------------------------------------------------- 
!  Sum the partition submatrix-vec A.p's across rows 
!  Exchange and sum piece of w with procs identified in reduce_exch_proc 
!--------------------------------------------------------------------- 
            DO  I = L2NPCOLS,1,(-(1))
               IF (TIMERON)  CALL TIMER_START(T_RCOMM)

               CALL MPI_IRECV(Q(REDUCE_RECV_STARTS(I)),REDUCE_RECV_LENGT
     &HS(I),DP_TYPE,REDUCE_EXCH_PROC(I),I,MPI_COMM_WORLD,REQUEST,IERR)
!DVM$          ACTUAL (Q(REDUCE_RECV_STARTS(I):REDUCE_RECV_STARTS(I)+
!DVM$& REDUCE_RECV_LENGTHS(I)))
!DVM$          GET_ACTUAL (W(REDUCE_SEND_STARTS(I):REDUCE_SEND_STARTS(I)
!DVM$& +REDUCE_SEND_LENGTHS(I)))
               CALL MPI_SEND(W(REDUCE_SEND_STARTS(I)),REDUCE_SEND_LENGTH
     &S(I),DP_TYPE,REDUCE_EXCH_PROC(I),I,MPI_COMM_WORLD,IERR)
               CALL MPI_WAIT(REQUEST,STATUS,IERR)
               IF (TIMERON)  CALL TIMER_STOP(T_RCOMM)
!DVM$          REGION 
!DVM$          PARALLEL (J), PRIVATE (J),TIE (Q(J),W(J))
               DO  J = SEND_START,SEND_START + REDUCE_RECV_LENGTHS(I) - 
     &1
                  W(J) = W(J) + Q(J)
               ENDDO  
!DVM$          END REGION
            ENDDO  

!--------------------------------------------------------------------- 
!  Exchange piece of q with transpose processor: 
!--------------------------------------------------------------------- 
            IF (L2NPCOLS .NE. 0) THEN
               IF (TIMERON)  CALL TIMER_START(T_RCOMM)

               CALL MPI_IRECV(Q,EXCH_RECV_LENGTH,DP_TYPE,EXCH_PROC,1,MPI
     &_COMM_WORLD,REQUEST,IERR)
!DVM$          ACTUAL (Q(1:EXCH_RECV_LENGTH))
!DVM$          GET_ACTUAL (W(SEND_START:SEND_START+SEND_LEN))
               CALL MPI_SEND(W(SEND_START),SEND_LEN,DP_TYPE,EXCH_PROC,1,
     &MPI_COMM_WORLD,IERR)
               CALL MPI_WAIT(REQUEST,STATUS,IERR)
               IF (TIMERON)  CALL TIMER_STOP(T_RCOMM)
            ELSE  
!DVM$          REGION 
!DVM$          PARALLEL (J), PRIVATE (J),TIE (Q(J),W(J))
               DO  J = 1,EXCH_RECV_LENGTH
                  Q(J) = W(J)
               ENDDO  
!DVM$          END REGION
            ENDIF  

!--------------------------------------------------------------------- 
!  Clear w for reuse... 
!--------------------------------------------------------------------- 
!DVM$       REGION 
!DVM$       PARALLEL (J), PRIVATE (J),TIE (W(J))
            DO  J = 1,MAX (LASTROW - FIRSTROW + 1,LASTCOL - FIRSTCOL + 1
     &)
               W(J) = 0.0D0
            ENDDO  
!DVM$       END REGION

!--------------------------------------------------------------------- 
!  Obtain p.q 
!--------------------------------------------------------------------- 
            SUM = 0.0D0
!DVM$       REGION 
!DVM$       PARALLEL (J), PRIVATE (J),TIE (P(J),Q(J)),REDUCTION (SUM (SU
!DVM$&M))
            DO  J = 1,LASTCOL - FIRSTCOL + 1
               SUM = SUM + P(J) * Q(J)
            ENDDO  
!DVM$       END REGION

!--------------------------------------------------------------------- 
!  Obtain d with a sum-reduce 
!--------------------------------------------------------------------- 
            DO  I = 1,L2NPCOLS
               IF (TIMERON)  CALL TIMER_START(T_RCOMM)
               CALL MPI_IRECV(D,1,DP_TYPE,REDUCE_EXCH_PROC(I),I,MPI_COMM
     &_WORLD,REQUEST,IERR)
               CALL MPI_SEND(SUM,1,DP_TYPE,REDUCE_EXCH_PROC(I),I,MPI_COM
     &M_WORLD,IERR)
               CALL MPI_WAIT(REQUEST,STATUS,IERR)
               IF (TIMERON)  CALL TIMER_STOP(T_RCOMM)
               SUM = SUM + D
            ENDDO  
            D = SUM

!--------------------------------------------------------------------- 
!  Obtain alpha = rho / (p.q) 
!--------------------------------------------------------------------- 
            ALPHA = RHO / D

!--------------------------------------------------------------------- 
!  Save a temporary of rho 
!--------------------------------------------------------------------- 
            RHO0 = RHO

!--------------------------------------------------------------------- 
!  Obtain z = z + alpha*p 
!  and    r = r - alpha*q 
!--------------------------------------------------------------------- 
!DVM$       REGION 
!DVM$       PARALLEL (J), PRIVATE (J),TIE (P(J),Q(J),R(J),Z(J))
            DO  J = 1,LASTCOL - FIRSTCOL + 1
               Z(J) = Z(J) + ALPHA * P(J)
               R(J) = R(J) - ALPHA * Q(J)
            ENDDO  
!DVM$       END REGION

!--------------------------------------------------------------------- 
!  rho = r.r 
!  Now, obtain the norm of r: First, sum squares of r elements locally... 
!--------------------------------------------------------------------- 
            SUM = 0.0D0
!DVM$       REGION 
!DVM$       PARALLEL (J), PRIVATE (J),TIE (R(J)),REDUCTION (SUM (SUM))
            DO  J = 1,LASTCOL - FIRSTCOL + 1
               SUM = SUM + R(J) * R(J)
            ENDDO  
!DVM$       END REGION

!--------------------------------------------------------------------- 
!  Obtain rho with a sum-reduce 
!--------------------------------------------------------------------- 
            DO  I = 1,L2NPCOLS
               IF (TIMERON)  CALL TIMER_START(T_RCOMM)
               CALL MPI_IRECV(RHO,1,DP_TYPE,REDUCE_EXCH_PROC(I),I,MPI_CO
     &MM_WORLD,REQUEST,IERR)
               CALL MPI_SEND(SUM,1,DP_TYPE,REDUCE_EXCH_PROC(I),I,MPI_COM
     &M_WORLD,IERR)
               CALL MPI_WAIT(REQUEST,STATUS,IERR)
               IF (TIMERON)  CALL TIMER_STOP(T_RCOMM)
               SUM = SUM + RHO
            ENDDO  
            RHO = SUM

!--------------------------------------------------------------------- 
!  Obtain beta: 
!--------------------------------------------------------------------- 
            BETA = RHO / RHO0

!--------------------------------------------------------------------- 
!  p = r + beta*p 
!--------------------------------------------------------------------- 
!DVM$       REGION 
!DVM$       PARALLEL (J), PRIVATE (J),TIE (P(J),R(J))
            DO  J = 1,LASTCOL - FIRSTCOL + 1
               P(J) = R(J) + BETA * P(J)
            ENDDO  
!DVM$       END REGION

! end of do cgit=1,cgitmax 
         ENDDO  

!--------------------------------------------------------------------- 
!  Compute residual norm explicitly:  ||r|| = ||x - A.z|| 
!  First, form A.z 
!  The partition submatrix-vector multiply 
!--------------------------------------------------------------------- 
!DVM$    REGION 
!DVM$    PARALLEL (J), PRIVATE (J,K,SUM),TIE (W(J))
         DO  J = 1,LASTROW - FIRSTROW + 1
            SUM = 0.D0
            DO  K = ROWSTR(J),ROWSTR(J + 1) - 1
               SUM = SUM + A(K) * Z(COLIDX(K))
            ENDDO  
            W(J) = SUM
         ENDDO  
!DVM$    END REGION

!--------------------------------------------------------------------- 
!  Sum the partition submatrix-vec A.z's across rows 
!--------------------------------------------------------------------- 
         DO  I = L2NPCOLS,1,(-(1))
            IF (TIMERON)  CALL TIMER_START(T_RCOMM)

            CALL MPI_IRECV(R(REDUCE_RECV_STARTS(I)),REDUCE_RECV_LENGTHS(
     &I),DP_TYPE,REDUCE_EXCH_PROC(I),I,MPI_COMM_WORLD,REQUEST,IERR)
!DVM$       ACTUAL (R(REDUCE_RECV_STARTS(I):REDUCE_RECV_STARTS(I)+
!DVM$& REDUCE_RECV_LENGTHS(I)))

!DVM$       GET_ACTUAL (W(REDUCE_SEND_STARTS(I):REDUCE_SEND_STARTS(I)+
!DVM$& REDUCE_SEND_LENGTHS(I)))
            CALL MPI_SEND(W(REDUCE_SEND_STARTS(I)),REDUCE_SEND_LENGTHS(I
     &),DP_TYPE,REDUCE_EXCH_PROC(I),I,MPI_COMM_WORLD,IERR)
            CALL MPI_WAIT(REQUEST,STATUS,IERR)
            IF (TIMERON)  CALL TIMER_STOP(T_RCOMM)
!DVM$       REGION 
!DVM$       PARALLEL (J), PRIVATE (J),TIE (R(J),W(J))
            DO  J = SEND_START,SEND_START + REDUCE_RECV_LENGTHS(I) - 1
               W(J) = W(J) + R(J)
            ENDDO  
!DVM$       END REGION
         ENDDO  

!--------------------------------------------------------------------- 
!  Exchange piece of q with transpose processor: 
!--------------------------------------------------------------------- 
         IF (L2NPCOLS .NE. 0) THEN
            IF (TIMERON)  CALL TIMER_START(T_RCOMM)

            CALL MPI_IRECV(R,EXCH_RECV_LENGTH,DP_TYPE,EXCH_PROC,1,MPI_CO
     &MM_WORLD,REQUEST,IERR)
!DVM$       ACTUAL (R(1:EXCH_RECV_LENGTH))
!DVM$       GET_ACTUAL (W(SEND_START:SEND_START+SEND_LEN))
            CALL MPI_SEND(W(SEND_START),SEND_LEN,DP_TYPE,EXCH_PROC,1,MPI
     &_COMM_WORLD,IERR)
            CALL MPI_WAIT(REQUEST,STATUS,IERR)
            IF (TIMERON)  CALL TIMER_STOP(T_RCOMM)
         ELSE  
!DVM$       REGION 
!DVM$       PARALLEL (J), PRIVATE (J),TIE (R(J),W(J))
            DO  J = 1,EXCH_RECV_LENGTH
               R(J) = W(J)
            ENDDO  
!DVM$       END REGION
         ENDIF  

!--------------------------------------------------------------------- 
!  At this point, r contains A.z 
!--------------------------------------------------------------------- 
         SUM = 0.0D0
!DVM$    REGION 
!DVM$    PARALLEL (J), PRIVATE (D,J),TIE (R(J),X(J)),REDUCTION (SUM (SUM
!DVM$&))
         DO  J = 1,LASTCOL - FIRSTCOL + 1
            D = X(J) - R(J)
            SUM = SUM + D * D
         ENDDO  
!DVM$    END REGION

!--------------------------------------------------------------------- 
!  Obtain d with a sum-reduce 
!--------------------------------------------------------------------- 
         DO  I = 1,L2NPCOLS
            IF (TIMERON)  CALL TIMER_START(T_RCOMM)
            CALL MPI_IRECV(D,1,DP_TYPE,REDUCE_EXCH_PROC(I),I,MPI_COMM_WO
     &RLD,REQUEST,IERR)
            CALL MPI_SEND(SUM,1,DP_TYPE,REDUCE_EXCH_PROC(I),I,MPI_COMM_W
     &ORLD,IERR)
            CALL MPI_WAIT(REQUEST,STATUS,IERR)
            IF (TIMERON)  CALL TIMER_STOP(T_RCOMM)
            SUM = SUM + D
         ENDDO  
         D = SUM
         IF (ME .EQ. ROOT)  RNORM = SQRT (D)
         IF (TIMERON)  CALL TIMER_STOP(T_CONJG)
         RETURN 

! end of routine conj_grad 
         END


!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         SUBROUTINE MAKEA (N, NZ, A, COLIDX, ROWSTR, NONZER, FIRSTROW, L
     &ASTROW, FIRSTCOL, LASTCOL, RCOND, AROW, ACOL, AELT, V, IV, SHIFT)

!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         IMPLICIT NONE
         INTEGER ::  N,NZ,NONZER
         INTEGER ::  FIRSTROW,LASTROW,FIRSTCOL,LASTCOL
         INTEGER ::  COLIDX(NZ),ROWSTR(N + 1)
         INTEGER ::  IV(2 * N + 1),AROW(NZ),ACOL(NZ)
         DOUBLE PRECISION ::  V(N + 1),AELT(NZ)
         DOUBLE PRECISION ::  RCOND,A(NZ),SHIFT

!--------------------------------------------------------------------- 
!       generate the test problem for benchmark 6 
!       makea generates a sparse matrix with a 
!       prescribed sparsity distribution 
!      
!       parameter    type        usage 
!      
!       input 
!      
!       n            i           number of cols/rows of matrix 
!       nz           i           nonzeros as declared array size 
!       rcond        r*8         condition number 
!       shift        r*8         main diagonal shift 
!      
!       output 
!      
!       a            r*8         array for nonzeros 
!       colidx       i           col indices 
!       rowstr       i           row pointers 
!      
!       workspace 
!      
!       iv, arow, acol i 
!       v, aelt        r*8 
!--------------------------------------------------------------------- 
         INTEGER ::  I,NNZA,IOUTER,IVELT,IVELT1,IROW,NZV,JCOL

!--------------------------------------------------------------------- 
!      nonzer is approximately  (int(sqrt(nnza /n))); 
!--------------------------------------------------------------------- 
         DOUBLE PRECISION ::  SIZE,RATIO,SCALE
         EXTERNAL SPARSE,SPRNVC,VECSET
         INTENT(INOUT) IV,V,AELT,ACOL,AROW,ROWSTR,COLIDX,A
         INTENT(IN) SHIFT,RCOND,LASTCOL,FIRSTCOL,LASTROW,FIRSTROW,NONZER
     &,NZ,N
         SIZE = 1.0D0
         RATIO = RCOND** (1.0D0 / DFLOAT (N))
         NNZA = 0

!--------------------------------------------------------------------- 
!  Initialize iv(n+1 .. 2n) to zero. 
!  Used by sprnvc to mark nonzero positions 
!--------------------------------------------------------------------- 
!DVM$    GET_ACTUAL (IV)
         DO  I = 1,N
            IV(N + I) = 0
         ENDDO  
!DVM$    ACTUAL (IV)
         DO  IOUTER = 1,N
            NZV = NONZER
!DVM$       GET_ACTUAL (IV)
            CALL SPRNVC(N,NZV,V,COLIDX,IV(1),IV(N + 1))
!DVM$       ACTUAL (IV)
            CALL VECSET(N,V,COLIDX,NZV,IOUTER,.5D0)
!DVM$       GET_ACTUAL (ACOL,AELT,AROW,COLIDX,V)
            DO  IVELT = 1,NZV
               JCOL = COLIDX(IVELT)
               IF (JCOL .GE. FIRSTCOL .AND. JCOL .LE. LASTCOL) THEN
                  SCALE = SIZE * V(IVELT)
                  DO  IVELT1 = 1,NZV
                     IROW = COLIDX(IVELT1)
                     IF (IROW .GE. FIRSTROW .AND. IROW .LE. LASTROW) THE
     &N
                        NNZA = NNZA + 1
                        IF (NNZA .GT. NZ)  GOTO 9999
                        ACOL(NNZA) = JCOL
                        AROW(NNZA) = IROW
                        AELT(NNZA) = V(IVELT1) * SCALE
                     ENDIF  
                  ENDDO  
               ENDIF  
            ENDDO  
!DVM$       ACTUAL (ACOL,AELT,AROW)
            SIZE = SIZE * RATIO
         ENDDO  

!--------------------------------------------------------------------- 
!       ... add the identity * rcond to the generated matrix to bound 
!           the smallest eigenvalue from below by rcond 
!--------------------------------------------------------------------- 
!DVM$    GET_ACTUAL (ACOL,AELT,AROW)
         DO  I = FIRSTROW,LASTROW
            IF (I .GE. FIRSTCOL .AND. I .LE. LASTCOL) THEN
               IOUTER = N + I
               NNZA = NNZA + 1
               IF (NNZA .GT. NZ)  GOTO 9999
               ACOL(NNZA) = I
               AROW(NNZA) = I
               AELT(NNZA) = RCOND - SHIFT
            ENDIF  
         ENDDO  
!DVM$    ACTUAL (ACOL,AELT,AROW)

!--------------------------------------------------------------------- 
!       ... make the sparse matrix from list of elements with duplicates 
!           (v and iv are used as  workspace) 
!--------------------------------------------------------------------- 
!DVM$    GET_ACTUAL (IV)
         CALL SPARSE(A,COLIDX,ROWSTR,N,AROW,ACOL,AELT,FIRSTROW,LASTROW,V
     &,IV(1),IV(N + 1),NNZA)
!DVM$    ACTUAL (IV)
         RETURN 
9999     CONTINUE
         WRITE (UNIT = *,FMT = *) 'Space for matrix elements exceeded in
     & makea'
         WRITE (UNIT = *,FMT = *) 'nnza, nzmax = ',NNZA,NZ
         WRITE (UNIT = *,FMT = *) ' iouter = ',IOUTER
         STOP
         END


!-------end   of makea------------------------------ 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         SUBROUTINE SPARSE (A, COLIDX, ROWSTR, N, AROW, ACOL, AELT, FIRS
     &TROW, LASTROW, X, MARK, NZLOC, NNZA)

!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         IMPLICIT NONE

      INCLUDE 'npbparams.h' 
         INTEGER ::  COLIDX(*),ROWSTR(NA + 1)
         INTEGER ::  FIRSTROW,LASTROW
         INTEGER ::  N,AROW(*),ACOL(*),NNZA
         DOUBLE PRECISION ::  A(*),AELT(*)

!--------------------------------------------------------------------- 
!       rows range from firstrow to lastrow 
!       the rowstr pointers are defined for nrows = lastrow-firstrow+1 values 
!--------------------------------------------------------------------- 
         INTEGER ::  NZLOC(N),NROWS
         DOUBLE PRECISION ::  X(N)
         LOGICAL ::  MARK(N)

!--------------------------------------------------- 
!       generate a sparse matrix from a list of 
!       [col, row, element] tri 
!--------------------------------------------------- 
         INTEGER ::  I,J,JAJP1,NZA,K,NZROW
         DOUBLE PRECISION ::  XI
         INTENT(INOUT) NZLOC,MARK,X,ROWSTR,COLIDX,A
         INTENT(IN) NNZA,LASTROW,FIRSTROW,AELT,ACOL,AROW,N

!--------------------------------------------------------------------- 
!    how many rows of result 
!--------------------------------------------------------------------- 
         NROWS = LASTROW - FIRSTROW + 1

!--------------------------------------------------------------------- 
!     ...count the number of triples in each row 
!--------------------------------------------------------------------- 
!DVM$    REGION 
!DVM$    PARALLEL (J), PRIVATE (J),TIE (ROWSTR(J))
         DO  J = 1,N
            ROWSTR(J) = 0
            MARK(J) = .FALSE.
         ENDDO  
!DVM$    END REGION
         ROWSTR(N + 1) = 0
!DVM$    ACTUAL (ROWSTR(N + 1))
!DVM$    GET_ACTUAL (AROW,ROWSTR)
         DO  NZA = 1,NNZA
            J = AROW(NZA) - FIRSTROW + 1 + 1
            ROWSTR(J) = ROWSTR(J) + 1
         ENDDO  
!DVM$    ACTUAL (ROWSTR)
         ROWSTR(1) = 1
!DVM$    ACTUAL (ROWSTR(1))
!DVM$    REGION 
!DVM$    PARALLEL (J), PRIVATE (J),TIE (ROWSTR(J)),ACROSS (ROWSTR(1:0))
         DO  J = 2,NROWS + 1
            ROWSTR(J) = ROWSTR(J) + ROWSTR(J - 1)
         ENDDO  
!DVM$    END REGION

!--------------------------------------------------------------------- 
!     ... rowstr(j) now is the location of the first nonzero 
!           of row j of a 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!     ... do a bucket sort of the triples on the row index 
!--------------------------------------------------------------------- 
!DVM$    GET_ACTUAL (A,ACOL,AELT,AROW,COLIDX,ROWSTR)
         DO  NZA = 1,NNZA
            J = AROW(NZA) - FIRSTROW + 1
            K = ROWSTR(J)
            A(K) = AELT(NZA)
            COLIDX(K) = ACOL(NZA)
            ROWSTR(J) = ROWSTR(J) + 1
         ENDDO  
!DVM$    ACTUAL (A,COLIDX,ROWSTR)

!--------------------------------------------------------------------- 
!       ... rowstr(j) now points to the first element of row j+1 
!--------------------------------------------------------------------- 
!DVM$    REGION 
!DVM$    PARALLEL (J), PRIVATE (J),TIE (ROWSTR(J)),ACROSS (ROWSTR(1:0))
         DO  J = NROWS,1,(-(1))
            ROWSTR(J + 1) = ROWSTR(J)
         ENDDO  
!DVM$    END REGION
         ROWSTR(1) = 1
!DVM$    ACTUAL (ROWSTR(1))

!--------------------------------------------------------------------- 
!       ... generate the actual output rows by adding elements 
!--------------------------------------------------------------------- 
         NZA = 0
!DVM$    GET_ACTUAL (MARK,X)
         DO  I = 1,N
            X(I) = 0.0
            MARK(I) = .FALSE.
         ENDDO  
!DVM$    ACTUAL (MARK,X)
!DVM$    GET_ACTUAL (ROWSTR)
         JAJP1 = ROWSTR(1)
!DVM$    GET_ACTUAL (A,COLIDX,MARK,NZLOC,ROWSTR,X)
         DO  J = 1,NROWS
            NZROW = 0

!--------------------------------------------------------------------- 
!          ...loop over the jth row of a 
!--------------------------------------------------------------------- 
            DO  K = JAJP1,ROWSTR(J + 1) - 1
               I = COLIDX(K)
               X(I) = X(I) + A(K)
               IF (.NOT.(MARK(I)) .AND. X(I) .NE. 0.D0) THEN
                  MARK(I) = .TRUE.
                  NZROW = NZROW + 1
                  NZLOC(NZROW) = I
               ENDIF  
            ENDDO  

!--------------------------------------------------------------------- 
!          ... extract the nonzeros of this row 
!--------------------------------------------------------------------- 
            DO  K = 1,NZROW
               I = NZLOC(K)
               MARK(I) = .FALSE.
               XI = X(I)
               X(I) = 0.D0
               IF (XI .NE. 0.D0) THEN
                  NZA = NZA + 1
                  A(NZA) = XI
                  COLIDX(NZA) = I
               ENDIF  
            ENDDO  
            JAJP1 = ROWSTR(J + 1)
            ROWSTR(J + 1) = NZA + ROWSTR(1)
         ENDDO  
!DVM$    ACTUAL (A,COLIDX,MARK,NZLOC,ROWSTR,X)

!C       write (*, 11000) nza 
         RETURN 
11000      FORMAT ( //,'final nonzero count in sparse ',            /,'n
     &umber of nonzeros       = ', I16 ) 
         END


!-------end   of sparse----------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         SUBROUTINE SPRNVC (N, NZ, V, IV, NZLOC, MARK)

!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         IMPLICIT NONE
         DOUBLE PRECISION ::  V(*)
         INTEGER ::  N,NZ,IV(*),NZLOC(N),NN1
         INTEGER ::  MARK(N)
         COMMON /URANDO/AMULT,TRAN
         DOUBLE PRECISION ::  AMULT,TRAN

!--------------------------------------------------------------------- 
!       generate a sparse n-vector (v, iv) 
!       having nzv nonzeros 
!      
!       mark(i) is set to 1 if position i is nonzero. 
!       mark is all zero on entry and is reset to all zero before exit 
!       this corrects a performance bug found by John G. Lewis, caused by 
!       reinitialization of mark on every one of the n calls to sprnvc 
!--------------------------------------------------------------------- 
         INTEGER ::  NZROW,NZV,II,I,ICNVRT
         EXTERNAL RANDLC,ICNVRT
         DOUBLE PRECISION ::  RANDLC,VECELT,VECLOC
         INTENT(INOUT) MARK,NZLOC
         INTENT(IN) NZ,N
         INTENT(OUT) IV,V
         NZV = 0
         NZROW = 0
         NN1 = 1
50       CONTINUE
         NN1 = 2 * NN1
         IF (NN1 .LT. N)  GOTO 50

!--------------------------------------------------------------------- 
!    nn1 is the smallest power of two not less than n 
!--------------------------------------------------------------------- 
100      CONTINUE
         IF (NZV .GE. NZ)  GOTO 110
         VECELT = RANDLC (TRAN,AMULT)

!--------------------------------------------------------------------- 
!   generate an integer between 1 and n in a portable manner 
!--------------------------------------------------------------------- 
         VECLOC = RANDLC (TRAN,AMULT)
         I = ICNVRT (VECLOC,NN1) + 1
         IF (I .GT. N)  GOTO 100

!--------------------------------------------------------------------- 
!  was this integer generated already? 
!--------------------------------------------------------------------- 
         IF (MARK(I) .EQ. 0) THEN
            MARK(I) = 1
!DVM$       ACTUAL (MARK(I))
            NZROW = NZROW + 1
            NZLOC(NZROW) = I
!DVM$       ACTUAL (NZLOC(NZROW))
            NZV = NZV + 1
            V(NZV) = VECELT
            IV(NZV) = I
         ENDIF  
         GOTO 100
110      CONTINUE
!DVM$    GET_ACTUAL (MARK,NZLOC)
         DO  II = 1,NZROW
            I = NZLOC(II)
            MARK(I) = 0
         ENDDO  
!DVM$    ACTUAL (MARK)
         RETURN 
         END


!-------end   of sprnvc----------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         FUNCTION ICNVRT (X, IPWR2) 

!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         IMPLICIT NONE
         DOUBLE PRECISION ::  X
         INTEGER ::  IPWR2,ICNVRT
         INTENT(IN) IPWR2,X

!--------------------------------------------------------------------- 
!    scale a double precision number x in (0,1) by a power of 2 and chop it 
!--------------------------------------------------------------------- 
         ICNVRT = INT (IPWR2 * X)
         RETURN 
         END


!-------end   of icnvrt----------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         SUBROUTINE VECSET (N, V, IV, NZV, I, VAL)

!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
         IMPLICIT NONE
         INTEGER ::  N,IV(*),NZV,I,K
         DOUBLE PRECISION ::  V(*),VAL

!--------------------------------------------------------------------- 
!       set ith element of sparse vector (v, iv) with 
!       nzv nonzeros to val 
!--------------------------------------------------------------------- 
         LOGICAL ::  SET
         INTENT(INOUT) NZV,IV
         INTENT(IN) VAL,I
         INTENT(OUT) V
         SET = .FALSE.
!DVM$    GET_ACTUAL (IV,V)
         DO  K = 1,NZV
            IF (IV(K) .EQ. I) THEN
               V(K) = VAL
               SET = .TRUE.
            ENDIF  
         ENDDO  
!DVM$    ACTUAL (V)
         IF (.NOT.(SET)) THEN
            NZV = NZV + 1
            V(NZV) = VAL
            IV(NZV) = I
         ENDIF  
         RETURN 

!-------end   of vecset----------------------------- 
         END

