LAPACK: TESTING/EIG/dget39.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DGET39( RMAX, LMAX, NINFO, KNT )
00002 *
00003 *  -- LAPACK test routine (version 3.1) --
00004 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00005 *     November 2006
00006 *
00007 *     .. Scalar Arguments ..
00008       INTEGER            KNT, LMAX, NINFO
00009       DOUBLE PRECISION   RMAX
00010 *     ..
00011 *
00012 *  Purpose
00013 *  =======
00014 *
00015 *  DGET39 tests DLAQTR, a routine for solving the real or
00016 *  special complex quasi upper triangular system
00017 *
00018 *       op(T)*p = scale*c,
00019 *  or
00020 *       op(T + iB)*(p+iq) = scale*(c+id),
00021 *
00022 *  in real arithmetic. T is upper quasi-triangular.
00023 *  If it is complex, then the first diagonal block of T must be
00024 *  1 by 1, B has the special structure
00025 *
00026 *                 B = [ b(1) b(2) ... b(n) ]
00027 *                     [       w            ]
00028 *                     [           w        ]
00029 *                     [              .     ]
00030 *                     [                 w  ]
00031 *
00032 *  op(A) = A or A', where A' denotes the conjugate transpose of
00033 *  the matrix A.
00034 *
00035 *  On input, X = [ c ].  On output, X = [ p ].
00036 *                [ d ]                  [ q ]
00037 *
00038 *  Scale is an output less than or equal to 1, chosen to avoid
00039 *  overflow in X.
00040 *  This subroutine is specially designed for the condition number
00041 *  estimation in the eigenproblem routine DTRSNA.
00042 *
00043 *  The test code verifies that the following residual is order 1:
00044 *
00045 *       ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)||
00046 *     -----------------------------------------
00047 *         max(ulp*(||T||+||B||)*(||x1||+||x2||),
00048 *             (||T||+||B||)*smlnum/ulp,
00049 *             smlnum)
00050 *
00051 *  (The (||T||+||B||)*smlnum/ulp term accounts for possible
00052 *   (gradual or nongradual) underflow in x1 and x2.)
00053 *
00054 *  Arguments
00055 *  ==========
00056 *
00057 *  RMAX    (output) DOUBLE PRECISION
00058 *          Value of the largest test ratio.
00059 *
00060 *  LMAX    (output) INTEGER
00061 *          Example number where largest test ratio achieved.
00062 *
00063 *  NINFO   (output) INTEGER
00064 *          Number of examples where INFO is nonzero.
00065 *
00066 *  KNT     (output) INTEGER
00067 *          Total number of examples tested.
00068 *
00069 *  =====================================================================
00070 *
00071 *     .. Parameters ..
00072       INTEGER            LDT, LDT2
00073       PARAMETER          ( LDT = 10, LDT2 = 2*LDT )
00074       DOUBLE PRECISION   ZERO, ONE
00075       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00076 *     ..
00077 *     .. Local Scalars ..
00078       INTEGER            I, INFO, IVM1, IVM2, IVM3, IVM4, IVM5, J, K, N,
00079      $                   NDIM
00080       DOUBLE PRECISION   BIGNUM, DOMIN, DUMM, EPS, NORM, NORMTB, RESID,
00081      $                   SCALE, SMLNUM, W, XNORM
00082 *     ..
00083 *     .. External Functions ..
00084       INTEGER            IDAMAX
00085       DOUBLE PRECISION   DASUM, DDOT, DLAMCH, DLANGE
00086       EXTERNAL           IDAMAX, DASUM, DDOT, DLAMCH, DLANGE
00087 *     ..
00088 *     .. External Subroutines ..
00089       EXTERNAL           DCOPY, DGEMV, DLABAD, DLAQTR
00090 *     ..
00091 *     .. Intrinsic Functions ..
00092       INTRINSIC          ABS, COS, DBLE, MAX, SIN, SQRT
00093 *     ..
00094 *     .. Local Arrays ..
00095       INTEGER            IDIM( 6 ), IVAL( 5, 5, 6 )
00096       DOUBLE PRECISION   B( LDT ), D( LDT2 ), DUM( 1 ), T( LDT, LDT ),
00097      $                   VM1( 5 ), VM2( 5 ), VM3( 5 ), VM4( 5 ),
00098      $                   VM5( 3 ), WORK( LDT ), X( LDT2 ), Y( LDT2 )
00099 *     ..
00100 *     .. Data statements ..
00101       DATA               IDIM / 4, 5*5 /
00102       DATA               IVAL / 3, 4*0, 1, 1, -1, 0, 0, 3, 2, 1, 0, 0,
00103      $                   4, 3, 2, 2, 0, 5*0, 1, 4*0, 2, 2, 3*0, 3, 3, 4,
00104      $                   0, 0, 4, 2, 2, 3, 0, 4*1, 5, 1, 4*0, 2, 4, -2,
00105      $                   0, 0, 3, 3, 4, 0, 0, 4, 2, 2, 3, 0, 5*1, 1,
00106      $                   4*0, 2, 1, -1, 0, 0, 9, 8, 1, 0, 0, 4, 9, 1, 2,
00107      $                   -1, 5*2, 9, 4*0, 6, 4, 0, 0, 0, 3, 2, 1, 1, 0,
00108      $                   5, 1, -1, 1, 0, 5*2, 4, 4*0, 2, 2, 0, 0, 0, 1,
00109      $                   4, 4, 0, 0, 2, 4, 2, 2, -1, 5*2 /
00110 *     ..
00111 *     .. Executable Statements ..
00112 *
00113 *     Get machine parameters
00114 *
00115       EPS = DLAMCH( 'P' )
00116       SMLNUM = DLAMCH( 'S' )
00117       BIGNUM = ONE / SMLNUM
00118       CALL DLABAD( SMLNUM, BIGNUM )
00119 *
00120 *     Set up test case parameters
00121 *
00122       VM1( 1 ) = ONE
00123       VM1( 2 ) = SQRT( SMLNUM )
00124       VM1( 3 ) = SQRT( VM1( 2 ) )
00125       VM1( 4 ) = SQRT( BIGNUM )
00126       VM1( 5 ) = SQRT( VM1( 4 ) )
00127 *
00128       VM2( 1 ) = ONE
00129       VM2( 2 ) = SQRT( SMLNUM )
00130       VM2( 3 ) = SQRT( VM2( 2 ) )
00131       VM2( 4 ) = SQRT( BIGNUM )
00132       VM2( 5 ) = SQRT( VM2( 4 ) )
00133 *
00134       VM3( 1 ) = ONE
00135       VM3( 2 ) = SQRT( SMLNUM )
00136       VM3( 3 ) = SQRT( VM3( 2 ) )
00137       VM3( 4 ) = SQRT( BIGNUM )
00138       VM3( 5 ) = SQRT( VM3( 4 ) )
00139 *
00140       VM4( 1 ) = ONE
00141       VM4( 2 ) = SQRT( SMLNUM )
00142       VM4( 3 ) = SQRT( VM4( 2 ) )
00143       VM4( 4 ) = SQRT( BIGNUM )
00144       VM4( 5 ) = SQRT( VM4( 4 ) )
00145 *
00146       VM5( 1 ) = ONE
00147       VM5( 2 ) = EPS
00148       VM5( 3 ) = SQRT( SMLNUM )
00149 *
00150 *     Initalization
00151 *
00152       KNT = 0
00153       RMAX = ZERO
00154       NINFO = 0
00155       SMLNUM = SMLNUM / EPS
00156 *
00157 *     Begin test loop
00158 *
00159       DO 140 IVM5 = 1, 3
00160          DO 130 IVM4 = 1, 5
00161             DO 120 IVM3 = 1, 5
00162                DO 110 IVM2 = 1, 5
00163                   DO 100 IVM1 = 1, 5
00164                      DO 90 NDIM = 1, 6
00165 *
00166                         N = IDIM( NDIM )
00167                         DO 20 I = 1, N
00168                            DO 10 J = 1, N
00169                               T( I, J ) = DBLE( IVAL( I, J, NDIM ) )*
00170      $                                    VM1( IVM1 )
00171                               IF( I.GE.J )
00172      $                           T( I, J ) = T( I, J )*VM5( IVM5 )
00173    10                      CONTINUE
00174    20                   CONTINUE
00175 *
00176                         W = ONE*VM2( IVM2 )
00177 *
00178                         DO 30 I = 1, N
00179                            B( I ) = COS( DBLE( I ) )*VM3( IVM3 )
00180    30                   CONTINUE
00181 *
00182                         DO 40 I = 1, 2*N
00183                            D( I ) = SIN( DBLE( I ) )*VM4( IVM4 )
00184    40                   CONTINUE
00185 *
00186                         NORM = DLANGE( '1', N, N, T, LDT, WORK )
00187                         K = IDAMAX( N, B, 1 )
00188                         NORMTB = NORM + ABS( B( K ) ) + ABS( W )
00189 *
00190                         CALL DCOPY( N, D, 1, X, 1 )
00191                         KNT = KNT + 1
00192                         CALL DLAQTR( .FALSE., .TRUE., N, T, LDT, DUM,
00193      $                               DUMM, SCALE, X, WORK, INFO )
00194                         IF( INFO.NE.0 )
00195      $                     NINFO = NINFO + 1
00196 *
00197 *                       || T*x - scale*d || /
00198 *                         max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
00199 *
00200                         CALL DCOPY( N, D, 1, Y, 1 )
00201                         CALL DGEMV( 'No transpose', N, N, ONE, T, LDT,
00202      $                              X, 1, -SCALE, Y, 1 )
00203                         XNORM = DASUM( N, X, 1 )
00204                         RESID = DASUM( N, Y, 1 )
00205                         DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM,
00206      $                          ( NORM*EPS )*XNORM )
00207                         RESID = RESID / DOMIN
00208                         IF( RESID.GT.RMAX ) THEN
00209                            RMAX = RESID
00210                            LMAX = KNT
00211                         END IF
00212 *
00213                         CALL DCOPY( N, D, 1, X, 1 )
00214                         KNT = KNT + 1
00215                         CALL DLAQTR( .TRUE., .TRUE., N, T, LDT, DUM,
00216      $                               DUMM, SCALE, X, WORK, INFO )
00217                         IF( INFO.NE.0 )
00218      $                     NINFO = NINFO + 1
00219 *
00220 *                       || T*x - scale*d || /
00221 *                         max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
00222 *
00223                         CALL DCOPY( N, D, 1, Y, 1 )
00224                         CALL DGEMV( 'Transpose', N, N, ONE, T, LDT, X,
00225      $                              1, -SCALE, Y, 1 )
00226                         XNORM = DASUM( N, X, 1 )
00227                         RESID = DASUM( N, Y, 1 )
00228                         DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM,
00229      $                          ( NORM*EPS )*XNORM )
00230                         RESID = RESID / DOMIN
00231                         IF( RESID.GT.RMAX ) THEN
00232                            RMAX = RESID
00233                            LMAX = KNT
00234                         END IF
00235 *
00236                         CALL DCOPY( 2*N, D, 1, X, 1 )
00237                         KNT = KNT + 1
00238                         CALL DLAQTR( .FALSE., .FALSE., N, T, LDT, B, W,
00239      $                               SCALE, X, WORK, INFO )
00240                         IF( INFO.NE.0 )
00241      $                     NINFO = NINFO + 1
00242 *
00243 *                       ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
00244 *                          max(ulp*(||T||+||B||)*(||x1||+||x2||),
00245 *                                  smlnum/ulp * (||T||+||B||), smlnum )
00246 *
00247 *
00248                         CALL DCOPY( 2*N, D, 1, Y, 1 )
00249                         Y( 1 ) = DDOT( N, B, 1, X( 1+N ), 1 ) +
00250      $                           SCALE*Y( 1 )
00251                         DO 50 I = 2, N
00252                            Y( I ) = W*X( I+N ) + SCALE*Y( I )
00253    50                   CONTINUE
00254                         CALL DGEMV( 'No transpose', N, N, ONE, T, LDT,
00255      $                              X, 1, -ONE, Y, 1 )
00256 *
00257                         Y( 1+N ) = DDOT( N, B, 1, X, 1 ) -
00258      $                             SCALE*Y( 1+N )
00259                         DO 60 I = 2, N
00260                            Y( I+N ) = W*X( I ) - SCALE*Y( I+N )
00261    60                   CONTINUE
00262                         CALL DGEMV( 'No transpose', N, N, ONE, T, LDT,
00263      $                              X( 1+N ), 1, ONE, Y( 1+N ), 1 )
00264 *
00265                         RESID = DASUM( 2*N, Y, 1 )
00266                         DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB,
00267      $                          EPS*( NORMTB*DASUM( 2*N, X, 1 ) ) )
00268                         RESID = RESID / DOMIN
00269                         IF( RESID.GT.RMAX ) THEN
00270                            RMAX = RESID
00271                            LMAX = KNT
00272                         END IF
00273 *
00274                         CALL DCOPY( 2*N, D, 1, X, 1 )
00275                         KNT = KNT + 1
00276                         CALL DLAQTR( .TRUE., .FALSE., N, T, LDT, B, W,
00277      $                               SCALE, X, WORK, INFO )
00278                         IF( INFO.NE.0 )
00279      $                     NINFO = NINFO + 1
00280 *
00281 *                       ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
00282 *                          max(ulp*(||T||+||B||)*(||x1||+||x2||),
00283 *                                  smlnum/ulp * (||T||+||B||), smlnum )
00284 *
00285                         CALL DCOPY( 2*N, D, 1, Y, 1 )
00286                         Y( 1 ) = B( 1 )*X( 1+N ) - SCALE*Y( 1 )
00287                         DO 70 I = 2, N
00288                            Y( I ) = B( I )*X( 1+N ) + W*X( I+N ) -
00289      $                              SCALE*Y( I )
00290    70                   CONTINUE
00291                         CALL DGEMV( 'Transpose', N, N, ONE, T, LDT, X,
00292      $                              1, ONE, Y, 1 )
00293 *
00294                         Y( 1+N ) = B( 1 )*X( 1 ) + SCALE*Y( 1+N )
00295                         DO 80 I = 2, N
00296                            Y( I+N ) = B( I )*X( 1 ) + W*X( I ) +
00297      $                                SCALE*Y( I+N )
00298    80                   CONTINUE
00299                         CALL DGEMV( 'Transpose', N, N, ONE, T, LDT,
00300      $                              X( 1+N ), 1, -ONE, Y( 1+N ), 1 )
00301 *
00302                         RESID = DASUM( 2*N, Y, 1 )
00303                         DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB,
00304      $                          EPS*( NORMTB*DASUM( 2*N, X, 1 ) ) )
00305                         RESID = RESID / DOMIN
00306                         IF( RESID.GT.RMAX ) THEN
00307                            RMAX = RESID
00308                            LMAX = KNT
00309                         END IF
00310 *
00311    90                CONTINUE
00312   100             CONTINUE
00313   110          CONTINUE
00314   120       CONTINUE
00315   130    CONTINUE
00316   140 CONTINUE
00317 *
00318       RETURN
00319 *
00320 *     End of DGET39
00321 *
00322       END