ddrgsx.f

famous linear algebra library (LAPACK) ports to windows

         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 = -17
      ELSE IF( LIWORK.LT.NSIZE+6 ) 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 = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2 )
*
*        workspace for sggesx
*
         MAXWRK = 9*( NSIZE+1 ) + NSIZE*
     $            ILAENV( 1, 'DGEQRF', ' ', NSIZE, 1, NSIZE, 0 )
         MAXWRK = MAX( MAXWRK, 9*( NSIZE+1 )+NSIZE*
     $            ILAENV( 1, 'DORGQR', ' ', NSIZE, 1, NSIZE, -1 ) )
*
*        workspace for dgesvd
*
         BDSPAC = 5*NSIZE*NSIZE / 2
         MAXWRK = MAX( MAXWRK, 3*NSIZE*NSIZE / 2+NSIZE*NSIZE*
     $            ILAENV( 1, 'DGEBRD', ' ', 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 = -19
*
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DDRGSX', -INFO )
         RETURN
      END IF
*
*     Important constants
*
      ULP = DLAMCH( 'P' )
      ULPINV = ONE / ULP
      SMLNUM = DLAMCH( 'S' ) / ULP
      BIGNUM = ONE / SMLNUM
      CALL DLABAD( 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 DGGESX, 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 DLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, AI,
     $                         LDA )
                  CALL DLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, BI,
     $                         LDA )
*
                  CALL DLATM5( 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 DLCTSX
*                 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 DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
                  CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
*
                  CALL DGGESX( 'V', 'V', 'S', DLCTSX, SENSE, 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 = 9999 )'DGGESX', LINFO, MPLUSN,
     $                  PRTYPE
                     INFO = LINFO
                     GO TO 30
                  END IF
*
*                 Compute the norm(A, B)
*
                  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 ) )
                  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.
                     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, PRTYPE
                              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, 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 causes, there are
*                    almost same number of test cases fail the test.
*
                     CALL DLAKF2( MM, MPLUSN-MM, AI, LDA,
     $                            AI( MM+1, MM+1 ), BI,
     $                            BI( MM+1, MM+1 ), C, LDC )
*
                     CALL DGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK,
     $                            1, WORK( 2 ), 1, WORK( 3 ), LWORK-2,
     $                            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
                           WRITE( NOUT, FMT = 9995 )'SGX'
*
*                          Matrix types
*
                           WRITE( NOUT, FMT = 9993 )
*
*                          Tests performed
*
                           WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
     $                        'transpose', ( '''', I = 1, 4 )
*
                        END IF