64 DOUBLE PRECISION thresh
70 DOUBLE PRECISION zero, one, ten
71 parameter ( zero = 0.0d0, one = 1.0d+0, ten = 1.0d1 )
73 parameter ( nsz = 5, nszb = 3*nsz-2 )
75 parameter ( nszp = ( nsz*( nsz+1 ) ) / 2,
81 INTEGER i, info, j, kl, ku, m, n
82 DOUBLE PRECISION ccond, eps, norm, ratio, rcmax, rcmin, rcond
85 DOUBLE PRECISION a( nsz, nsz ), ab( nszb, nsz ), ap( nszp ),
86 $ c( nsz ), pow( npow ), r( nsz ), reslts( 5 ),
97 INTRINSIC abs, max, min
101 path( 1: 1 ) =
'Double precision'
109 pow( i ) = ten**( i-1 )
110 rpow( i ) = one / pow( i )
120 IF( i.LE.m .AND. j.LE.n )
THEN
121 a( i, j ) = pow( i+j+1 )*( -1 )**( i+j )
128 CALL dgeequ( m, n, a, nsz, r, c, rcond, ccond, norm, info )
133 IF( n.NE.0 .AND. m.NE.0 )
THEN
134 reslts( 1 ) = max( reslts( 1 ),
135 $ abs( ( rcond-rpow( m ) ) / rpow( m ) ) )
136 reslts( 1 ) = max( reslts( 1 ),
137 $ abs( ( ccond-rpow( n ) ) / rpow( n ) ) )
138 reslts( 1 ) = max( reslts( 1 ),
139 $ abs( ( norm-pow( n+m+1 ) ) / pow( n+m+
142 reslts( 1 ) = max( reslts( 1 ),
143 $ abs( ( r( i )-rpow( i+n+1 ) ) /
147 reslts( 1 ) = max( reslts( 1 ),
148 $ abs( ( c( j )-pow( n-j+1 ) ) /
160 a( max( nsz-1, 1 ), j ) = zero
162 CALL dgeequ( nsz, nsz, a, nsz, r, c, rcond, ccond, norm, info )
163 IF( info.NE.max( nsz-1, 1 ) )
167 a( max( nsz-1, 1 ), j ) = one
170 a( i, max( nsz-1, 1 ) ) = zero
172 CALL dgeequ( nsz, nsz, a, nsz, r, c, rcond, ccond, norm, info )
173 IF( info.NE.nsz+max( nsz-1, 1 ) )
175 reslts( 1 ) = reslts( 1 ) / eps
181 DO 230 kl = 0, max( m-1, 0 )
182 DO 220 ku = 0, max( n-1, 0 )
191 IF( i.LE.min( m, j+kl ) .AND. i.GE.
192 $ max( 1, j-ku ) .AND. j.LE.n )
THEN
193 ab( ku+1+i-j, j ) = pow( i+j+1 )*
199 CALL dgbequ( m, n, kl, ku, ab, nszb, r, c, rcond,
200 $ ccond, norm, info )
203 IF( .NOT.( ( n+kl.LT.m .AND. info.EQ.n+kl+1 ) .OR.
204 $ ( m+ku.LT.n .AND. info.EQ.2*m+ku+1 ) ) )
THEN
208 IF( n.NE.0 .AND. m.NE.0 )
THEN
213 rcmin = min( rcmin, r( i ) )
214 rcmax = max( rcmax, r( i ) )
216 ratio = rcmin / rcmax
217 reslts( 2 ) = max( reslts( 2 ),
218 $ abs( ( rcond-ratio ) / ratio ) )
223 rcmin = min( rcmin, c( j ) )
224 rcmax = max( rcmax, c( j ) )
226 ratio = rcmin / rcmax
227 reslts( 2 ) = max( reslts( 2 ),
228 $ abs( ( ccond-ratio ) / ratio ) )
230 reslts( 2 ) = max( reslts( 2 ),
231 $ abs( ( norm-pow( n+m+1 ) ) /
236 IF( i.LE.j+kl .AND. i.GE.j-ku )
THEN
237 ratio = abs( r( i )*pow( i+j+1 )*
239 rcmax = max( rcmax, ratio )
242 reslts( 2 ) = max( reslts( 2 ),
249 IF( i.LE.j+kl .AND. i.GE.j-ku )
THEN
250 ratio = abs( r( i )*pow( i+j+1 )*
252 rcmax = max( rcmax, ratio )
255 reslts( 2 ) = max( reslts( 2 ),
265 reslts( 2 ) = reslts( 2 ) / eps
273 IF( i.LE.n .AND. j.EQ.i )
THEN
274 a( i, j ) = pow( i+j+1 )*( -1 )**( i+j )
281 CALL dpoequ( n, a, nsz, r, rcond, norm, info )
287 reslts( 3 ) = max( reslts( 3 ),
288 $ abs( ( rcond-rpow( n ) ) / rpow( n ) ) )
289 reslts( 3 ) = max( reslts( 3 ),
290 $ abs( ( norm-pow( 2*n+1 ) ) / pow( 2*n+
293 reslts( 3 ) = max( reslts( 3 ),
294 $ abs( ( r( i )-rpow( i+1 ) ) / rpow( i+
300 a( max( nsz-1, 1 ), max( nsz-1, 1 ) ) = -one
301 CALL dpoequ( nsz, a, nsz, r, rcond, norm, info )
302 IF( info.NE.max( nsz-1, 1 ) )
304 reslts( 3 ) = reslts( 3 ) / eps
312 DO 300 i = 1, ( n*( n+1 ) ) / 2
316 ap( ( i*( i+1 ) ) / 2 ) = pow( 2*i+1 )
319 CALL dppequ(
'U', n, ap, r, rcond, norm, info )
325 reslts( 4 ) = max( reslts( 4 ),
326 $ abs( ( rcond-rpow( n ) ) / rpow( n ) ) )
327 reslts( 4 ) = max( reslts( 4 ),
328 $ abs( ( norm-pow( 2*n+1 ) ) / pow( 2*n+
331 reslts( 4 ) = max( reslts( 4 ),
332 $ abs( ( r( i )-rpow( i+1 ) ) / rpow( i+
340 DO 330 i = 1, ( n*( n+1 ) ) / 2
345 ap( j ) = pow( 2*i+1 )
349 CALL dppequ(
'L', n, ap, r, rcond, norm, info )
355 reslts( 4 ) = max( reslts( 4 ),
356 $ abs( ( rcond-rpow( n ) ) / rpow( n ) ) )
357 reslts( 4 ) = max( reslts( 4 ),
358 $ abs( ( norm-pow( 2*n+1 ) ) / pow( 2*n+
361 reslts( 4 ) = max( reslts( 4 ),
362 $ abs( ( r( i )-rpow( i+1 ) ) / rpow( i+
369 i = ( nsz*( nsz+1 ) ) / 2 - 2
371 CALL dppequ(
'L', nsz, ap, r, rcond, norm, info )
372 IF( info.NE.max( nsz-1, 1 ) )
374 reslts( 4 ) = reslts( 4 ) / eps
379 DO 450 kl = 0, max( n-1, 0 )
389 ab( kl+1, j ) = pow( 2*j+1 )
392 CALL dpbequ(
'U', n, kl, ab, nszb, r, rcond, norm, info )
398 reslts( 5 ) = max( reslts( 5 ),
399 $ abs( ( rcond-rpow( n ) ) / rpow( n ) ) )
400 reslts( 5 ) = max( reslts( 5 ),
401 $ abs( ( norm-pow( 2*n+1 ) ) / pow( 2*n+
404 reslts( 5 ) = max( reslts( 5 ),
405 $ abs( ( r( i )-rpow( i+1 ) ) /
411 ab( kl+1, max( n-1, 1 ) ) = -one
412 CALL dpbequ(
'U', n, kl, ab, nszb, r, rcond, norm, info )
413 IF( info.NE.max( n-1, 1 ) )
425 ab( 1, j ) = pow( 2*j+1 )
428 CALL dpbequ(
'L', n, kl, ab, nszb, r, rcond, norm, info )
434 reslts( 5 ) = max( reslts( 5 ),
435 $ abs( ( rcond-rpow( n ) ) / rpow( n ) ) )
436 reslts( 5 ) = max( reslts( 5 ),
437 $ abs( ( norm-pow( 2*n+1 ) ) / pow( 2*n+
440 reslts( 5 ) = max( reslts( 5 ),
441 $ abs( ( r( i )-rpow( i+1 ) ) /
447 ab( 1, max( n-1, 1 ) ) = -one
448 CALL dpbequ(
'L', n, kl, ab, nszb, r, rcond, norm, info )
449 IF( info.NE.max( n-1, 1 ) )
454 reslts( 5 ) = reslts( 5 ) / eps
455 ok = ( reslts( 1 ).LE.thresh ) .AND.
456 $ ( reslts( 2 ).LE.thresh ) .AND.
457 $ ( reslts( 3 ).LE.thresh ) .AND.
458 $ ( reslts( 4 ).LE.thresh ) .AND. ( reslts( 5 ).LE.thresh )
459 WRITE( nout, fmt = * )
461 WRITE( nout, fmt = 9999 )path
463 IF( reslts( 1 ).GT.thresh )
464 $
WRITE( nout, fmt = 9998 )reslts( 1 ), thresh
465 IF( reslts( 2 ).GT.thresh )
466 $
WRITE( nout, fmt = 9997 )reslts( 2 ), thresh
467 IF( reslts( 3 ).GT.thresh )
468 $
WRITE( nout, fmt = 9996 )reslts( 3 ), thresh
469 IF( reslts( 4 ).GT.thresh )
470 $
WRITE( nout, fmt = 9995 )reslts( 4 ), thresh
471 IF( reslts( 5 ).GT.thresh )
472 $
WRITE( nout, fmt = 9994 )reslts( 5 ), thresh
474 9999
FORMAT( 1x,
'All tests for ', a3,
475 $
' routines passed the threshold' )
476 9998
FORMAT(
' DGEEQU failed test with value ', d10.3,
' exceeding',
477 $
' threshold ', d10.3 )
478 9997
FORMAT(
' DGBEQU failed test with value ', d10.3,
' exceeding',
479 $
' threshold ', d10.3 )
480 9996
FORMAT(
' DPOEQU failed test with value ', d10.3,
' exceeding',
481 $
' threshold ', d10.3 )
482 9995
FORMAT(
' DPPEQU failed test with value ', d10.3,
' exceeding',
483 $
' threshold ', d10.3 )
484 9994
FORMAT(
' DPBEQU failed test with value ', d10.3,
' exceeding',
485 $
' threshold ', d10.3 )
double precision function dlamch(CMACH)
DLAMCH
subroutine dgeequ(M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, INFO)
DGEEQU
subroutine dpbequ(UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, INFO)
DPBEQU
subroutine dppequ(UPLO, N, AP, S, SCOND, AMAX, INFO)
DPPEQU
subroutine dgbequ(M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, AMAX, INFO)
DGBEQU
subroutine dpoequ(N, A, LDA, S, SCOND, AMAX, INFO)
DPOEQU