cdrgsx.f

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

*     ..
*     .. Scalars in Common ..
      LOGICAL            FS
      INTEGER            K, M, MPLUSN, N
*     ..
*     .. Common blocks ..
      COMMON             / MN / M, N, MPLUSN, K, FS
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, AIMAG, MAX, REAL, SQRT
*     ..
*     .. Statement Functions ..
      REAL               ABS1
*     ..
*     .. Statement Function definitions ..
      ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
*     ..
*     .. Executable Statements ..
*
*     Check for errors
*
      IF( NSIZE.LT.0 ) THEN
         INFO = -1
      ELSE IF( THRESH.LT.ZERO ) THEN
         INFO = -2
      ELSE IF( NIN.LE.0 ) THEN
         INFO = -3
      ELSE IF( NOUT.LE.0 ) THEN
         INFO = -4
      ELSE IF( LDA.LT.1 .OR. LDA.LT.NSIZE ) THEN
         INFO = -6
      ELSE IF( LDC.LT.1 .OR. LDC.LT.NSIZE*NSIZE / 2 ) THEN
         INFO = -15
      ELSE IF( LIWORK.LT.NSIZE+2 ) THEN
         INFO = -21
      END IF
*
*     Compute workspace
*      (Note: Comments in the code beginning "Workspace:" describe the
*       minimal amount of workspace needed at that point in the code,
*       as well as the preferred amount for good performance.
*       NB refers to the optimal block size for the immediately
*       following subroutine, as returned by ILAENV.)
*
      MINWRK = 1
      IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
         MINWRK = 3*NSIZE*NSIZE / 2
*
*        workspace for cggesx
*
         MAXWRK = NSIZE*( 1+ILAENV( 1, 'CGEQRF', ' ', NSIZE, 1, NSIZE,
     $            0 ) )
         MAXWRK = MAX( MAXWRK, NSIZE*( 1+ILAENV( 1, 'CUNGQR', ' ',
     $            NSIZE, 1, NSIZE, -1 ) ) )
*
*        workspace for cgesvd
*
         BDSPAC = 3*NSIZE*NSIZE / 2
         MAXWRK = MAX( MAXWRK, NSIZE*NSIZE*
     $            ( 1+ILAENV( 1, 'CGEBRD', ' ', NSIZE*NSIZE / 2,
     $            NSIZE*NSIZE / 2, -1, -1 ) ) )
         MAXWRK = MAX( MAXWRK, BDSPAC )
*
         MAXWRK = MAX( MAXWRK, MINWRK )
*
         WORK( 1 ) = MAXWRK
      END IF
*
      IF( LWORK.LT.MINWRK )
     $   INFO = -18
*
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'CDRGSX', -INFO )
         RETURN
      END IF
*
*     Important constants
*
      ULP = SLAMCH( 'P' )
      ULPINV = ONE / ULP
      SMLNUM = SLAMCH( 'S' ) / ULP
      BIGNUM = ONE / SMLNUM
      CALL SLABAD( SMLNUM, BIGNUM )
      THRSH2 = TEN*THRESH
      NTESTT = 0
      NERRS = 0
*
*     Go to the tests for read-in matrix pairs
*
      IFUNC = 0
      IF( NSIZE.EQ.0 )
     $   GO TO 70
*
*     Test the built-in matrix pairs.
*     Loop over different functions (IFUNC) of CGGESX, types (PRTYPE)
*     of test matrices, different size (M+N)
*
      PRTYPE = 0
      QBA = 3
      QBB = 4
      WEIGHT = SQRT( ULP )
*
      DO 60 IFUNC = 0, 3
         DO 50 PRTYPE = 1, 5
            DO 40 M = 1, NSIZE - 1
               DO 30 N = 1, NSIZE - M
*
                  WEIGHT = ONE / WEIGHT
                  MPLUSN = M + N
*
*                 Generate test matrices
*
                  FS = .TRUE.
                  K = 0
*
                  CALL CLASET( 'Full', MPLUSN, MPLUSN, CZERO, CZERO, AI,
     $                         LDA )
                  CALL CLASET( 'Full', MPLUSN, MPLUSN, CZERO, CZERO, BI,
     $                         LDA )
*
                  CALL CLATM5( PRTYPE, M, N, AI, LDA, AI( M+1, M+1 ),
     $                         LDA, AI( 1, M+1 ), LDA, BI, LDA,
     $                         BI( M+1, M+1 ), LDA, BI( 1, M+1 ), LDA,
     $                         Q, LDA, Z, LDA, WEIGHT, QBA, QBB )
*
*                 Compute the Schur factorization and swapping the
*                 m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
*                 Swapping is accomplished via the function CLCTSX
*                 which is supplied below.
*
                  IF( IFUNC.EQ.0 ) THEN
                     SENSE = 'N'
                  ELSE IF( IFUNC.EQ.1 ) THEN
                     SENSE = 'E'
                  ELSE IF( IFUNC.EQ.2 ) THEN
                     SENSE = 'V'
                  ELSE IF( IFUNC.EQ.3 ) THEN
                     SENSE = 'B'
                  END IF
*
                  CALL CLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
                  CALL CLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
*
                  CALL CGGESX( 'V', 'V', 'S', CLCTSX, SENSE, MPLUSN, AI,
     $                         LDA, BI, LDA, MM, ALPHA, BETA, Q, LDA, Z,
     $                         LDA, PL, DIFEST, WORK, LWORK, RWORK,
     $                         IWORK, LIWORK, BWORK, LINFO )
*
                  IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
                     RESULT( 1 ) = ULPINV
                     WRITE( NOUT, FMT = 9999 )'CGGESX', LINFO, MPLUSN,
     $                  PRTYPE
                     INFO = LINFO
                     GO TO 30
                  END IF
*
*                 Compute the norm(A, B)
*
                  CALL CLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK,
     $                         MPLUSN )
                  CALL CLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
     $                         WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
                  ABNRM = CLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN,
     $                    RWORK )
*
*                 Do tests (1) to (4)
*
                  RESULT( 2 ) = ZERO
                  CALL CGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z,
     $                         LDA, WORK, RWORK, RESULT( 1 ) )
                  CALL CGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z,
     $                         LDA, WORK, RWORK, RESULT( 2 ) )
                  CALL CGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q,
     $                         LDA, WORK, RWORK, RESULT( 3 ) )
                  CALL CGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z,
     $                         LDA, WORK, RWORK, RESULT( 4 ) )
                  NTEST = 4
*
*                 Do tests (5) and (6): check Schur form of A and
*                 compare eigenvalues with diagonals.
*
                  TEMP1 = ZERO
                  RESULT( 5 ) = ZERO
                  RESULT( 6 ) = ZERO
*
                  DO 10 J = 1, MPLUSN
                     ILABAD = .FALSE.
                     TEMP2 = ( ABS1( ALPHA( J )-AI( J, J ) ) /
     $                       MAX( SMLNUM, ABS1( ALPHA( J ) ),
     $                       ABS1( AI( J, J ) ) )+
     $                       ABS1( BETA( J )-BI( J, J ) ) /
     $                       MAX( SMLNUM, ABS1( BETA( J ) ),
     $                       ABS1( 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
                     TEMP1 = MAX( TEMP1, TEMP2 )
                     IF( ILABAD ) THEN
                        WRITE( NOUT, FMT = 9997 )J, MPLUSN, PRTYPE
                     END IF
   10             CONTINUE
                  RESULT( 6 ) = TEMP1
                  NTEST = NTEST + 2
*
*                 Test (7) (if sorting worked)
*
                  RESULT( 7 ) = ZERO
                  IF( LINFO.EQ.MPLUSN+3 ) THEN
                     RESULT( 7 ) = ULPINV
                  ELSE IF( MM.NE.N ) THEN
                     RESULT( 7 ) = ULPINV
                  END IF
                  NTEST = NTEST + 1
*
*                 Test (8): compare the estimated value DIF and its
*                 value. first, compute the exact DIF.
*
                  RESULT( 8 ) = ZERO
                  MN2 = MM*( MPLUSN-MM )*2
                  IF( IFUNC.GE.2 .AND. MN2.LE.NCMAX*NCMAX ) THEN
*
*                    Note: for either following two cases, there are
*                    almost same number of test cases fail the test.
*
                     CALL CLAKF2( MM, MPLUSN-MM, AI, LDA,
     $                            AI( MM+1, MM+1 ), BI,
     $                            BI( MM+1, MM+1 ), C, LDC )
*
                     CALL CGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK,
     $                            1, WORK( 2 ), 1, WORK( 3 ), LWORK-2,
     $                            RWORK, INFO )
                     DIFTRU = S( MN2 )
*
                     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
                     NTEST = NTEST + 1
                  END IF
*
*                 Test (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
                     NTEST = NTEST + 1
                  END IF
*
                  NTESTT = NTESTT + NTEST
*
*                 Print out tests which fail.
*
                  DO 20 J = 1, 9
                     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