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
            IF (intallc.LT.imem(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 )
            IF (info.NE.0) 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
            IF( i.GT.1 ) 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 )
            IF( i.LT.n-1 ) THEN
               CALL pdelget( 'All', '1-Tree', elem3, mem(ipa), i+2, i+1,
     $              desca )
            ELSE
               elem3 = zero
            END IF
            IF( elem2.NE.zero .AND. abs(elem1)+abs(elem2)+abs(elem3).GT.
     $         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 )
            IF( anorm.GT.zero )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.
*
           IF( (orth.GT.thresh).OR.(rnorm.GT.thresh) ) THEN
              passed = 'FAILED'
           ELSE
              passed = 'PASSED'
           END IF
         IF( debug ) WRITE(*,*) '% #', iam, ': Print results...'
         IF( iam.EQ.0 ) 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