LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dchkgk.f
Go to the documentation of this file.
00001 *> \brief \b DCHKGK
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 DCHKGK( NIN, NOUT )
00012 * 
00013 *       .. Scalar Arguments ..
00014 *       INTEGER            NIN, NOUT
00015 *       ..
00016 *  
00017 *
00018 *> \par Purpose:
00019 *  =============
00020 *>
00021 *> \verbatim
00022 *>
00023 *> DCHKGK tests DGGBAK, a routine for backward balancing  of
00024 *> a matrix pair (A, B).
00025 *> \endverbatim
00026 *
00027 *  Arguments:
00028 *  ==========
00029 *
00030 *> \param[in] NIN
00031 *> \verbatim
00032 *>          NIN is INTEGER
00033 *>          The logical unit number for input.  NIN > 0.
00034 *> \endverbatim
00035 *>
00036 *> \param[in] NOUT
00037 *> \verbatim
00038 *>          NOUT is INTEGER
00039 *>          The logical unit number for output.  NOUT > 0.
00040 *> \endverbatim
00041 *
00042 *  Authors:
00043 *  ========
00044 *
00045 *> \author Univ. of Tennessee 
00046 *> \author Univ. of California Berkeley 
00047 *> \author Univ. of Colorado Denver 
00048 *> \author NAG Ltd. 
00049 *
00050 *> \date November 2011
00051 *
00052 *> \ingroup double_eig
00053 *
00054 *  =====================================================================
00055       SUBROUTINE DCHKGK( NIN, NOUT )
00056 *
00057 *  -- LAPACK test routine (version 3.4.0) --
00058 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00059 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00060 *     November 2011
00061 *
00062 *     .. Scalar Arguments ..
00063       INTEGER            NIN, NOUT
00064 *     ..
00065 *
00066 *  =====================================================================
00067 *
00068 *     .. Parameters ..
00069       INTEGER            LDA, LDB, LDVL, LDVR
00070       PARAMETER          ( LDA = 50, LDB = 50, LDVL = 50, LDVR = 50 )
00071       INTEGER            LDE, LDF, LDWORK
00072       PARAMETER          ( LDE = 50, LDF = 50, LDWORK = 50 )
00073       DOUBLE PRECISION   ZERO, ONE
00074       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00075 *     ..
00076 *     .. Local Scalars ..
00077       INTEGER            I, IHI, ILO, INFO, J, KNT, M, N, NINFO
00078       DOUBLE PRECISION   ANORM, BNORM, EPS, RMAX, VMAX
00079 *     ..
00080 *     .. Local Arrays ..
00081       INTEGER            LMAX( 4 )
00082       DOUBLE PRECISION   A( LDA, LDA ), AF( LDA, LDA ), B( LDB, LDB ),
00083      $                   BF( LDB, LDB ), E( LDE, LDE ), F( LDF, LDF ),
00084      $                   LSCALE( LDA ), RSCALE( LDA ), VL( LDVL, LDVL ),
00085      $                   VLF( LDVL, LDVL ), VR( LDVR, LDVR ),
00086      $                   VRF( LDVR, LDVR ), WORK( LDWORK, LDWORK )
00087 *     ..
00088 *     .. External Functions ..
00089       DOUBLE PRECISION   DLAMCH, DLANGE
00090       EXTERNAL           DLAMCH, DLANGE
00091 *     ..
00092 *     .. External Subroutines ..
00093       EXTERNAL           DGEMM, DGGBAK, DGGBAL, DLACPY
00094 *     ..
00095 *     .. Intrinsic Functions ..
00096       INTRINSIC          ABS, MAX
00097 *     ..
00098 *     .. Executable Statements ..
00099 *
00100 *     Initialization
00101 *
00102       LMAX( 1 ) = 0
00103       LMAX( 2 ) = 0
00104       LMAX( 3 ) = 0
00105       LMAX( 4 ) = 0
00106       NINFO = 0
00107       KNT = 0
00108       RMAX = ZERO
00109 *
00110       EPS = DLAMCH( 'Precision' )
00111 *
00112    10 CONTINUE
00113       READ( NIN, FMT = * )N, M
00114       IF( N.EQ.0 )
00115      $   GO TO 100
00116 *
00117       DO 20 I = 1, N
00118          READ( NIN, FMT = * )( A( I, J ), J = 1, N )
00119    20 CONTINUE
00120 *
00121       DO 30 I = 1, N
00122          READ( NIN, FMT = * )( B( I, J ), J = 1, N )
00123    30 CONTINUE
00124 *
00125       DO 40 I = 1, N
00126          READ( NIN, FMT = * )( VL( I, J ), J = 1, M )
00127    40 CONTINUE
00128 *
00129       DO 50 I = 1, N
00130          READ( NIN, FMT = * )( VR( I, J ), J = 1, M )
00131    50 CONTINUE
00132 *
00133       KNT = KNT + 1
00134 *
00135       ANORM = DLANGE( 'M', N, N, A, LDA, WORK )
00136       BNORM = DLANGE( 'M', N, N, B, LDB, WORK )
00137 *
00138       CALL DLACPY( 'FULL', N, N, A, LDA, AF, LDA )
00139       CALL DLACPY( 'FULL', N, N, B, LDB, BF, LDB )
00140 *
00141       CALL DGGBAL( 'B', N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE,
00142      $             WORK, INFO )
00143       IF( INFO.NE.0 ) THEN
00144          NINFO = NINFO + 1
00145          LMAX( 1 ) = KNT
00146       END IF
00147 *
00148       CALL DLACPY( 'FULL', N, M, VL, LDVL, VLF, LDVL )
00149       CALL DLACPY( 'FULL', N, M, VR, LDVR, VRF, LDVR )
00150 *
00151       CALL DGGBAK( 'B', 'L', N, ILO, IHI, LSCALE, RSCALE, M, VL, LDVL,
00152      $             INFO )
00153       IF( INFO.NE.0 ) THEN
00154          NINFO = NINFO + 1
00155          LMAX( 2 ) = KNT
00156       END IF
00157 *
00158       CALL DGGBAK( 'B', 'R', N, ILO, IHI, LSCALE, RSCALE, M, VR, LDVR,
00159      $             INFO )
00160       IF( INFO.NE.0 ) THEN
00161          NINFO = NINFO + 1
00162          LMAX( 3 ) = KNT
00163       END IF
00164 *
00165 *     Test of DGGBAK
00166 *
00167 *     Check tilde(VL)'*A*tilde(VR) - VL'*tilde(A)*VR
00168 *     where tilde(A) denotes the transformed matrix.
00169 *
00170       CALL DGEMM( 'N', 'N', N, M, N, ONE, AF, LDA, VR, LDVR, ZERO, WORK,
00171      $            LDWORK )
00172       CALL DGEMM( 'T', 'N', M, M, N, ONE, VL, LDVL, WORK, LDWORK, ZERO,
00173      $            E, LDE )
00174 *
00175       CALL DGEMM( 'N', 'N', N, M, N, ONE, A, LDA, VRF, LDVR, ZERO, WORK,
00176      $            LDWORK )
00177       CALL DGEMM( 'T', 'N', M, M, N, ONE, VLF, LDVL, WORK, LDWORK, ZERO,
00178      $            F, LDF )
00179 *
00180       VMAX = ZERO
00181       DO 70 J = 1, M
00182          DO 60 I = 1, M
00183             VMAX = MAX( VMAX, ABS( E( I, J )-F( I, J ) ) )
00184    60    CONTINUE
00185    70 CONTINUE
00186       VMAX = VMAX / ( EPS*MAX( ANORM, BNORM ) )
00187       IF( VMAX.GT.RMAX ) THEN
00188          LMAX( 4 ) = KNT
00189          RMAX = VMAX
00190       END IF
00191 *
00192 *     Check tilde(VL)'*B*tilde(VR) - VL'*tilde(B)*VR
00193 *
00194       CALL DGEMM( 'N', 'N', N, M, N, ONE, BF, LDB, VR, LDVR, ZERO, WORK,
00195      $            LDWORK )
00196       CALL DGEMM( 'T', 'N', M, M, N, ONE, VL, LDVL, WORK, LDWORK, ZERO,
00197      $            E, LDE )
00198 *
00199       CALL DGEMM( 'N', 'N', N, M, N, ONE, B, LDB, VRF, LDVR, ZERO, WORK,
00200      $            LDWORK )
00201       CALL DGEMM( 'T', 'N', M, M, N, ONE, VLF, LDVL, WORK, LDWORK, ZERO,
00202      $            F, LDF )
00203 *
00204       VMAX = ZERO
00205       DO 90 J = 1, M
00206          DO 80 I = 1, M
00207             VMAX = MAX( VMAX, ABS( E( I, J )-F( I, J ) ) )
00208    80    CONTINUE
00209    90 CONTINUE
00210       VMAX = VMAX / ( EPS*MAX( ANORM, BNORM ) )
00211       IF( VMAX.GT.RMAX ) THEN
00212          LMAX( 4 ) = KNT
00213          RMAX = VMAX
00214       END IF
00215 *
00216       GO TO 10
00217 *
00218   100 CONTINUE
00219 *
00220       WRITE( NOUT, FMT = 9999 )
00221  9999 FORMAT( 1X, '.. test output of DGGBAK .. ' )
00222 *
00223       WRITE( NOUT, FMT = 9998 )RMAX
00224  9998 FORMAT( ' value of largest test error                  =', D12.3 )
00225       WRITE( NOUT, FMT = 9997 )LMAX( 1 )
00226  9997 FORMAT( ' example number where DGGBAL info is not 0    =', I4 )
00227       WRITE( NOUT, FMT = 9996 )LMAX( 2 )
00228  9996 FORMAT( ' example number where DGGBAK(L) info is not 0 =', I4 )
00229       WRITE( NOUT, FMT = 9995 )LMAX( 3 )
00230  9995 FORMAT( ' example number where DGGBAK(R) info is not 0 =', I4 )
00231       WRITE( NOUT, FMT = 9994 )LMAX( 4 )
00232  9994 FORMAT( ' example number having largest error          =', I4 )
00233       WRITE( NOUT, FMT = 9993 )NINFO
00234  9993 FORMAT( ' number of examples where info is not 0       =', I4 )
00235       WRITE( NOUT, FMT = 9992 )KNT
00236  9992 FORMAT( ' total number of examples tested              =', I4 )
00237 *
00238       RETURN
00239 *
00240 *     End of DCHKGK
00241 *
00242       END
 All Files Functions