pdhseqrdriver.f

Go to the documentation of this file.

***********************************************************************
*     Test program for ScaLAPACK-style routine PDHSEQR                *
***********************************************************************
*
*     Contributor: Robert Granat and Meiyue Shao
*     This version is of Feb 2011.
*
      PROGRAM pdhseqrdriver
*
*     Declarations
*
      IMPLICIT NONE
*     ...Parameters...
      LOGICAL           balance, comphess, compresi,
     $                  comporth
      LOGICAL           DEBUG, prn, timesteps, barr,
     $                  uni_lapack
      INTEGER           slv_min, slv_max
      parameter( debug = .false.,
     $                    prn = .false.,
     $                    timesteps = .true.,
     $                    comphess = .true.,
     $                    compresi = .true.,
     $                    comporth = .true.,
     $                    balance = .true.,
     $                    barr = .false.,
*     Solver: 1-PDLAQR1, 2-PDHSEQR.
     $                    slv_min = 2, slv_max = 2,
     $                    uni_lapack = .true. )
      INTEGER           n, nb, arsrc, acsrc
      parameter(
*     Problem size.
     $                    n = 500, nb = 50,
*     What processor should hold the first element in A?
     $                    arsrc = 0, acsrc = 0 )
      INTEGER           block_cyclic_2d, csrc_, ctxt_, dlen_, dt_,
     $                  lld_, mb_, m_, nb_, n_, rsrc_
      parameter( block_cyclic_2d = 1, dlen_ = 9, dt_ = 1,
     $                    ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
     $                    rsrc_ = 7, csrc_ = 8, lld_ = 9 )
      INTEGER           dpalloc, intallc
      INTEGER           dpsiz, intsz, nout, izero
      parameter( dpsiz = 8, dpalloc = 8 000 000,
     $                    intsz = 4, intallc = 8 000 000,
     $                    nout = 6, izero = 0)
      DOUBLE PRECISION  zero, one, two
      parameter( zero = 0.0d+00, one = 1.0d+00, two = 2.0d+00 )
*
*     ...Local Scalars...
      INTEGER           ictxt, iam, nprow, npcol, myrow, mycol,
     $                  sys_nprocs, nprocs, arows, acols, temp_ictxt
      INTEGER           threads
      INTEGER           info, ktop, kbot, ilo, ihi, i
      INTEGER           ipa, ipacpy, ipq, wr1, wi1, wr2, wi2, ipw1,
     $                  ipw2, ipiw
      INTEGER           totit, sweeps, totns, hess
      DOUBLE PRECISION  eps, thresh
      DOUBLE PRECISION  stamp, tottime, t_ba, t_gen, t_hs, t_sch, t_qr,
     $                  t_res, itpereig, swpspeig, nspeig, speedup, 
     $                  efficiency
      DOUBLE PRECISION  rnorm, anorm, r1, orth, o1, o2, dpdum, elem1,
     $                  elem2, elem3, ediff
      INTEGER           solver
      CHARACTER*6       passed
*
*     ...Local Arrays...
      INTEGER           desca( dlen_ ), descq( dlen_ ), descvec( dlen_ )
      DOUBLE PRECISION  scale( n )
      DOUBLE PRECISION, ALLOCATABLE :: mem(:)
      INTEGER, ALLOCATABLE :: imem(:)
*
*     ...Intrinsic Functions...
      INTRINSIC         int, dble, sqrt, max, min
*
*     ...External Functions...
      INTEGER           numroc
      DOUBLE PRECISION  pdlamch, pdlange, mpi_wtime
      EXTERNAL          blacs_pinfo, blacs_get, blacs_gridinit,
     $                  blacs_gridinfo, blacs_gridexit, blacs_exit
      EXTERNAL          numroc, pdlamch, pdlaset, pdgehrd, pdlange
      EXTERNAL          dgebal, dgehrd
      EXTERNAL          mpi_wtime
      EXTERNAL          pdgebal
      EXTERNAL          pdmatgen2
*
*     ...Executable statements...
*
      CALL blacs_pinfo( iam, sys_nprocs )
      nprow = int( sqrt( dble(sys_nprocs) ) )
      npcol = sys_nprocs / nprow
      CALL blacs_get( 0, 0, ictxt )
      CALL blacs_gridinit( ictxt, '2D', nprow, npcol )
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
c      print*, iam, ictxt, myrow, mycol
c      IF ( ( MYROW.GE.NPROW ) .OR. ( MYCOL.GE.NPCOL ) ) GO TO 777
      IF ( ictxt.LT.0 ) GO TO 777
*
*     Read out the number of underlying threads and set stack size in
*     kilobytes.
*
    thresh = 30.0
      tottime = mpi_wtime()
      t_gen = 0.0d+00
      t_res = 0.0d+00
      t_sch = 0.0d+00
*
*     Allocate and Init memory with zeros.
*
      info = 0
      ALLOCATE ( mem( dpalloc ), stat = info )
      IF( info.NE.0 ) THEN
         WRITE(*,*) '% Could not allocate MEM. INFO = ', info
         GO TO 777
      END IF
      ALLOCATE ( imem( intallc ), stat = info )
      IF( info.NE.0 ) THEN
         WRITE(*,*) '% Could not allocate IMEM. INFO = ', info
         GO TO 777
      END IF
      mem( 1:dpalloc ) = zero
      imem( 1:intallc ) = izero
*
*     Get machine epsilon.
*
      eps = pdlamch( ictxt, 'Epsilon' )      
*
*     Print welcoming message.
*
      IF( iam.EQ.0 ) THEN
         WRITE(*,*)
         WRITE(*,*) 'ScaLAPACK Test for PDHSEQR'
         WRITE(*,*) 
         WRITE(*,*) 'epsilon   = ', eps
         WRITE(*,*) 'threshold = ', thresh
         WRITE(*,*)
         WRITE(*,*) 'Residual and Orthogonality Residual computed by:'
         WRITE(*,*)
         WRITE(*,*) 'Residual      = ',
     $   ' || T - Q^T*A*Q ||_F / ( ||A||_F * eps * sqrt(N) )'
           WRITE(*,*)
         WRITE(*,*) 'Orthogonality = ',
     $   ' MAX( || I - Q^T*Q ||_F, || I - Q*Q^T ||_F ) / ',
     $   ' (eps * N)'
           WRITE(*,*) 
           WRITE(*,*) 
     $  'Test passes if both residuals are less then threshold'        
       WRITE( nout, * )
       WRITE( nout, fmt = 9995 )
       WRITE( nout, fmt = 9994 )
      END IF
*
*     Loop over problem parameters.
*
      DO ktop = 1, 1
      DO kbot = n, n
      DO solver = slv_max, slv_min, -1
*
*        Set INFO to zero for this run.
*
         info = 0
         nprocs = nprow*npcol
         temp_ictxt = ictxt
*
*        Count the number of rows and columns of current problem
*        for the current block sizes and grid properties.
*
         stamp = mpi_wtime()
         arows = numroc( n, nb, myrow, 0, nprow )
         acols = numroc( n, nb, mycol, 0, npcol )
*
*        Set up matrix descriptors.
*
         IF( debug ) WRITE(*,*) '% #', iam, ': Set up descriptors...'
         IF( barr ) CALL blacs_barrier(ictxt, 'A')
         CALL descinit( desca, n, n, nb, nb, min(arsrc,nprow-1),
     $        min(npcol-1,acsrc), temp_ictxt, max(1, arows), info )
         IF ( info.NE.0 ) THEN
            WRITE(*,*) "% DESCINIT DESCA failed, INFO =", info
            GO TO 999
         END IF
         CALL descinit( descq, n, n, nb, nb, min(arsrc,nprow-1),
     $        min(npcol-1,acsrc), temp_ictxt, max(1, arows), info )
         IF ( info.NE.0 ) THEN
            WRITE(*,*) "% DESCINIT DESCQ failed, INFO =", info
            GO TO 999
         END IF
         CALL descinit( descvec, n, 1, n, 1, min(arsrc,nprow-1),
     $        min(npcol-1,acsrc), temp_ictxt, n, info )
         IF ( info.NE.0 ) THEN
            WRITE(*,*) "% DESCINIT DESCVEC failed, INFO =", info
            GO TO 999
         END IF
*
*        Assign pointer for ScaLAPACK arrays - first set DP memory.
*
         IF( debug ) WRITE(*,*) '% #', iam, ': Assign pointers...'
         ipa    = 1
         ipacpy = ipa + desca( lld_ ) * acols
         ipq    = ipacpy + desca( lld_ ) * acols
         wr1    = ipq + descq( lld_ ) * acols
         wi1    = wr1 + n
         wr2    = wi1 + n
         wi2    = wr2 + n
         ipw1   = wi2 + n
         ipw2   = ipw1 + desca( lld_ ) * acols
         IF( debug ) WRITE(*,*) '% (IPW2,DPALLOC):', ipw2, dpalloc
*         PRINT*, '%', IPA, IPACPY, IPQ, WR1, WI1, WR2, WI2, IPW1, IPW2
         IF( ipw2+desca(lld_)*acols .GT. dpalloc+1 ) THEN
            WRITE(*,*) '% Not enough DP memory!'
            GO TO 999
         END IF
*
*        Then set integer memory pointers.
*
         ipiw = 1
*
*        Generate testproblem.
*
         IF( barr ) CALL blacs_barrier(ictxt, 'A')
         CALL pdlaset( 'All over', n, n, zero, one, mem(ipq), 1, 1,
     $        descq )
         CALL pdmatgen2( temp_ictxt, 'Random', 'NoDiagDominant',
     $        n, n, nb, nb, mem(ipa), desca( lld_ ), 0, 0, 7, 0,
     $        arows, 0, acols, myrow, mycol, nprow, npcol )
         IF( .NOT. comphess ) THEN
            CALL pdlaset( 'Lower triangular', n-2, n-2, zero, zero,
     $           mem(ipa), 3, 1, desca )
            CALL pdlaset( 'All over', n, n, zero, one, mem(ipq),
     $           1, 1, descq )
            IF( ktop.GT.1 )
     $           CALL pdlaset( 'Lower triangular', ktop-1, ktop-1,
     $           zero, zero, mem(ipa), 2, 1, descq )
            IF( kbot.LT.n )
     $           CALL pdlaset( 'Lower triangular', n-kbot, n-kbot,
     $           zero, zero, mem(ipa), kbot+1, kbot, descq )
         END IF
*
*        Do balancing if general matrix.
*
         IF( barr ) CALL blacs_barrier(ictxt, 'A')
         t_ba = mpi_wtime()
         IF( comphess .AND. balance ) THEN
            IF( nprocs.EQ.1 .AND. solver.NE.2 .AND. uni_lapack ) THEN
               IF( debug ) WRITE(*,*) '% #', iam, ': == dgebal =='
               CALL dgebal( 'Both', n, mem(ipa), desca(lld_), ilo,
     $              ihi, scale, info )
               IF ( info.NE.0 ) THEN
                  WRITE(*,*) "% DGEBAL failed, INFO =", info
                  GO TO 999
               END IF
            ELSE
               IF( debug ) WRITE(*,*) '% #', iam, ': == pdgebal =='
               CALL pdgebal( 'Both', n, mem(ipa), desca, ilo, ihi,
     $              scale, info )
               IF ( info.NE.0 ) THEN
                  WRITE(*,*) "% PDGEBAL failed, INFO =", info
                  GO TO 999
               END IF
            END IF
         ELSEIF( comphess ) THEN
            ilo = 1
            ihi = n
         ELSE
            ilo = ktop
            ihi = kbot
         END IF
         t_ba = mpi_wtime() - t_ba
c         IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
c     $      ' %%% Balancing took in seconds:',T_BA
         IF( debug ) WRITE(*,*) '% #', iam, ': ILO,IHI=',ilo,ihi
*
*        Make a copy of A.
*
         IF( barr ) CALL blacs_barrier(ictxt, 'A')
         IF( debug ) WRITE(*,*) '% #', iam, ': Copy matrix A'
         CALL pdlacpy( 'All', n, n, mem(ipa), 1, 1, desca, mem(ipacpy),
     $                 1, 1, desca )
*
*        Print matrices to screen in debugging mode.
*
         IF( prn )
     $      CALL pdlaprnt( n, n, mem(ipacpy), 1, 1, desca, 0, 0,
     $           'A', nout, mem(ipw1) )
         t_gen = t_gen + mpi_wtime() - stamp - t_ba
c         IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
c     $      ' %%% Generation took in seconds:',T_GEN
*
*        Only compute the Hessenberg form if necessary.
*
         t_hs = mpi_wtime()
         IF( .NOT. comphess ) GO TO 30
*
*        Reduce A to Hessenberg form.
*
         IF( barr ) CALL blacs_barrier(ictxt, 'A')
         IF( debug ) WRITE(*,*) '% #', iam,
     $      ': Reduce to Hessenberg form...N=',n, ilo,ihi
*         PRINT*, '% PDGEHRD: IPW2,MEM(IPW2)', IPW2, MEM(IPW2)
         IF( nprocs.EQ.1 .AND. solver.NE.2 .AND. uni_lapack ) THEN
            IF( debug ) WRITE(*,*) '% #', iam, ': == dgehrd =='
            CALL dgehrd( n, ilo, ihi, mem(ipa), desca(lld_),
     $           mem(ipw1), mem(ipw2), -1, info )
            IF (dpalloc-ipw2.LT.mem(ipw2)) THEN
               WRITE(*,*) "% Not enough memory for DGEHRD"
               GO TO 999
            END IF
            CALL dgehrd( n, ilo, ihi, mem(ipa), desca(lld_),
     $           mem(ipw1), mem(ipw2), dpalloc-ipw2, info )
            IF ( info.NE.0 ) THEN
               WRITE(*,*) "% DGEHRD failed, INFO =", info
               GO TO 999
            END IF
         ELSE
            IF( debug ) WRITE(*,*) '% #', iam, ': == pdgehrd =='
            CALL pdgehrd( n, ilo, ihi, mem(ipa), 1, 1, desca, mem(ipw1),
     $           mem(ipw2), -1, info )
            IF (dpalloc-ipw2.LT.mem(ipw2)) THEN
               WRITE(*,*) "% Not enough memory for PDGEHRD"
               GO TO 999
            END IF
            CALL pdgehrd( n, ilo, ihi, mem(ipa), 1, 1, desca, mem(ipw1),
     $           mem(ipw2), dpalloc-ipw2, info )
            IF ( info.NE.0 ) THEN
               WRITE(*,*) "% PDGEHRD failed, INFO =", info
               GO TO 999
            END IF
         END IF
*
*        Form Q explicitly.
*
         IF( barr ) CALL blacs_barrier(ictxt, 'A')
         IF( debug ) WRITE(*,*) '% #', iam, ':Form Q explicitly'
*         PRINT*, '% PDORMHR: IPW2,MEM(IPW2)', IPW2, MEM(IPW2)
         IF( debug ) WRITE(*,*) '% #', iam, ': == pdormhr =='
         CALL pdormhr( 'L', 'N', n, n, ilo, ihi, mem(ipa), 1, 1,
     $        desca, mem(ipw1), mem(ipq), 1, 1, descq, mem(ipw2),
     $        -1, info )
         IF (dpalloc-ipw2.LT.mem(ipw2)) THEN
            WRITE(*,*) "% Not enough memory for PDORMHR"
            GO TO 999
         END IF
         CALL pdormhr( 'L', 'N', n, n, ilo, ihi, mem(ipa), 1, 1,
     $        desca, mem(ipw1), mem(ipq), 1, 1, descq, mem(ipw2),
     $        dpalloc-ipw2, info )
         IF ( info.NE.0 ) THEN
            WRITE(*,*) "% PDORMHR failed, INFO =", info
            GO TO 999
         END IF
*
*        Extract the upper Hessenberg part of A.
*
         CALL pdlaset( 'Lower triangular', n-2, n-2, zero, zero,
     $        mem(ipa), 3, 1, desca )
*
*        Print reduced matrix A in debugging mode.
*
         IF( prn ) THEN
            CALL pdlaprnt( n, n, mem(ipa), 1, 1, desca, 0, 0, 'H', nout,
     $           mem(ipw1) )
            CALL pdlaprnt( n, n, mem(ipq), 1, 1, descq, 0, 0, 'Q', nout,
     $           mem(ipw1) )
         END IF
*
 30      CONTINUE
         t_hs = mpi_wtime() - t_hs
c         IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
c     $      ' %%% Hessenberg took in seconds:',T_HS
*
*        Compute the real Schur form of the Hessenberg matrix A.
*
         IF( barr ) CALL blacs_barrier(ictxt, 'A')
         t_qr = mpi_wtime()
         IF( solver.EQ.1 ) THEN
            IF( debug ) WRITE(*,*) '% #', iam, ': == pdlaqr1 =='
*            PRINT*, '% PDLAQR1: IPW1,MEM(IPW1)', IPW1, MEM(IPW1)
            CALL pdlaqr1( .true., .true., n, ilo, ihi, mem(ipa), desca,
     $           mem(wr1), mem(wi1), ilo, ihi, mem(ipq), descq,
     $           mem(ipw1), -1, imem, -1, info )
            IF (dpalloc-ipw1.LT.mem(ipw1)) THEN
               WRITE(*,*) "% Not enough DP memory for PDLAQR1"
               GO TO 999
            END IF
            IF (intallc.LT.imem(1)) THEN
               WRITE(*,*) "% Not enough INT memory for PDLAQR1"
               GO TO 999
            END IF
            CALL pdlaqr1( .true., .true., n, ilo, ihi, mem(ipa), desca,
     $           mem(wr1), mem(wi1), ilo, ihi, mem(ipq), descq,
     $           mem(ipw1), dpalloc-ipw1+1, imem, intallc, info )
            IF (info.NE.0) THEN
               WRITE(*,*) "% PDLAQR1: INFO =", info
            END IF
         ELSEIF( solver.EQ.2 ) THEN
            IF( debug ) WRITE(*,*) '% #', iam, ': == pdhseqr =='
*            PRINT*, '% PDHSEQR: IPW1,MEM(IPW1)', IPW1, MEM(IPW1)
         IF( barr ) CALL blacs_barrier(ictxt, 'A')
            CALL pdhseqr( 'Schur', 'Vectors', n, ilo, ihi, mem(ipa),
     $           desca, mem(wr2), mem(wi2), mem(ipq), descq, mem(ipw1),
     $           -1, imem, -1, info )
         IF( barr ) CALL blacs_barrier(ictxt, 'A')
            IF (dpalloc-ipw1.LT.mem(ipw1)) THEN
               WRITE(*,*) "% Not enough dp memory for pdhseqr"
               GO TO 999
            END IF
.LT.            IF (INTALLCIMEM(1)) THEN
               WRITE(*,*) "% Not enough int memory for pdhseqr"
               GO TO 999
            END IF
         IF( BARR ) CALL BLACS_BARRIER(ICTXT, 'A')
            CALL PDHSEQR( 'Schur', 'Vectors', N, ILO, IHI, MEM(IPA),
     $           DESCA, MEM(WR2), MEM(WI2), MEM(IPQ), DESCQ, MEM(IPW1),
     $           DPALLOC-IPW1+1, IMEM, INTALLC, INFO )
.NE.            IF (INFO0) THEN
               WRITE(*,*) "% PDHSEQR: info =", INFO
            END IF
         ELSE
             WRITE(*,*) '% ERROR: Illegal SOLVER number!'
             GO TO 999
         END IF
         T_QR = MPI_WTIME() - T_QR
c         IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
c     $      ' %%% QR-algorithm took in seconds:',T_QR
         T_SCH = T_SCH + T_QR + T_HS + T_BA
*         TOTIT = IMEM(1)
*         SWEEPS = IMEM(2)
*         TOTNS = IMEM(3)
         ITPEREIG = DBLE(TOTIT) / DBLE(N)
         SWPSPEIG = DBLE(SWEEPS) / DBLE(N)
         NSPEIG = DBLE(TOTNS) / DBLE(N)
*
*        Print reduced matrix A in debugging mode.
*
         IF( PRN ) THEN
            CALL PDLAPRNT( N, N, MEM(IPA), 1, 1, DESCA, 0, 0, 'T',
     $           NOUT, MEM(IPW1) )
            CALL PDLAPRNT( N, N, MEM(IPQ), 1, 1, DESCQ, 0, 0, 'Z',
     $           NOUT, MEM(IPW1) )
         END IF
*
*        Check that returned Schur form is really a quasi-triangular
*        matrix.
*
         HESS = 0
         DO I = 1, N-1
.GT.            IF( I1 ) THEN
               CALL PDELGET( 'All', '1-Tree', ELEM1, MEM(IPA), I, I-1,
     $              DESCA )
            ELSE
               ELEM1 = ZERO
            END IF
            CALL PDELGET( 'All', '1-Tree', ELEM2, MEM(IPA), I+1, I,
     $           DESCA )
.LT.            IF( IN-1 ) THEN
               CALL PDELGET( 'All', '1-Tree', ELEM3, MEM(IPA), I+2, I+1,
     $              DESCA )
            ELSE
               ELEM3 = ZERO
            END IF
.NE..AND..GT.            IF( ELEM2ZERO  ABS(ELEM1)+ABS(ELEM2)+ABS(ELEM3)
     $         ABS(ELEM2) ) HESS = HESS + 1
         END DO
*
*        Compute residual norms and other results:
*
*           1) RNORM = || T - Q'*A*Q ||_F / ||A||_F
*           2) ORTH  = MAX( || I - Q'*Q ||_F, || I - Q*Q' ||_F ) /
*                  (epsilon*N)
*
         STAMP = MPI_WTIME()
         IF( COMPRESI ) THEN
            IF( DEBUG ) WRITE(*,*) '% #', IAM, ': Compute residuals 1'
            IF( DEBUG ) WRITE(*,*) '% #', IAM, ': pdgemm 3'
            CALL PDGEMM( 'N', 'N', N, N, N, ONE, MEM(IPACPY), 1, 1,
     $           DESCA, MEM(IPQ), 1, 1, DESCQ, ZERO, MEM(IPW1), 1, 1,
     $           DESCA )
            IF( DEBUG ) WRITE(*,*) '% #', IAM, ': pdgemm 4'
            IF( DEBUG ) WRITE(*,*) '% #', IAM, ': N=',N
            IF( DEBUG ) WRITE(*,*) '% #', IAM, ': DESCA=',DESCA(1:DLEN_)
            IF( DEBUG ) WRITE(*,*) '% #', IAM, ': DESCQ=',DESCQ(1:DLEN_)
            CALL PDGEMM( 'T', 'N', N, N, N, -ONE, MEM(IPQ), 1, 1,
     $           DESCQ, MEM(IPW1), 1, 1, DESCA, ONE, MEM(IPA), 1, 1,
     $           DESCA )
            R1 = PDLANGE( 'Frobenius', N, N, MEM(IPA), 1, 1, DESCA,
     $           DPDUM )
            ANORM = PDLANGE( 'Frobenius', N, N, MEM(IPACPY), 1, 1,
     $           DESCA, DPDUM )
.GT.            IF( ANORMZERO )THEN
               RNORM = R1 / (ANORM*EPS*SQRT(DBLE(N)))
            ELSE
               RNORM = R1
            END IF
         ELSE
            RNORM = 0.0D0
         END IF
*
         IF( COMPORTH ) THEN
            IF( DEBUG ) WRITE(*,*) '% #', IAM, ': Compute residuals 2'
            CALL PDLASET( 'All', N, N, ZERO, ONE, MEM(IPW1), 1, 1,
     $           DESCQ )
            CALL PDLACPY( 'All', N, N, MEM(IPQ), 1, 1, DESCQ, MEM(IPW2),
     $           1, 1, DESCQ )
            CALL PDGEMM( 'T', 'N', N, N, N, -ONE, MEM(IPQ), 1, 1, DESCQ,
     $           MEM(IPW2), 1, 1, DESCQ, ONE, MEM(IPW1), 1, 1, DESCQ )
            O1 = PDLANGE( 'Frobenius', N, N, MEM(IPW1), 1, 1, DESCQ,
     $           DPDUM )
            CALL PDLASET( 'All', N, N, ZERO, ONE, MEM(IPW1), 1, 1,
     $           DESCQ )
            CALL PDGEMM( 'N', 'T', N, N, N, -ONE, MEM(IPQ), 1, 1, DESCQ,
     $           MEM(IPW2), 1, 1, DESCQ, ONE, MEM(IPW1), 1, 1, DESCQ )
            O2 = PDLANGE( 'Frobenius', N, N, MEM(IPW1), 1, 1, DESCQ,
     $           DPDUM )
            ORTH = MAX(O1,O2) / (EPS*DBLE(N))
         ELSE
            ORTH = 0.0D0
         END IF
*
         T_RES = T_RES + MPI_WTIME() - STAMP
c         IF( TIMESTEPS.AND.IAM.EQ.0 ) WRITE(*,*)
c     $      ' %%% Residuals took in seconds:',T_RES
         TOTTIME = MPI_WTIME() - TOTTIME
c         IF( IAM.EQ.0 ) WRITE(*,*)
c     $      ' %%% Total execution time in seconds:', TOTTIME
*
*
*        Print results to screen.
*
.GT..OR..GT.       IF( (ORTHTHRESH)(RNORMTHRESH) ) THEN
          PASSED = 'FAILED'
       ELSE
          PASSED = 'PASSED'
       END IF
         IF( DEBUG ) WRITE(*,*) '% #', IAM, ': Print results...'
.EQ.         IF( IAM0 ) THEN
            WRITE( NOUT, FMT = 9993 ) N, NB, NPROW, NPCOL, T_QR, PASSED
         END IF
         CALL BLACS_BARRIER( ICTXT, 'All' )
      END DO
      END DO
      END DO
 999  CONTINUE
*
*     Deallocate MEM and IMEM.
*
      DEALLOCATE( MEM, IMEM )
*
      CALL BLACS_GRIDEXIT( ICTXT )
*
 777  CONTINUE
*
      CALL BLACS_EXIT( 0 )
*
*     Format specifications.
*
 6666 FORMAT(A2,A3,A6,A4,A5,A6,A3,A3,A3,A9,A9,A9,A8,A8,A9,A9,A9,A9,A9,
     $       A9,A9,A9,A9,A9,A9,A5,A5,A8,A5,A5)
 7777 FORMAT(A2,I3,I6,I4,I5,I6,I3,I3,I3,F9.2,F9.2,F9.2,F8.2,F8.2,F9.2,
     $       F9.2,F9.2,F9.2,F9.2,F9.2,F9.2,F9.2,E9.2,E9.2,E9.2,I5,I5,
     $       F8.4,I5,I5,A2)
 9995 FORMAT( '    N  NB    P    Q  QR Time  CHECK' )
 9994 FORMAT( '----- --- ---- ---- -------- ------' )
 9993 FORMAT( I5, 1X, I3, 1X, I4, 1X, I4, 1X, F8.2, 1X, A6 )
          
*
      END