sdrges.f

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

                  CALL SORM2R( 'R', 'T', N, N, N-1, Z, LDQ, WORK( N+1 ),
     $                         B, LDA, WORK( 2*N+1 ), IINFO )
                  IF( IINFO.NE.0 )
     $               GO TO 100
               END IF
            ELSE
*
*              Random matrices
*
               DO 90 JC = 1, N
                  DO 80 JR = 1, N
                     A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )*
     $                             SLARND( 2, ISEED )
                     B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )*
     $                             SLARND( 2, ISEED )
   80             CONTINUE
   90          CONTINUE
            END IF
*
  100       CONTINUE
*
            IF( IINFO.NE.0 ) THEN
               WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
     $            IOLDSD
               INFO = ABS( IINFO )
               RETURN
            END IF
*
  110       CONTINUE
*
            DO 120 I = 1, 13
               RESULT( I ) = -ONE
  120       CONTINUE
*
*           Test with and without sorting of eigenvalues
*
            DO 150 ISORT = 0, 1
               IF( ISORT.EQ.0 ) THEN
                  SORT = 'N'
                  RSUB = 0
               ELSE
                  SORT = 'S'
                  RSUB = 5
               END IF
*
*              Call SGGES to compute H, T, Q, Z, alpha, and beta.
*
               CALL SLACPY( 'Full', N, N, A, LDA, S, LDA )
               CALL SLACPY( 'Full', N, N, B, LDA, T, LDA )
               NTEST = 1 + RSUB + ISORT
               RESULT( 1+RSUB+ISORT ) = ULPINV
               CALL SGGES( 'V', 'V', SORT, SLCTES, N, S, LDA, T, LDA,
     $                     SDIM, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDQ,
     $                     WORK, LWORK, BWORK, IINFO )
               IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
                  RESULT( 1+RSUB+ISORT ) = ULPINV
                  WRITE( NOUNIT, FMT = 9999 )'SGGES', IINFO, N, JTYPE,
     $               IOLDSD
                  INFO = ABS( IINFO )
                  GO TO 160
               END IF
*
               NTEST = 4 + RSUB
*
*              Do tests 1--4 (or tests 7--9 when reordering )
*
               IF( ISORT.EQ.0 ) THEN
                  CALL SGET51( 1, N, A, LDA, S, LDA, Q, LDQ, Z, LDQ,
     $                         WORK, RESULT( 1 ) )
                  CALL SGET51( 1, N, B, LDA, T, LDA, Q, LDQ, Z, LDQ,
     $                         WORK, RESULT( 2 ) )
               ELSE
                  CALL SGET54( N, A, LDA, B, LDA, S, LDA, T, LDA, Q,
     $                         LDQ, Z, LDQ, WORK, RESULT( 7 ) )
               END IF
               CALL SGET51( 3, N, A, LDA, T, LDA, Q, LDQ, Q, LDQ, WORK,
     $                      RESULT( 3+RSUB ) )
               CALL SGET51( 3, N, B, LDA, T, LDA, Z, LDQ, Z, LDQ, WORK,
     $                      RESULT( 4+RSUB ) )
*
*              Do test 5 and 6 (or Tests 10 and 11 when reordering):
*              check Schur form of A and compare eigenvalues with
*              diagonals.
*
               NTEST = 6 + RSUB
               TEMP1 = ZERO
*
               DO 130 J = 1, N
                  ILABAD = .FALSE.
                  IF( ALPHAI( J ).EQ.ZERO ) THEN
                     TEMP2 = ( ABS( ALPHAR( J )-S( J, J ) ) /
     $                       MAX( SAFMIN, ABS( ALPHAR( J ) ), ABS( S( J,
     $                       J ) ) )+ABS( BETA( J )-T( J, J ) ) /
     $                       MAX( SAFMIN, ABS( BETA( J ) ), ABS( T( J,
     $                       J ) ) ) ) / ULP
*
                     IF( J.LT.N ) THEN
                        IF( S( J+1, J ).NE.ZERO ) THEN
                           ILABAD = .TRUE.
                           RESULT( 5+RSUB ) = ULPINV
                        END IF
                     END IF
                     IF( J.GT.1 ) THEN
                        IF( S( J, J-1 ).NE.ZERO ) THEN
                           ILABAD = .TRUE.
                           RESULT( 5+RSUB ) = 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.N ) THEN
                        ILABAD = .TRUE.
                     ELSE IF( I1.LT.N-1 ) THEN
                        IF( S( I1+2, I1+1 ).NE.ZERO ) THEN
                           ILABAD = .TRUE.
                           RESULT( 5+RSUB ) = ULPINV
                        END IF
                     ELSE IF( I1.GT.1 ) THEN
                        IF( S( I1, I1-1 ).NE.ZERO ) THEN
                           ILABAD = .TRUE.
                           RESULT( 5+RSUB ) = ULPINV
                        END IF
                     END IF
                     IF( .NOT.ILABAD ) THEN
                        CALL SGET53( S( I1, I1 ), LDA, T( I1, I1 ), LDA,
     $                               BETA( J ), ALPHAR( J ),
     $                               ALPHAI( J ), TEMP2, IERR )
                        IF( IERR.GE.3 ) THEN
                           WRITE( NOUNIT, FMT = 9998 )IERR, J, N,
     $                        JTYPE, IOLDSD
                           INFO = ABS( IERR )
                        END IF
                     ELSE
                        TEMP2 = ULPINV
                     END IF
*
                  END IF
                  TEMP1 = MAX( TEMP1, TEMP2 )
                  IF( ILABAD ) THEN
                     WRITE( NOUNIT, FMT = 9997 )J, N, JTYPE, IOLDSD
                  END IF
  130          CONTINUE
               RESULT( 6+RSUB ) = TEMP1
*
               IF( ISORT.GE.1 ) THEN
*
*                 Do test 12
*
                  NTEST = 12
                  RESULT( 12 ) = ZERO
                  KNTEIG = 0
                  DO 140 I = 1, N
                     IF( SLCTES( ALPHAR( I ), ALPHAI( I ),
     $                   BETA( I ) ) .OR. SLCTES( ALPHAR( I ),
     $                   -ALPHAI( I ), BETA( I ) ) ) THEN
                        KNTEIG = KNTEIG + 1
                     END IF
                     IF( I.LT.N ) THEN
                        IF( ( SLCTES( ALPHAR( I+1 ), ALPHAI( I+1 ),
     $                      BETA( I+1 ) ) .OR. SLCTES( ALPHAR( I+1 ),
     $                      -ALPHAI( I+1 ), BETA( I+1 ) ) ) .AND.
     $                      ( .NOT.( SLCTES( ALPHAR( I ), ALPHAI( I ),
     $                      BETA( I ) ) .OR. SLCTES( ALPHAR( I ),
     $                      -ALPHAI( I ), BETA( I ) ) ) ) .AND.
     $                      IINFO.NE.N+2 ) THEN
                           RESULT( 12 ) = ULPINV
                        END IF
                     END IF
  140             CONTINUE
                  IF( SDIM.NE.KNTEIG ) THEN
                     RESULT( 12 ) = ULPINV
                  END IF
               END IF
*
  150       CONTINUE
*
*           End of Loop -- Check for RESULT(j) > THRESH
*
  160       CONTINUE
*
            NTESTT = NTESTT + NTEST
*
*           Print out tests which fail.
*
            DO 170 JR = 1, NTEST
               IF( RESULT( JR ).GE.THRESH ) THEN
*
*                 If this is the first test to fail,
*                 print a header to the data file.
*
                  IF( NERRS.EQ.0 ) THEN
                     WRITE( NOUNIT, FMT = 9996 )'SGS'
*
*                    Matrix types
*
                     WRITE( NOUNIT, FMT = 9995 )
                     WRITE( NOUNIT, FMT = 9994 )
                     WRITE( NOUNIT, FMT = 9993 )'Orthogonal'
*
*                    Tests performed
*
                     WRITE( NOUNIT, FMT = 9992 )'orthogonal', '''',
     $                  'transpose', ( '''', J = 1, 8 )
*
                  END IF
                  NERRS = NERRS + 1
                  IF( RESULT( JR ).LT.10000.0 ) THEN
                     WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR,
     $                  RESULT( JR )
                  ELSE
                     WRITE( NOUNIT, FMT = 9990 )N, JTYPE, IOLDSD, JR,
     $                  RESULT( JR )
                  END IF
               END IF
  170       CONTINUE
*
  180    CONTINUE
  190 CONTINUE
*
*     Summary
*
      CALL ALASVM( 'SGS', NOUNIT, NERRS, NTESTT, 0 )
*
      WORK( 1 ) = MAXWRK
*
      RETURN
*
 9999 FORMAT( ' SDRGES: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
     $      I6, ', JTYPE=', I6, ', ISEED=(', 4( I4, ',' ), I5, ')' )
*
 9998 FORMAT( ' SDRGES: SGET53 returned INFO=', I1, ' for eigenvalue ',
     $      I6, '.', / 9X, 'N=', I6, ', JTYPE=', I6, ', ISEED=(',
     $      4( I4, ',' ), I5, ')' )
*
 9997 FORMAT( ' SDRGES: S not in Schur form at eigenvalue ', I6, '.',
     $      / 9X, 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ),
     $      I5, ')' )
*
 9996 FORMAT( / 1X, A3, ' -- Real Generalized Schur form driver' )
*
 9995 FORMAT( ' Matrix types (see SDRGES for details): ' )
*
 9994 FORMAT( ' Special Matrices:', 23X,
     $      '(J''=transposed Jordan block)',
     $      / '   1=(0,0)  2=(I,0)  3=(0,I)  4=(I,I)  5=(J'',J'')  ',
     $      '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices:  ( ',
     $      'D=diag(0,1,2,...) )', / '   7=(D,I)   9=(large*D, small*I',
     $      ')  11=(large*I, small*D)  13=(large*D, large*I)', /
     $      '   8=(I,D)  10=(small*D, large*I)  12=(small*I, large*D) ',
     $      ' 14=(small*D, small*I)', / '  15=(D, reversed D)' )
 9993 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:',
     $      / '  16=Transposed Jordan Blocks             19=geometric ',
     $      'alpha, beta=0,1', / '  17=arithm. alpha&beta             ',
     $      '      20=arithmetic alpha, beta=0,1', / '  18=clustered ',
     $      'alpha, beta=0,1            21=random alpha, beta=0,1',
     $      / ' Large & Small Matrices:', / '  22=(large, small)   ',
     $      '23=(small,large)    24=(small,small)    25=(large,large)',
     $      / '  26=random O(1) matrices.' )
*
 9992 FORMAT( / ' Tests performed:  (S is Schur, T is triangular, ',
     $      'Q and Z are ', A, ',', / 19X,
     $      'l and r are the appropriate left and right', / 19X,
     $      'eigenvectors, resp., a is alpha, b is beta, and', / 19X, A,
     $      ' means ', A, '.)', / ' Without ordering: ',
     $      / '  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 = A is in Schur form S',
     $      / '  6 = difference between (alpha,beta)',
     $      ' and diagonals of (S,T)', / ' With ordering: ',
     $      / '  7 = | (A,B) - Q (S,T) Z', A,
     $      ' | / ( |(A,B)| n ulp )  ', / '  8 = | I - QQ', A,
     $      ' | / ( n ulp )            9 = | I - ZZ', A,
     $      ' | / ( n ulp )', / ' 10 = A is in Schur form S',
     $      / ' 11 = difference between (alpha,beta) and diagonals',
     $      ' of (S,T)', / ' 12 = SDIM is the correct number of ',
     $      'selected eigenvalues', / )
 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
     $      4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 )
 9990 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
     $      4( I4, ',' ), ' result ', I2, ' is', 1P, E10.3 )
*
*     End of SDRGES
*
      END