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