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 REAL ZERO, ONE
71 parameter( zero = 0.0e+0, one = 1.0e+0 )
72
73
74 INTEGER I, IHI, ILO, INFO, J, KNT, M, N, NINFO
75 REAL ANORM, BNORM, EPS, RMAX, VMAX
76
77
78 INTEGER LMAX( 4 )
79 REAL 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 REAL SLAMCH, SLANGE
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 =
slamch(
'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 =
slange(
'M', n, n, a, lda, work )
133 bnorm =
slange(
'M', n, n, b, ldb, work )
134
135 CALL slacpy(
'FULL', n, n, a, lda, af, lda )
136 CALL slacpy(
'FULL', n, n, b, ldb, bf, ldb )
137
138 CALL sggbal(
'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 slacpy(
'FULL', n, m, vl, ldvl, vlf, ldvl )
146 CALL slacpy(
'FULL', n, m, vr, ldvr, vrf, ldvr )
147
148 CALL sggbak(
'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 sggbak(
'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 sgemm(
'N',
'N', n, m, n, one, af, lda, vr, ldvr, zero, work,
168 $ ldwork )
169 CALL sgemm(
'T',
'N', m, m, n, one, vl, ldvl, work, ldwork, zero,
171
172 CALL sgemm(
'N',
'N', n, m, n, one, a, lda, vrf, ldvr, zero, work,
173 $ ldwork )
174 CALL sgemm(
'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 sgemm(
'N',
'N', n, m, n, one, bf, ldb, vr, ldvr, zero, work,
192 $ ldwork )
193 CALL sgemm(
'T',
'N', m, m, n, one, vl, ldvl, work, ldwork, zero,
195
196 CALL sgemm(
'N',
'N', n, m, n, one, b, ldb, vrf, ldvr, zero, work,
197 $ ldwork )
198 CALL sgemm(
'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 SGGBAK .. ' )
219
220 WRITE( nout, fmt = 9998 )rmax
221 9998 FORMAT( ' value of largest test error =', e12.3 )
222 WRITE( nout, fmt = 9997 )lmax( 1 )
223 9997 FORMAT( ' example number where SGGBAL info is not 0 =', i4 )
224 WRITE( nout, fmt = 9996 )lmax( 2 )
225 9996 FORMAT( ' example number where SGGBAK(L) info is not 0 =', i4 )
226 WRITE( nout, fmt = 9995 )lmax( 3 )
227 9995 FORMAT( ' example number where SGGBAK(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 = 9992 )ninfo
231 9992 FORMAT( ' number of examples where info is not 0 =', i4 )
232 WRITE( nout, fmt = 9991 )knt
233 9991 FORMAT( ' total number of examples tested =', i4 )
234
235 RETURN
236
237
238
logical function lde(ri, rj, lr)
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
subroutine sggbak(job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
SGGBAK
subroutine sggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
SGGBAL
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
real function slamch(cmach)
SLAMCH
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...