54
55
56
57
58
59
60 INTEGER NIN, NOUT
61
62
63
64
65
66 INTEGER LDA, LDB, LDVL, LDVR
67 parameter( lda = 50, ldb = 50, ldvl = 50, ldvr = 50 )
68 INTEGER LDE, LDF, LDWORK
69 parameter(
lde = 50, ldf = 50, ldwork = 50 )
70 DOUBLE PRECISION ZERO, ONE
71 parameter( zero = 0.0d+0, one = 1.0d+0 )
72
73
74 INTEGER I, IHI, ILO, INFO, J, KNT, M, N, NINFO
75 DOUBLE PRECISION ANORM, BNORM, EPS, RMAX, VMAX
76
77
78 INTEGER LMAX( 4 )
79 DOUBLE PRECISION A( LDA, LDA ), AF( LDA, LDA ), B( LDB, LDB ),
80 $ BF( LDB, LDB ), E( LDE, LDE ), F( LDF, LDF ),
81 $ LSCALE( LDA ), RSCALE( LDA ), VL( LDVL, LDVL ),
82 $ VLF( LDVL, LDVL ), VR( LDVR, LDVR ),
83 $ VRF( LDVR, LDVR ), WORK( LDWORK, LDWORK )
84
85
86 DOUBLE PRECISION DLAMCH, DLANGE
88
89
91
92
93 INTRINSIC abs, max
94
95
96
97
98
99 lmax( 1 ) = 0
100 lmax( 2 ) = 0
101 lmax( 3 ) = 0
102 lmax( 4 ) = 0
103 ninfo = 0
104 knt = 0
105 rmax = zero
106
107 eps =
dlamch(
'Precision' )
108
109 10 CONTINUE
110 READ( nin, fmt = * )n, m
111 IF( n.EQ.0 )
112 $ GO TO 100
113
114 DO 20 i = 1, n
115 READ( nin, fmt = * )( a( i, j ), j = 1, n )
116 20 CONTINUE
117
118 DO 30 i = 1, n
119 READ( nin, fmt = * )( b( i, j ), j = 1, n )
120 30 CONTINUE
121
122 DO 40 i = 1, n
123 READ( nin, fmt = * )( vl( i, j ), j = 1, m )
124 40 CONTINUE
125
126 DO 50 i = 1, n
127 READ( nin, fmt = * )( vr( i, j ), j = 1, m )
128 50 CONTINUE
129
130 knt = knt + 1
131
132 anorm =
dlange(
'M', n, n, a, lda, work )
133 bnorm =
dlange(
'M', n, n, b, ldb, work )
134
135 CALL dlacpy(
'FULL', n, n, a, lda, af, lda )
136 CALL dlacpy(
'FULL', n, n, b, ldb, bf, ldb )
137
138 CALL dggbal(
'B', n, a, lda, b, ldb, ilo, ihi, lscale, rscale,
139 $ work, info )
140 IF( info.NE.0 ) THEN
141 ninfo = ninfo + 1
142 lmax( 1 ) = knt
143 END IF
144
145 CALL dlacpy(
'FULL', n, m, vl, ldvl, vlf, ldvl )
146 CALL dlacpy(
'FULL', n, m, vr, ldvr, vrf, ldvr )
147
148 CALL dggbak(
'B',
'L', n, ilo, ihi, lscale, rscale, m, vl, ldvl,
149 $ info )
150 IF( info.NE.0 ) THEN
151 ninfo = ninfo + 1
152 lmax( 2 ) = knt
153 END IF
154
155 CALL dggbak(
'B',
'R', n, ilo, ihi, lscale, rscale, m, vr, ldvr,
156 $ info )
157 IF( info.NE.0 ) THEN
158 ninfo = ninfo + 1
159 lmax( 3 ) = knt
160 END IF
161
162
163
164
165
166
167 CALL dgemm(
'N',
'N', n, m, n, one, af, lda, vr, ldvr, zero, work,
168 $ ldwork )
169 CALL dgemm(
'T',
'N', m, m, n, one, vl, ldvl, work, ldwork, zero,
171
172 CALL dgemm(
'N',
'N', n, m, n, one, a, lda, vrf, ldvr, zero, work,
173 $ ldwork )
174 CALL dgemm(
'T',
'N', m, m, n, one, vlf, ldvl, work, ldwork, zero,
175 $ f, ldf )
176
177 vmax = zero
178 DO 70 j = 1, m
179 DO 60 i = 1, m
180 vmax = max( vmax, abs( e( i, j )-f( i, j ) ) )
181 60 CONTINUE
182 70 CONTINUE
183 vmax = vmax / ( eps*max( anorm, bnorm ) )
184 IF( vmax.GT.rmax ) THEN
185 lmax( 4 ) = knt
186 rmax = vmax
187 END IF
188
189
190
191 CALL dgemm(
'N',
'N', n, m, n, one, bf, ldb, vr, ldvr, zero, work,
192 $ ldwork )
193 CALL dgemm(
'T',
'N', m, m, n, one, vl, ldvl, work, ldwork, zero,
195
196 CALL dgemm(
'N',
'N', n, m, n, one, b, ldb, vrf, ldvr, zero, work,
197 $ ldwork )
198 CALL dgemm(
'T',
'N', m, m, n, one, vlf, ldvl, work, ldwork, zero,
199 $ f, ldf )
200
201 vmax = zero
202 DO 90 j = 1, m
203 DO 80 i = 1, m
204 vmax = max( vmax, abs( e( i, j )-f( i, j ) ) )
205 80 CONTINUE
206 90 CONTINUE
207 vmax = vmax / ( eps*max( anorm, bnorm ) )
208 IF( vmax.GT.rmax ) THEN
209 lmax( 4 ) = knt
210 rmax = vmax
211 END IF
212
213 GO TO 10
214
215 100 CONTINUE
216
217 WRITE( nout, fmt = 9999 )
218 9999 FORMAT( 1x, '.. test output of DGGBAK .. ' )
219
220 WRITE( nout, fmt = 9998 )rmax
221 9998 FORMAT( ' value of largest test error =', d12.3 )
222 WRITE( nout, fmt = 9997 )lmax( 1 )
223 9997 FORMAT( ' example number where DGGBAL info is not 0 =', i4 )
224 WRITE( nout, fmt = 9996 )lmax( 2 )
225 9996 FORMAT( ' example number where DGGBAK(L) info is not 0 =', i4 )
226 WRITE( nout, fmt = 9995 )lmax( 3 )
227 9995 FORMAT( ' example number where DGGBAK(R) info is not 0 =', i4 )
228 WRITE( nout, fmt = 9994 )lmax( 4 )
229 9994 FORMAT( ' example number having largest error =', i4 )
230 WRITE( nout, fmt = 9993 )ninfo
231 9993 FORMAT( ' number of examples where info is not 0 =', i4 )
232 WRITE( nout, fmt = 9992 )knt
233 9992 FORMAT( ' total number of examples tested =', i4 )
234
235 RETURN
236
237
238
logical function lde(ri, rj, lr)
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
subroutine dggbak(job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
DGGBAK
subroutine dggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
DGGBAL
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
double precision function dlamch(cmach)
DLAMCH
double precision function dlange(norm, m, n, a, lda, work)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...