LAPACK: TESTING/EIG/schkgk.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SCHKGK( NIN, NOUT )
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            NIN, NOUT
00009 *     ..
00010 *
00011 *  Purpose
00012 *  =======
00013 *
00014 *  SCHKGK tests SGGBAK, a routine for backward balancing  of
00015 *  a matrix pair (A, B).
00016 *
00017 *  Arguments
00018 *  =========
00019 *
00020 *  NIN     (input) INTEGER
00021 *          The logical unit number for input.  NIN > 0.
00022 *
00023 *  NOUT    (input) INTEGER
00024 *          The logical unit number for output.  NOUT > 0.
00025 *
00026 *  =====================================================================
00027 *
00028 *     .. Parameters ..
00029       INTEGER            LDA, LDB, LDVL, LDVR
00030       PARAMETER          ( LDA = 50, LDB = 50, LDVL = 50, LDVR = 50 )
00031       INTEGER            LDE, LDF, LDWORK
00032       PARAMETER          ( LDE = 50, LDF = 50, LDWORK = 50 )
00033       REAL               ZERO, ONE
00034       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00035 *     ..
00036 *     .. Local Scalars ..
00037       INTEGER            I, IHI, ILO, INFO, J, KNT, M, N, NINFO
00038       REAL               ANORM, BNORM, EPS, RMAX, VMAX
00039 *     ..
00040 *     .. Local Arrays ..
00041       INTEGER            LMAX( 4 )
00042       REAL               A( LDA, LDA ), AF( LDA, LDA ), B( LDB, LDB ),
00043      $                   BF( LDB, LDB ), E( LDE, LDE ), F( LDF, LDF ),
00044      $                   LSCALE( LDA ), RSCALE( LDA ), VL( LDVL, LDVL ),
00045      $                   VLF( LDVL, LDVL ), VR( LDVR, LDVR ),
00046      $                   VRF( LDVR, LDVR ), WORK( LDWORK, LDWORK )
00047 *     ..
00048 *     .. External Functions ..
00049       REAL               SLAMCH, SLANGE
00050       EXTERNAL           SLAMCH, SLANGE
00051 *     ..
00052 *     .. External Subroutines ..
00053       EXTERNAL           SGEMM, SGGBAK, SGGBAL, SLACPY
00054 *     ..
00055 *     .. Intrinsic Functions ..
00056       INTRINSIC          ABS, MAX
00057 *     ..
00058 *     .. Executable Statements ..
00059 *
00060 *     Initialization
00061 *
00062       LMAX( 1 ) = 0
00063       LMAX( 2 ) = 0
00064       LMAX( 3 ) = 0
00065       LMAX( 4 ) = 0
00066       NINFO = 0
00067       KNT = 0
00068       RMAX = ZERO
00069 *
00070       EPS = SLAMCH( 'Precision' )
00071 *
00072    10 CONTINUE
00073       READ( NIN, FMT = * )N, M
00074       IF( N.EQ.0 )
00075      $   GO TO 100
00076 *
00077       DO 20 I = 1, N
00078          READ( NIN, FMT = * )( A( I, J ), J = 1, N )
00079    20 CONTINUE
00080 *
00081       DO 30 I = 1, N
00082          READ( NIN, FMT = * )( B( I, J ), J = 1, N )
00083    30 CONTINUE
00084 *
00085       DO 40 I = 1, N
00086          READ( NIN, FMT = * )( VL( I, J ), J = 1, M )
00087    40 CONTINUE
00088 *
00089       DO 50 I = 1, N
00090          READ( NIN, FMT = * )( VR( I, J ), J = 1, M )
00091    50 CONTINUE
00092 *
00093       KNT = KNT + 1
00094 *
00095       ANORM = SLANGE( 'M', N, N, A, LDA, WORK )
00096       BNORM = SLANGE( 'M', N, N, B, LDB, WORK )
00097 *
00098       CALL SLACPY( 'FULL', N, N, A, LDA, AF, LDA )
00099       CALL SLACPY( 'FULL', N, N, B, LDB, BF, LDB )
00100 *
00101       CALL SGGBAL( 'B', N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE,
00102      $             WORK, INFO )
00103       IF( INFO.NE.0 ) THEN
00104          NINFO = NINFO + 1
00105          LMAX( 1 ) = KNT
00106       END IF
00107 *
00108       CALL SLACPY( 'FULL', N, M, VL, LDVL, VLF, LDVL )
00109       CALL SLACPY( 'FULL', N, M, VR, LDVR, VRF, LDVR )
00110 *
00111       CALL SGGBAK( 'B', 'L', N, ILO, IHI, LSCALE, RSCALE, M, VL, LDVL,
00112      $             INFO )
00113       IF( INFO.NE.0 ) THEN
00114          NINFO = NINFO + 1
00115          LMAX( 2 ) = KNT
00116       END IF
00117 *
00118       CALL SGGBAK( 'B', 'R', N, ILO, IHI, LSCALE, RSCALE, M, VR, LDVR,
00119      $             INFO )
00120       IF( INFO.NE.0 ) THEN
00121          NINFO = NINFO + 1
00122          LMAX( 3 ) = KNT
00123       END IF
00124 *
00125 *     Test of SGGBAK
00126 *
00127 *     Check tilde(VL)'*A*tilde(VR) - VL'*tilde(A)*VR
00128 *     where tilde(A) denotes the transformed matrix.
00129 *
00130       CALL SGEMM( 'N', 'N', N, M, N, ONE, AF, LDA, VR, LDVR, ZERO, WORK,
00131      $            LDWORK )
00132       CALL SGEMM( 'T', 'N', M, M, N, ONE, VL, LDVL, WORK, LDWORK, ZERO,
00133      $            E, LDE )
00134 *
00135       CALL SGEMM( 'N', 'N', N, M, N, ONE, A, LDA, VRF, LDVR, ZERO, WORK,
00136      $            LDWORK )
00137       CALL SGEMM( 'T', 'N', M, M, N, ONE, VLF, LDVL, WORK, LDWORK, ZERO,
00138      $            F, LDF )
00139 *
00140       VMAX = ZERO
00141       DO 70 J = 1, M
00142          DO 60 I = 1, M
00143             VMAX = MAX( VMAX, ABS( E( I, J )-F( I, J ) ) )
00144    60    CONTINUE
00145    70 CONTINUE
00146       VMAX = VMAX / ( EPS*MAX( ANORM, BNORM ) )
00147       IF( VMAX.GT.RMAX ) THEN
00148          LMAX( 4 ) = KNT
00149          RMAX = VMAX
00150       END IF
00151 *
00152 *     Check tilde(VL)'*B*tilde(VR) - VL'*tilde(B)*VR
00153 *
00154       CALL SGEMM( 'N', 'N', N, M, N, ONE, BF, LDB, VR, LDVR, ZERO, WORK,
00155      $            LDWORK )
00156       CALL SGEMM( 'T', 'N', M, M, N, ONE, VL, LDVL, WORK, LDWORK, ZERO,
00157      $            E, LDE )
00158 *
00159       CALL SGEMM( 'N', 'N', N, M, N, ONE, B, LDB, VRF, LDVR, ZERO, WORK,
00160      $            LDWORK )
00161       CALL SGEMM( 'T', 'N', M, M, N, ONE, VLF, LDVL, WORK, LDWORK, ZERO,
00162      $            F, LDF )
00163 *
00164       VMAX = ZERO
00165       DO 90 J = 1, M
00166          DO 80 I = 1, M
00167             VMAX = MAX( VMAX, ABS( E( I, J )-F( I, J ) ) )
00168    80    CONTINUE
00169    90 CONTINUE
00170       VMAX = VMAX / ( EPS*MAX( ANORM, BNORM ) )
00171       IF( VMAX.GT.RMAX ) THEN
00172          LMAX( 4 ) = KNT
00173          RMAX = VMAX
00174       END IF
00175 *
00176       GO TO 10
00177 *
00178   100 CONTINUE
00179 *
00180       WRITE( NOUT, FMT = 9999 )
00181  9999 FORMAT( 1X, '.. test output of SGGBAK .. ' )
00182 *
00183       WRITE( NOUT, FMT = 9998 )RMAX
00184  9998 FORMAT( ' value of largest test error                  =', E12.3 )
00185       WRITE( NOUT, FMT = 9997 )LMAX( 1 )
00186  9997 FORMAT( ' example number where SGGBAL info is not 0    =', I4 )
00187       WRITE( NOUT, FMT = 9996 )LMAX( 2 )
00188  9996 FORMAT( ' example number where SGGBAK(L) info is not 0 =', I4 )
00189       WRITE( NOUT, FMT = 9995 )LMAX( 3 )
00190  9995 FORMAT( ' example number where SGGBAK(R) info is not 0 =', I4 )
00191       WRITE( NOUT, FMT = 9994 )LMAX( 4 )
00192  9994 FORMAT( ' example number having largest error          =', I4 )
00193       WRITE( NOUT, FMT = 9992 )NINFO
00194  9992 FORMAT( ' number of examples where info is not 0       =', I4 )
00195       WRITE( NOUT, FMT = 9991 )KNT
00196  9991 FORMAT( ' total number of examples tested              =', I4 )
00197 *
00198       RETURN
00199 *
00200 *     End of SCHKGK
00201 *
00202       END