sdrves.f

来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 843 行 · 第 1/3 页

     $                         LDA )
               END IF
*
            ELSE IF( ITYPE.EQ.10 ) THEN
*
*              Triangular, random eigenvalues
*
               CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE,
     $                      'T', 'N', WORK( N+1 ), 1, ONE,
     $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0,
     $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
*
            ELSE
*
               IINFO = 1
            END IF
*
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9992 )'Generator', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               RETURN
            END IF
*
   90       CONTINUE
*
*           Test for minimal and generous workspace
*
            DO 250 IWK = 1, 2
               IF( IWK.EQ.1 ) THEN
                  NNWORK = 3*N
               ELSE
                  NNWORK = 5*N + 2*N**2
               END IF
               NNWORK = MAX( NNWORK, 1 )
*
*              Initialize RESULT
*
               DO 100 J = 1, 13
                  RESULT( J ) = -ONE
  100          CONTINUE
*
*              Test with and without sorting of eigenvalues
*
               DO 210 ISORT = 0, 1
                  IF( ISORT.EQ.0 ) THEN
                     SORT = 'N'
                     RSUB = 0
                  ELSE
                     SORT = 'S'
                     RSUB = 6
                  END IF
*
*                 Compute Schur form and Schur vectors, and test them
*
                  CALL SLACPY( 'F', N, N, A, LDA, H, LDA )
                  CALL SGEES( 'V', SORT, SSLECT, N, H, LDA, SDIM, WR,
     $                        WI, VS, LDVS, WORK, NNWORK, BWORK, IINFO )
                  IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
                     RESULT( 1+RSUB ) = ULPINV
                     WRITE( NOUNIT, FMT = 9992 )'SGEES1', IINFO, N,
     $                  JTYPE, IOLDSD
                     INFO = ABS( IINFO )
                     GO TO 220
                  END IF
*
*                 Do Test (1) or Test (7)
*
                  RESULT( 1+RSUB ) = ZERO
                  DO 120 J = 1, N - 2
                     DO 110 I = J + 2, N
                        IF( H( I, J ).NE.ZERO )
     $                     RESULT( 1+RSUB ) = ULPINV
  110                CONTINUE
  120             CONTINUE
                  DO 130 I = 1, N - 2
                     IF( H( I+1, I ).NE.ZERO .AND. H( I+2, I+1 ).NE.
     $                   ZERO )RESULT( 1+RSUB ) = ULPINV
  130             CONTINUE
                  DO 140 I = 1, N - 1
                     IF( H( I+1, I ).NE.ZERO ) THEN
                        IF( H( I, I ).NE.H( I+1, I+1 ) .OR.
     $                      H( I, I+1 ).EQ.ZERO .OR.
     $                      SIGN( ONE, H( I+1, I ) ).EQ.
     $                      SIGN( ONE, H( I, I+1 ) ) )RESULT( 1+RSUB )
     $                      = ULPINV
                     END IF
  140             CONTINUE
*
*                 Do Tests (2) and (3) or Tests (8) and (9)
*
                  LWORK = MAX( 1, 2*N*N )
                  CALL SHST01( N, 1, N, A, LDA, H, LDA, VS, LDVS, WORK,
     $                         LWORK, RES )
                  RESULT( 2+RSUB ) = RES( 1 )
                  RESULT( 3+RSUB ) = RES( 2 )
*
*                 Do Test (4) or Test (10)
*
                  RESULT( 4+RSUB ) = ZERO
                  DO 150 I = 1, N
                     IF( H( I, I ).NE.WR( I ) )
     $                  RESULT( 4+RSUB ) = ULPINV
  150             CONTINUE
                  IF( N.GT.1 ) THEN
                     IF( H( 2, 1 ).EQ.ZERO .AND. WI( 1 ).NE.ZERO )
     $                  RESULT( 4+RSUB ) = ULPINV
                     IF( H( N, N-1 ).EQ.ZERO .AND. WI( N ).NE.ZERO )
     $                  RESULT( 4+RSUB ) = ULPINV
                  END IF
                  DO 160 I = 1, N - 1
                     IF( H( I+1, I ).NE.ZERO ) THEN
                        TMP = SQRT( ABS( H( I+1, I ) ) )*
     $                        SQRT( ABS( H( I, I+1 ) ) )
                        RESULT( 4+RSUB ) = MAX( RESULT( 4+RSUB ),
     $                                     ABS( WI( I )-TMP ) /
     $                                     MAX( ULP*TMP, UNFL ) )
                        RESULT( 4+RSUB ) = MAX( RESULT( 4+RSUB ),
     $                                     ABS( WI( I+1 )+TMP ) /
     $                                     MAX( ULP*TMP, UNFL ) )
                     ELSE IF( I.GT.1 ) THEN
                        IF( H( I+1, I ).EQ.ZERO .AND. H( I, I-1 ).EQ.
     $                      ZERO .AND. WI( I ).NE.ZERO )RESULT( 4+RSUB )
     $                       = ULPINV
                     END IF
  160             CONTINUE
*
*                 Do Test (5) or Test (11)
*
                  CALL SLACPY( 'F', N, N, A, LDA, HT, LDA )
                  CALL SGEES( 'N', SORT, SSLECT, N, HT, LDA, SDIM, WRT,
     $                        WIT, VS, LDVS, WORK, NNWORK, BWORK,
     $                        IINFO )
                  IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
                     RESULT( 5+RSUB ) = ULPINV
                     WRITE( NOUNIT, FMT = 9992 )'SGEES2', IINFO, N,
     $                  JTYPE, IOLDSD
                     INFO = ABS( IINFO )
                     GO TO 220
                  END IF
*
                  RESULT( 5+RSUB ) = ZERO
                  DO 180 J = 1, N
                     DO 170 I = 1, N
                        IF( H( I, J ).NE.HT( I, J ) )
     $                     RESULT( 5+RSUB ) = ULPINV
  170                CONTINUE
  180             CONTINUE
*
*                 Do Test (6) or Test (12)
*
                  RESULT( 6+RSUB ) = ZERO
                  DO 190 I = 1, N
                     IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
     $                  RESULT( 6+RSUB ) = ULPINV
  190             CONTINUE
*
*                 Do Test (13)
*
                  IF( ISORT.EQ.1 ) THEN
                     RESULT( 13 ) = ZERO
                     KNTEIG = 0
                     DO 200 I = 1, N
                        IF( SSLECT( WR( I ), WI( I ) ) .OR.
     $                      SSLECT( WR( I ), -WI( I ) ) )
     $                      KNTEIG = KNTEIG + 1
                        IF( I.LT.N ) THEN
                           IF( ( SSLECT( WR( I+1 ),
     $                         WI( I+1 ) ) .OR. SSLECT( WR( I+1 ),
     $                         -WI( I+1 ) ) ) .AND.
     $                         ( .NOT.( SSLECT( WR( I ),
     $                         WI( I ) ) .OR. SSLECT( WR( I ),
     $                         -WI( I ) ) ) ) .AND. IINFO.NE.N+2 )
     $                         RESULT( 13 ) = ULPINV
                        END IF
  200                CONTINUE
                     IF( SDIM.NE.KNTEIG ) THEN
                        RESULT( 13 ) = ULPINV
                     END IF
                  END IF
*
  210          CONTINUE
*
*              End of Loop -- Check for RESULT(j) > THRESH
*
  220          CONTINUE
*
               NTEST = 0
               NFAIL = 0
               DO 230 J = 1, 13
                  IF( RESULT( J ).GE.ZERO )
     $               NTEST = NTEST + 1
                  IF( RESULT( J ).GE.THRESH )
     $               NFAIL = NFAIL + 1
  230          CONTINUE
*
               IF( NFAIL.GT.0 )
     $            NTESTF = NTESTF + 1
               IF( NTESTF.EQ.1 ) THEN
                  WRITE( NOUNIT, FMT = 9999 )PATH
                  WRITE( NOUNIT, FMT = 9998 )
                  WRITE( NOUNIT, FMT = 9997 )
                  WRITE( NOUNIT, FMT = 9996 )
                  WRITE( NOUNIT, FMT = 9995 )THRESH
                  WRITE( NOUNIT, FMT = 9994 )
                  NTESTF = 2
               END IF
*
               DO 240 J = 1, 13
                  IF( RESULT( J ).GE.THRESH ) THEN
                     WRITE( NOUNIT, FMT = 9993 )N, IWK, IOLDSD, JTYPE,
     $                  J, RESULT( J )
                  END IF
  240          CONTINUE
*
               NERRS = NERRS + NFAIL
               NTESTT = NTESTT + NTEST
*
  250       CONTINUE
  260    CONTINUE
  270 CONTINUE
*
*     Summary
*
      CALL SLASUM( PATH, NOUNIT, NERRS, NTESTT )
*
 9999 FORMAT( / 1X, A3, ' -- Real Schur Form Decomposition Driver',
     $      / ' Matrix types (see SDRVES for details): ' )
*
 9998 FORMAT( / ' Special Matrices:', / '  1=Zero matrix.             ',
     $      '           ', '  5=Diagonal: geometr. spaced entries.',
     $      / '  2=Identity matrix.                    ', '  6=Diagona',
     $      'l: clustered entries.', / '  3=Transposed Jordan block.  ',
     $      '          ', '  7=Diagonal: large, evenly spaced.', / '  ',
     $      '4=Diagonal: evenly spaced entries.    ', '  8=Diagonal: s',
     $      'mall, evenly spaced.' )
 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / '  9=Well-cond., ev',
     $      'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
     $      'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
     $      ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
     $      'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
     $      'lex ', / ' 12=Well-cond., random complex ', 6X, '   ',
     $      ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
     $      'tioned, evenly spaced.     ', ' 18=Ill-cond., small rand.',
     $      ' complx ' )
 9996 FORMAT( ' 19=Matrix with random O(1) entries.    ', ' 21=Matrix ',
     $      'with small random entries.', / ' 20=Matrix with large ran',
     $      'dom entries.   ', / )
 9995 FORMAT( ' Tests performed with test threshold =', F8.2,
     $      / ' ( A denotes A on input and T denotes A on output)',
     $      / / ' 1 = 0 if T in Schur form (no sort), ',
     $      '  1/ulp otherwise', /
     $      ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
     $      / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ', /
     $      ' 4 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (no sort),',
     $      '  1/ulp otherwise', /
     $      ' 5 = 0 if T same no matter if VS computed (no sort),',
     $      '  1/ulp otherwise', /
     $      ' 6 = 0 if WR, WI same no matter if VS computed (no sort)',
     $      ',  1/ulp otherwise' )
 9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', '  1/ulp otherwise',
     $      / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
     $      / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
     $      / ' 10 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (sort),',
     $      '  1/ulp otherwise', /
     $      ' 11 = 0 if T same no matter if VS computed (sort),',
     $      '  1/ulp otherwise', /
     $      ' 12 = 0 if WR, WI same no matter if VS computed (sort),',
     $      '  1/ulp otherwise', /
     $      ' 13 = 0 if sorting succesful, 1/ulp otherwise', / )
 9993 FORMAT( ' N=', I5, ', IWK=', I2, ', seed=', 4( I4, ',' ),
     $      ' type ', I2, ', test(', I2, ')=', G10.3 )
 9992 FORMAT( ' SDRVES: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
     $      I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
*
      RETURN
*
*     End of SDRVES
*
      END