ddrgsx.f

famous linear algebra library (LAPACK) ports to windows

                        NERRS = NERRS + 1
                        IF( RESULT( J ).LT.10000.0D0 ) THEN
                           WRITE( NOUT, FMT = 9991 )MPLUSN, PRTYPE,
     $                        WEIGHT, M, J, RESULT( J )
                        ELSE
                           WRITE( NOUT, FMT = 9990 )MPLUSN, PRTYPE,
     $                        WEIGHT, M, J, RESULT( J )
                        END IF
                     END IF
   20             CONTINUE
*
   30          CONTINUE
   40       CONTINUE
   50    CONTINUE
   60 CONTINUE
*
      GO TO 150
*
   70 CONTINUE
*
*     Read in data from file to check accuracy of condition estimation
*     Read input data until N=0
*
      NPTKNT = 0
*
   80 CONTINUE
      READ( NIN, FMT = *, END = 140 )MPLUSN
      IF( MPLUSN.EQ.0 )
     $   GO TO 140
      READ( NIN, FMT = *, END = 140 )N
      DO 90 I = 1, MPLUSN
         READ( NIN, FMT = * )( AI( I, J ), J = 1, MPLUSN )
   90 CONTINUE
      DO 100 I = 1, MPLUSN
         READ( NIN, FMT = * )( BI( I, J ), J = 1, MPLUSN )
  100 CONTINUE
      READ( NIN, FMT = * )PLTRU, DIFTRU
*
      NPTKNT = NPTKNT + 1
      FS = .TRUE.
      K = 0
      M = MPLUSN - N
*
      CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
      CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
*
*     Compute the Schur factorization while swaping the
*     m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
*
      CALL DGGESX( 'V', 'V', 'S', DLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA,
     $             MM, ALPHAR, ALPHAI, BETA, Q, LDA, Z, LDA, PL, DIFEST,
     $             WORK, LWORK, IWORK, LIWORK, BWORK, LINFO )
*
      IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
         RESULT( 1 ) = ULPINV
         WRITE( NOUT, FMT = 9998 )'DGGESX', LINFO, MPLUSN, NPTKNT
         GO TO 130
      END IF
*
*     Compute the norm(A, B)
*        (should this be norm of (A,B) or (AI,BI)?)
*
      CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, MPLUSN )
      CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
     $             WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
      ABNRM = DLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, WORK )
*
*     Do tests (1) to (4)
*
      CALL DGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, LDA, WORK,
     $             RESULT( 1 ) )
      CALL DGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, LDA, WORK,
     $             RESULT( 2 ) )
      CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, LDA, WORK,
     $             RESULT( 3 ) )
      CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, LDA, WORK,
     $             RESULT( 4 ) )
*
*     Do tests (5) and (6): check Schur form of A and compare
*     eigenvalues with diagonals.
*
      NTEST = 6
      TEMP1 = ZERO
      RESULT( 5 ) = ZERO
      RESULT( 6 ) = ZERO
*
      DO 110 J = 1, MPLUSN
         ILABAD = .FALSE.
         IF( ALPHAI( J ).EQ.ZERO ) THEN
            TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
     $              MAX( SMLNUM, ABS( ALPHAR( J ) ), ABS( AI( J,
     $              J ) ) )+ABS( BETA( J )-BI( J, J ) ) /
     $              MAX( SMLNUM, ABS( BETA( J ) ), ABS( BI( J, J ) ) ) )
     $               / ULP
            IF( J.LT.MPLUSN ) THEN
               IF( AI( J+1, J ).NE.ZERO ) THEN
                  ILABAD = .TRUE.
                  RESULT( 5 ) = ULPINV
               END IF
            END IF
            IF( J.GT.1 ) THEN
               IF( AI( J, J-1 ).NE.ZERO ) THEN
                  ILABAD = .TRUE.
                  RESULT( 5 ) = ULPINV
               END IF
            END IF
         ELSE
            IF( ALPHAI( J ).GT.ZERO ) THEN
               I1 = J
            ELSE
               I1 = J - 1
            END IF
            IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
               ILABAD = .TRUE.
            ELSE IF( I1.LT.MPLUSN-1 ) THEN
               IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
                  ILABAD = .TRUE.
                  RESULT( 5 ) = ULPINV
               END IF
            ELSE IF( I1.GT.1 ) THEN
               IF( AI( I1, I1-1 ).NE.ZERO ) THEN
                  ILABAD = .TRUE.
                  RESULT( 5 ) = ULPINV
               END IF
            END IF
            IF( .NOT.ILABAD ) THEN
               CALL DGET53( AI( I1, I1 ), LDA, BI( I1, I1 ), LDA,
     $                      BETA( J ), ALPHAR( J ), ALPHAI( J ), TEMP2,
     $                      IINFO )
               IF( IINFO.GE.3 ) THEN
                  WRITE( NOUT, FMT = 9997 )IINFO, J, MPLUSN, NPTKNT
                  INFO = ABS( IINFO )
               END IF
            ELSE
               TEMP2 = ULPINV
            END IF
         END IF
         TEMP1 = MAX( TEMP1, TEMP2 )
         IF( ILABAD ) THEN
            WRITE( NOUT, FMT = 9996 )J, MPLUSN, NPTKNT
         END IF
  110 CONTINUE
      RESULT( 6 ) = TEMP1
*
*     Test (7) (if sorting worked)  <--------- need to be checked.
*
      NTEST = 7
      RESULT( 7 ) = ZERO
      IF( LINFO.EQ.MPLUSN+3 )
     $   RESULT( 7 ) = ULPINV
*
*     Test (8): compare the estimated value of DIF and its true value.
*
      NTEST = 8
      RESULT( 8 ) = ZERO
      IF( DIFEST( 2 ).EQ.ZERO ) THEN
         IF( DIFTRU.GT.ABNRM*ULP )
     $      RESULT( 8 ) = ULPINV
      ELSE IF( DIFTRU.EQ.ZERO ) THEN
         IF( DIFEST( 2 ).GT.ABNRM*ULP )
     $      RESULT( 8 ) = ULPINV
      ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
     $         ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
         RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), DIFEST( 2 ) / DIFTRU )
      END IF
*
*     Test (9)
*
      NTEST = 9
      RESULT( 9 ) = ZERO
      IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
         IF( DIFTRU.GT.ABNRM*ULP )
     $      RESULT( 9 ) = ULPINV
         IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
     $      RESULT( 9 ) = ULPINV
         IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
     $      RESULT( 9 ) = ULPINV
      END IF
*
*     Test (10): compare the estimated value of PL and it true value.
*
      NTEST = 10
      RESULT( 10 ) = ZERO
      IF( PL( 1 ).EQ.ZERO ) THEN
         IF( PLTRU.GT.ABNRM*ULP )
     $      RESULT( 10 ) = ULPINV
      ELSE IF( PLTRU.EQ.ZERO ) THEN
         IF( PL( 1 ).GT.ABNRM*ULP )
     $      RESULT( 10 ) = ULPINV
      ELSE IF( ( PLTRU.GT.THRESH*PL( 1 ) ) .OR.
     $         ( PLTRU*THRESH.LT.PL( 1 ) ) ) THEN
         RESULT( 10 ) = ULPINV
      END IF
*
      NTESTT = NTESTT + NTEST
*
*     Print out tests which fail.
*
      DO 120 J = 1, NTEST
         IF( RESULT( J ).GE.THRESH ) THEN
*
*           If this is the first test to fail,
*           print a header to the data file.
*
            IF( NERRS.EQ.0 ) THEN
               WRITE( NOUT, FMT = 9995 )'SGX'
*
*              Matrix types
*
               WRITE( NOUT, FMT = 9994 )
*
*              Tests performed
*
               WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
     $            'transpose', ( '''', I = 1, 4 )
*
            END IF
            NERRS = NERRS + 1
            IF( RESULT( J ).LT.10000.0D0 ) THEN
               WRITE( NOUT, FMT = 9989 )NPTKNT, MPLUSN, J, RESULT( J )
            ELSE
               WRITE( NOUT, FMT = 9988 )NPTKNT, MPLUSN, J, RESULT( J )
            END IF
         END IF
*
  120 CONTINUE
*
  130 CONTINUE
      GO TO 80
  140 CONTINUE
*
  150 CONTINUE
*
*     Summary
*
      CALL ALASVM( 'SGX', NOUT, NERRS, NTESTT, 0 )
*
      WORK( 1 ) = MAXWRK
*
      RETURN
*
 9999 FORMAT( ' DDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
     $      I6, ', JTYPE=', I6, ')' )
*
 9998 FORMAT( ' DDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
     $      I6, ', Input Example #', I2, ')' )
*
 9997 FORMAT( ' DDRGSX: DGET53 returned INFO=', I1, ' for eigenvalue ',
     $      I6, '.', / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
*
 9996 FORMAT( ' DDRGSX: S not in Schur form at eigenvalue ', I6, '.',
     $      / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
*
 9995 FORMAT( / 1X, A3, ' -- Real Expert Generalized Schur form',
     $      ' problem driver' )
*
 9994 FORMAT( 'Input Example' )
*
 9993 FORMAT( ' Matrix types: ', /
     $      '  1:  A is a block diagonal matrix of Jordan blocks ',
     $      'and B is the identity ', / '      matrix, ',
     $      / '  2:  A and B are upper triangular matrices, ',
     $      / '  3:  A and B are as type 2, but each second diagonal ',
     $      'block in A_11 and ', /
     $      '      each third diaongal block in A_22 are 2x2 blocks,',
     $      / '  4:  A and B are block diagonal matrices, ',
     $      / '  5:  (A,B) has potentially close or common ',
     $      'eigenvalues.', / )
*
 9992 FORMAT( / ' Tests performed:  (S is Schur, T is triangular, ',
     $      'Q and Z are ', A, ',', / 19X,
     $      ' a is alpha, b is beta, and ', A, ' means ', A, '.)',
     $      / '  1 = | A - Q S Z', A,
     $      ' | / ( |A| n ulp )      2 = | B - Q T Z', A,
     $      ' | / ( |B| n ulp )', / '  3 = | I - QQ', A,
     $      ' | / ( n ulp )             4 = | I - ZZ', A,
     $      ' | / ( n ulp )', / '  5 = 1/ULP  if A is not in ',
     $      'Schur form S', / '  6 = difference between (alpha,beta)',
     $      ' and diagonals of (S,T)', /
     $      '  7 = 1/ULP  if SDIM is not the correct number of ',
     $      'selected eigenvalues', /
     $      '  8 = 1/ULP  if DIFEST/DIFTRU > 10*THRESH or ',
     $      'DIFTRU/DIFEST > 10*THRESH',
     $      / '  9 = 1/ULP  if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
     $      'when reordering fails', /
     $      ' 10 = 1/ULP  if PLEST/PLTRU > THRESH or ',
     $      'PLTRU/PLEST > THRESH', /
     $      '    ( Test 10 is only for input examples )', / )
 9991 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.4,
     $      ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, F8.2 )
 9990 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.4,
     $      ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, D10.4 )
 9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
     $      ' result ', I2, ' is', 0P, F8.2 )
 9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
     $      ' result ', I2, ' is', 1P, D10.3 )
*
*     End of DDRGSX
*
      END