167 SUBROUTINE dget22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR,
176 CHARACTER transa, transe, transw
180 DOUBLE PRECISION a( lda, * ), e(
lde, * ), result( 2 ), wi( * ),
187 DOUBLE PRECISION zero, one
188 parameter( zero = 0.0d0, one = 1.0d0 )
191 CHARACTER norma, norme
192 INTEGER iecol, ierow, ince, ipair, itrnse,
j, jcol,
194 DOUBLE PRECISION anorm, enorm, enrmax, enrmin, errnrm, temp1,
198 DOUBLE PRECISION wmat( 2, 2 )
209 INTRINSIC abs, dble, max, min
220 unfl =
dlamch(
'Safe minimum' )
221 ulp =
dlamch(
'Precision' )
228 IF(
lsame( transa,
'T' ) .OR.
lsame( transa,
'C' ) )
THEN
231 IF(
lsame( transe,
'T' ) .OR.
lsame( transe,
'C' ) )
THEN
241 IF( itrnse.EQ.0 )
THEN
248 IF( ipair.EQ.0 .AND. jvec.LT.n .AND. wi( jvec ).NE.zero )
250 IF( ipair.EQ.1 )
THEN
255 temp1 = max( temp1, abs( e(
j, jvec ) )+
256 $ abs( e(
j, jvec+1 ) ) )
258 enrmin = min( enrmin, temp1 )
259 enrmax = max( enrmax, temp1 )
261 ELSE IF( ipair.EQ.2 )
THEN
268 temp1 = max( temp1, abs( e(
j, jvec ) ) )
270 enrmin = min( enrmin, temp1 )
271 enrmax = max( enrmax, temp1 )
287 IF( ipair.EQ.0 .AND. jvec.LT.n .AND. wi( jvec ).NE.zero )
289 IF( ipair.EQ.1 )
THEN
290 work( jvec ) = max( work( jvec ),
291 $ abs( e(
j, jvec ) )+abs( e(
j,
293 work( jvec+1 ) = work( jvec )
294 ELSE IF( ipair.EQ.2 )
THEN
297 work( jvec ) = max( work( jvec ),
298 $ abs( e(
j, jvec ) ) )
305 enrmin = min( enrmin, work( jvec ) )
306 enrmax = max( enrmax, work( jvec ) )
312 anorm = max(
dlange( norma, n, n, a, lda, work ), unfl )
316 enorm = max(
dlange( norme, n, n, e,
lde, work ), ulp )
322 CALL
dlaset(
'Full', n, n, zero, zero, work, n )
329 IF( itrnse.EQ.1 )
THEN
335 IF( ipair.EQ.0 .AND. wi( jcol ).NE.zero )
338 IF( ipair.EQ.1 )
THEN
339 wmat( 1, 1 ) = wr( jcol )
340 wmat( 2, 1 ) = -wi( jcol )
341 wmat( 1, 2 ) = wi( jcol )
342 wmat( 2, 2 ) = wr( jcol )
343 CALL
dgemm( transe, transw, n, 2, 2, one, e( ierow, iecol ),
344 $
lde, wmat, 2, zero, work( n*( jcol-1 )+1 ), n )
346 ELSE IF( ipair.EQ.2 )
THEN
351 CALL
daxpy( n, wr( jcol ), e( ierow, iecol ), ince,
352 $ work( n*( jcol-1 )+1 ), 1 )
358 CALL
dgemm( transa, transe, n, n, n, one, a, lda, e,
lde, -one,
361 errnrm =
dlange(
'One', n, n, work, n, work( n*n+1 ) ) / enorm
365 IF( anorm.GT.errnrm )
THEN
366 result( 1 ) = ( errnrm / anorm ) / ulp
368 IF( anorm.LT.one )
THEN
369 result( 1 ) = ( min( errnrm, anorm ) / anorm ) / ulp
371 result( 1 ) = min( errnrm / anorm, one ) / ulp
377 result( 2 ) = max( abs( enrmax-one ), abs( enrmin-one ) ) /
LOGICAL function lsame(CA, CB)
LSAME
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
LOGICAL function lde(RI, RJ, LR)
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 ...
subroutine dget22(TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR, WI, WORK, RESULT)
DGET22
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...