54
55
56
57
58
59
60 INTEGER NOUT
61 DOUBLE PRECISION THRESH
62
63
64
65
66
67 DOUBLE PRECISION ZERO, ONE, TEN
68 parameter( zero = 0.0d0, one = 1.0d+0, ten = 1.0d1 )
69 INTEGER NSZ, NSZB
70 parameter( nsz = 5, nszb = 3*nsz-2 )
71 INTEGER NSZP, NPOW
72 parameter( nszp = ( nsz*( nsz+1 ) ) / 2,
73 $ npow = 2*nsz+1 )
74
75
76 LOGICAL OK
77 CHARACTER*3 PATH
78 INTEGER I, INFO, J, KL, KU, M, N
79 DOUBLE PRECISION CCOND, EPS, NORM, RATIO, RCMAX, RCMIN, RCOND
80
81
82 DOUBLE PRECISION A( NSZ, NSZ ), AB( NSZB, NSZ ), AP( NSZP ),
83 $ C( NSZ ), POW( NPOW ), R( NSZ ), RESLTS( 5 ),
84 $ RPOW( NPOW )
85
86
87 DOUBLE PRECISION DLAMCH
89
90
92
93
94 INTRINSIC abs, max, min
95
96
97
98 path( 1: 1 ) = 'Double precision'
99 path( 2: 3 ) = 'EQ'
100
102 DO 10 i = 1, 5
103 reslts( i ) = zero
104 10 CONTINUE
105 DO 20 i = 1, npow
106 pow( i ) = ten**( i-1 )
107 rpow( i ) = one / pow( i )
108 20 CONTINUE
109
110
111
112 DO 80 n = 0, nsz
113 DO 70 m = 0, nsz
114
115 DO 40 j = 1, nsz
116 DO 30 i = 1, nsz
117 IF( i.LE.m .AND. j.LE.n ) THEN
118 a( i, j ) = pow( i+j+1 )*( -1 )**( i+j )
119 ELSE
120 a( i, j ) = zero
121 END IF
122 30 CONTINUE
123 40 CONTINUE
124
125 CALL dgeequ( m, n, a, nsz, r, c, rcond, ccond, norm, info )
126
127 IF( info.NE.0 ) THEN
128 reslts( 1 ) = one
129 ELSE
130 IF( n.NE.0 .AND. m.NE.0 ) THEN
131 reslts( 1 ) = max( reslts( 1 ),
132 $ abs( ( rcond-rpow( m ) ) / rpow( m ) ) )
133 reslts( 1 ) = max( reslts( 1 ),
134 $ abs( ( ccond-rpow( n ) ) / rpow( n ) ) )
135 reslts( 1 ) = max( reslts( 1 ),
136 $ abs( ( norm-pow( n+m+1 ) ) / pow( n+m+
137 $ 1 ) ) )
138 DO 50 i = 1, m
139 reslts( 1 ) = max( reslts( 1 ),
140 $ abs( ( r( i )-rpow( i+n+1 ) ) /
141 $ rpow( i+n+1 ) ) )
142 50 CONTINUE
143 DO 60 j = 1, n
144 reslts( 1 ) = max( reslts( 1 ),
145 $ abs( ( c( j )-pow( n-j+1 ) ) /
146 $ pow( n-j+1 ) ) )
147 60 CONTINUE
148 END IF
149 END IF
150
151 70 CONTINUE
152 80 CONTINUE
153
154
155
156 DO 90 j = 1, nsz
157 a( max( nsz-1, 1 ), j ) = zero
158 90 CONTINUE
159 CALL dgeequ( nsz, nsz, a, nsz, r, c, rcond, ccond, norm, info )
160 IF( info.NE.max( nsz-1, 1 ) )
161 $ reslts( 1 ) = one
162
163 DO 100 j = 1, nsz
164 a( max( nsz-1, 1 ), j ) = one
165 100 CONTINUE
166 DO 110 i = 1, nsz
167 a( i, max( nsz-1, 1 ) ) = zero
168 110 CONTINUE
169 CALL dgeequ( nsz, nsz, a, nsz, r, c, rcond, ccond, norm, info )
170 IF( info.NE.nsz+max( nsz-1, 1 ) )
171 $ reslts( 1 ) = one
172 reslts( 1 ) = reslts( 1 ) / eps
173
174
175
176 DO 250 n = 0, nsz
177 DO 240 m = 0, nsz
178 DO 230 kl = 0, max( m-1, 0 )
179 DO 220 ku = 0, max( n-1, 0 )
180
181 DO 130 j = 1, nsz
182 DO 120 i = 1, nszb
183 ab( i, j ) = zero
184 120 CONTINUE
185 130 CONTINUE
186 DO 150 j = 1, n
187 DO 140 i = 1, m
188 IF( i.LE.min( m, j+kl ) .AND. i.GE.
189 $ max( 1, j-ku ) .AND. j.LE.n ) THEN
190 ab( ku+1+i-j, j ) = pow( i+j+1 )*
191 $ ( -1 )**( i+j )
192 END IF
193 140 CONTINUE
194 150 CONTINUE
195
196 CALL dgbequ( m, n, kl, ku, ab, nszb, r, c, rcond,
197 $ ccond, norm, info )
198
199 IF( info.NE.0 ) THEN
200 IF( .NOT.( ( n+kl.LT.m .AND. info.EQ.n+kl+1 ) .OR.
201 $ ( m+ku.LT.n .AND. info.EQ.2*m+ku+1 ) ) ) THEN
202 reslts( 2 ) = one
203 END IF
204 ELSE
205 IF( n.NE.0 .AND. m.NE.0 ) THEN
206
207 rcmin = r( 1 )
208 rcmax = r( 1 )
209 DO 160 i = 1, m
210 rcmin = min( rcmin, r( i ) )
211 rcmax = max( rcmax, r( i ) )
212 160 CONTINUE
213 ratio = rcmin / rcmax
214 reslts( 2 ) = max( reslts( 2 ),
215 $ abs( ( rcond-ratio ) / ratio ) )
216
217 rcmin = c( 1 )
218 rcmax = c( 1 )
219 DO 170 j = 1, n
220 rcmin = min( rcmin, c( j ) )
221 rcmax = max( rcmax, c( j ) )
222 170 CONTINUE
223 ratio = rcmin / rcmax
224 reslts( 2 ) = max( reslts( 2 ),
225 $ abs( ( ccond-ratio ) / ratio ) )
226
227 reslts( 2 ) = max( reslts( 2 ),
228 $ abs( ( norm-pow( n+m+1 ) ) /
229 $ pow( n+m+1 ) ) )
230 DO 190 i = 1, m
231 rcmax = zero
232 DO 180 j = 1, n
233 IF( i.LE.j+kl .AND. i.GE.j-ku ) THEN
234 ratio = abs( r( i )*pow( i+j+1 )*
235 $ c( j ) )
236 rcmax = max( rcmax, ratio )
237 END IF
238 180 CONTINUE
239 reslts( 2 ) = max( reslts( 2 ),
240 $ abs( one-rcmax ) )
241 190 CONTINUE
242
243 DO 210 j = 1, n
244 rcmax = zero
245 DO 200 i = 1, m
246 IF( i.LE.j+kl .AND. i.GE.j-ku ) THEN
247 ratio = abs( r( i )*pow( i+j+1 )*
248 $ c( j ) )
249 rcmax = max( rcmax, ratio )
250 END IF
251 200 CONTINUE
252 reslts( 2 ) = max( reslts( 2 ),
253 $ abs( one-rcmax ) )
254 210 CONTINUE
255 END IF
256 END IF
257
258 220 CONTINUE
259 230 CONTINUE
260 240 CONTINUE
261 250 CONTINUE
262 reslts( 2 ) = reslts( 2 ) / eps
263
264
265
266 DO 290 n = 0, nsz
267
268 DO 270 i = 1, nsz
269 DO 260 j = 1, nsz
270 IF( i.LE.n .AND. j.EQ.i ) THEN
271 a( i, j ) = pow( i+j+1 )*( -1 )**( i+j )
272 ELSE
273 a( i, j ) = zero
274 END IF
275 260 CONTINUE
276 270 CONTINUE
277
278 CALL dpoequ( n, a, nsz, r, rcond, norm, info )
279
280 IF( info.NE.0 ) THEN
281 reslts( 3 ) = one
282 ELSE
283 IF( n.NE.0 ) THEN
284 reslts( 3 ) = max( reslts( 3 ),
285 $ abs( ( rcond-rpow( n ) ) / rpow( n ) ) )
286 reslts( 3 ) = max( reslts( 3 ),
287 $ abs( ( norm-pow( 2*n+1 ) ) / pow( 2*n+
288 $ 1 ) ) )
289 DO 280 i = 1, n
290 reslts( 3 ) = max( reslts( 3 ),
291 $ abs( ( r( i )-rpow( i+1 ) ) / rpow( i+
292 $ 1 ) ) )
293 280 CONTINUE
294 END IF
295 END IF
296 290 CONTINUE
297 a( max( nsz-1, 1 ), max( nsz-1, 1 ) ) = -one
298 CALL dpoequ( nsz, a, nsz, r, rcond, norm, info )
299 IF( info.NE.max( nsz-1, 1 ) )
300 $ reslts( 3 ) = one
301 reslts( 3 ) = reslts( 3 ) / eps
302
303
304
305 DO 360 n = 0, nsz
306
307
308
309 DO 300 i = 1, ( n*( n+1 ) ) / 2
310 ap( i ) = zero
311 300 CONTINUE
312 DO 310 i = 1, n
313 ap( ( i*( i+1 ) ) / 2 ) = pow( 2*i+1 )
314 310 CONTINUE
315
316 CALL dppequ(
'U', n, ap, r, rcond, norm, info )
317
318 IF( info.NE.0 ) THEN
319 reslts( 4 ) = one
320 ELSE
321 IF( n.NE.0 ) THEN
322 reslts( 4 ) = max( reslts( 4 ),
323 $ abs( ( rcond-rpow( n ) ) / rpow( n ) ) )
324 reslts( 4 ) = max( reslts( 4 ),
325 $ abs( ( norm-pow( 2*n+1 ) ) / pow( 2*n+
326 $ 1 ) ) )
327 DO 320 i = 1, n
328 reslts( 4 ) = max( reslts( 4 ),
329 $ abs( ( r( i )-rpow( i+1 ) ) / rpow( i+
330 $ 1 ) ) )
331 320 CONTINUE
332 END IF
333 END IF
334
335
336
337 DO 330 i = 1, ( n*( n+1 ) ) / 2
338 ap( i ) = zero
339 330 CONTINUE
340 j = 1
341 DO 340 i = 1, n
342 ap( j ) = pow( 2*i+1 )
343 j = j + ( n-i+1 )
344 340 CONTINUE
345
346 CALL dppequ(
'L', n, ap, r, rcond, norm, info )
347
348 IF( info.NE.0 ) THEN
349 reslts( 4 ) = one
350 ELSE
351 IF( n.NE.0 ) THEN
352 reslts( 4 ) = max( reslts( 4 ),
353 $ abs( ( rcond-rpow( n ) ) / rpow( n ) ) )
354 reslts( 4 ) = max( reslts( 4 ),
355 $ abs( ( norm-pow( 2*n+1 ) ) / pow( 2*n+
356 $ 1 ) ) )
357 DO 350 i = 1, n
358 reslts( 4 ) = max( reslts( 4 ),
359 $ abs( ( r( i )-rpow( i+1 ) ) / rpow( i+
360 $ 1 ) ) )
361 350 CONTINUE
362 END IF
363 END IF
364
365 360 CONTINUE
366 i = ( nsz*( nsz+1 ) ) / 2 - 2
367 ap( i ) = -one
368 CALL dppequ(
'L', nsz, ap, r, rcond, norm, info )
369 IF( info.NE.max( nsz-1, 1 ) )
370 $ reslts( 4 ) = one
371 reslts( 4 ) = reslts( 4 ) / eps
372
373
374
375 DO 460 n = 0, nsz
376 DO 450 kl = 0, max( n-1, 0 )
377
378
379
380 DO 380 j = 1, nsz
381 DO 370 i = 1, nszb
382 ab( i, j ) = zero
383 370 CONTINUE
384 380 CONTINUE
385 DO 390 j = 1, n
386 ab( kl+1, j ) = pow( 2*j+1 )
387 390 CONTINUE
388
389 CALL dpbequ(
'U', n, kl, ab, nszb, r, rcond, norm, info )
390
391 IF( info.NE.0 ) THEN
392 reslts( 5 ) = one
393 ELSE
394 IF( n.NE.0 ) THEN
395 reslts( 5 ) = max( reslts( 5 ),
396 $ abs( ( rcond-rpow( n ) ) / rpow( n ) ) )
397 reslts( 5 ) = max( reslts( 5 ),
398 $ abs( ( norm-pow( 2*n+1 ) ) / pow( 2*n+
399 $ 1 ) ) )
400 DO 400 i = 1, n
401 reslts( 5 ) = max( reslts( 5 ),
402 $ abs( ( r( i )-rpow( i+1 ) ) /
403 $ rpow( i+1 ) ) )
404 400 CONTINUE
405 END IF
406 END IF
407 IF( n.NE.0 ) THEN
408 ab( kl+1, max( n-1, 1 ) ) = -one
409 CALL dpbequ(
'U', n, kl, ab, nszb, r, rcond, norm, info )
410 IF( info.NE.max( n-1, 1 ) )
411 $ reslts( 5 ) = one
412 END IF
413
414
415
416 DO 420 j = 1, nsz
417 DO 410 i = 1, nszb
418 ab( i, j ) = zero
419 410 CONTINUE
420 420 CONTINUE
421 DO 430 j = 1, n
422 ab( 1, j ) = pow( 2*j+1 )
423 430 CONTINUE
424
425 CALL dpbequ(
'L', n, kl, ab, nszb, r, rcond, norm, info )
426
427 IF( info.NE.0 ) THEN
428 reslts( 5 ) = one
429 ELSE
430 IF( n.NE.0 ) THEN
431 reslts( 5 ) = max( reslts( 5 ),
432 $ abs( ( rcond-rpow( n ) ) / rpow( n ) ) )
433 reslts( 5 ) = max( reslts( 5 ),
434 $ abs( ( norm-pow( 2*n+1 ) ) / pow( 2*n+
435 $ 1 ) ) )
436 DO 440 i = 1, n
437 reslts( 5 ) = max( reslts( 5 ),
438 $ abs( ( r( i )-rpow( i+1 ) ) /
439 $ rpow( i+1 ) ) )
440 440 CONTINUE
441 END IF
442 END IF
443 IF( n.NE.0 ) THEN
444 ab( 1, max( n-1, 1 ) ) = -one
445 CALL dpbequ(
'L', n, kl, ab, nszb, r, rcond, norm, info )
446 IF( info.NE.max( n-1, 1 ) )
447 $ reslts( 5 ) = one
448 END IF
449 450 CONTINUE
450 460 CONTINUE
451 reslts( 5 ) = reslts( 5 ) / eps
452 ok = ( reslts( 1 ).LE.thresh ) .AND.
453 $ ( reslts( 2 ).LE.thresh ) .AND.
454 $ ( reslts( 3 ).LE.thresh ) .AND.
455 $ ( reslts( 4 ).LE.thresh ) .AND. ( reslts( 5 ).LE.thresh )
456 WRITE( nout, fmt = * )
457 IF( ok ) THEN
458 WRITE( nout, fmt = 9999 )path
459 ELSE
460 IF( reslts( 1 ).GT.thresh )
461 $ WRITE( nout, fmt = 9998 )reslts( 1 ), thresh
462 IF( reslts( 2 ).GT.thresh )
463 $ WRITE( nout, fmt = 9997 )reslts( 2 ), thresh
464 IF( reslts( 3 ).GT.thresh )
465 $ WRITE( nout, fmt = 9996 )reslts( 3 ), thresh
466 IF( reslts( 4 ).GT.thresh )
467 $ WRITE( nout, fmt = 9995 )reslts( 4 ), thresh
468 IF( reslts( 5 ).GT.thresh )
469 $ WRITE( nout, fmt = 9994 )reslts( 5 ), thresh
470 END IF
471 9999 FORMAT( 1x, 'All tests for ', a3,
472 $ ' routines passed the threshold' )
473 9998 FORMAT( ' DGEEQU failed test with value ', d10.3, ' exceeding',
474 $ ' threshold ', d10.3 )
475 9997 FORMAT( ' DGBEQU failed test with value ', d10.3, ' exceeding',
476 $ ' threshold ', d10.3 )
477 9996 FORMAT( ' DPOEQU failed test with value ', d10.3, ' exceeding',
478 $ ' threshold ', d10.3 )
479 9995 FORMAT( ' DPPEQU failed test with value ', d10.3, ' exceeding',
480 $ ' threshold ', d10.3 )
481 9994 FORMAT( ' DPBEQU failed test with value ', d10.3, ' exceeding',
482 $ ' threshold ', d10.3 )
483 RETURN
484
485
486
subroutine dgbequ(m, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd, amax, info)
DGBEQU
subroutine dgeequ(m, n, a, lda, r, c, rowcnd, colcnd, amax, info)
DGEEQU
double precision function dlamch(cmach)
DLAMCH
subroutine dpbequ(uplo, n, kd, ab, ldab, s, scond, amax, info)
DPBEQU
subroutine dpoequ(n, a, lda, s, scond, amax, info)
DPOEQU
subroutine dppequ(uplo, n, ap, s, scond, amax, info)
DPPEQU