90 INTEGER LWORK, M, N, L, NB, LDT
92 DOUBLE PRECISION RESULT(6)
98 COMPLEX*16,
ALLOCATABLE :: AF(:,:), Q(:,:),
99 $ r(:,:), rwork(:), work( : ), t(:,:),
100 $ cf(:,:), df(:,:), a(:,:), c(:,:), d(:,:)
103 DOUBLE PRECISION ZERO
104 COMPLEX*16 ONE, CZERO
105 parameter( zero = 0.0, one = (1.0,0.0), czero=(0.0,0.0) )
108 INTEGER INFO, J, K, M2, NP1
109 DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
115 DOUBLE PRECISION DLAMCH
116 DOUBLE PRECISION ZLANGE, ZLANSY
118 EXTERNAL dlamch, zlange, zlansy, lsame
121 DATA iseed / 1988, 1989, 1990, 1991 /
123 eps = dlamch(
'Epsilon' )
135 ALLOCATE(a(m2,n),af(m2,n),q(m2,m2),r(m2,m2),rwork(m2),
136 $ work(lwork),t(nb,n),c(m2,n),cf(m2,n),
142 CALL
zlaset(
'Full', m2, n, czero, czero, a, m2 )
143 CALL
zlaset(
'Full', nb, n, czero, czero, t, nb )
145 CALL
zlarnv( 2, iseed, j, a( 1, j ) )
149 CALL
zlarnv( 2, iseed, m-l, a( min(n+m,n+1), j ) )
154 CALL
zlarnv( 2, iseed, min(j,l), a( min(n+m,n+m-l+1), j ) )
160 CALL
zlacpy(
'Full', m2, n, a, m2, af, m2 )
164 CALL
ztpqrt( m,n,l,nb,af,m2,af(np1,1),m2,t,ldt,work,info)
168 CALL
zlaset(
'Full', m2, m2, czero, one, q, m2 )
169 CALL
zgemqrt(
'R',
'N', m2, m2, k, nb, af, m2, t, ldt, q, m2,
174 CALL
zlaset(
'Full', m2, n, czero, czero, r, m2 )
175 CALL
zlacpy(
'Upper', m2, n, af, m2, r, m2 )
179 CALL
zgemm(
'C',
'N', m2, n, m2, -one, q, m2, a, m2, one, r, m2 )
180 anorm = zlange(
'1', m2, n, a, m2, rwork )
181 resid = zlange(
'1', m2, n, r, m2, rwork )
182 IF( anorm.GT.zero )
THEN
183 result( 1 ) = resid / (eps*anorm*max(1,m2))
190 CALL
zlaset(
'Full', m2, m2, czero, one, r, m2 )
191 CALL
zherk(
'U',
'C', m2, m2, dreal(-one), q, m2, dreal(one),
193 resid = zlansy(
'1',
'Upper', m2, r, m2, rwork )
194 result( 2 ) = resid / (eps*max(1,m2))
199 CALL
zlarnv( 2, iseed, m2, c( 1, j ) )
201 cnorm = zlange(
'1', m2, n, c, m2, rwork)
202 CALL
zlacpy(
'Full', m2, n, c, m2, cf, m2 )
206 CALL
ztpmqrt(
'L',
'N', m,n,k,l,nb,af(np1,1),m2,t,ldt,cf,m2,
207 $ cf(np1,1),m2,work,info)
211 CALL
zgemm(
'N',
'N', m2, n, m2, -one, q, m2, c, m2, one, cf, m2 )
212 resid = zlange(
'1', m2, n, cf, m2, rwork )
213 IF( cnorm.GT.zero )
THEN
214 result( 3 ) = resid / (eps*max(1,m2)*cnorm)
221 CALL
zlacpy(
'Full', m2, n, c, m2, cf, m2 )
225 CALL
ztpmqrt(
'L',
'C',m,n,k,l,nb,af(np1,1),m2,t,ldt,cf,m2,
226 $ cf(np1,1),m2,work,info)
230 CALL
zgemm(
'C',
'N',m2,n,m2,-one,q,m2,c,m2,one,cf,m2)
231 resid = zlange(
'1', m2, n, cf, m2, rwork )
232 IF( cnorm.GT.zero )
THEN
233 result( 4 ) = resid / (eps*max(1,m2)*cnorm)
241 CALL
zlarnv( 2, iseed, n, d( 1, j ) )
243 dnorm = zlange(
'1', n, m2, d, n, rwork)
244 CALL
zlacpy(
'Full', n, m2, d, n, df, n )
248 CALL
ztpmqrt(
'R',
'N',n,m,n,l,nb,af(np1,1),m2,t,ldt,df,n,
249 $ df(1,np1),n,work,info)
253 CALL
zgemm(
'N',
'N',n,m2,m2,-one,d,n,q,m2,one,df,n)
254 resid = zlange(
'1',n, m2,df,n,rwork )
255 IF( cnorm.GT.zero )
THEN
256 result( 5 ) = resid / (eps*max(1,m2)*dnorm)
263 CALL
zlacpy(
'Full',n,m2,d,n,df,n )
267 CALL
ztpmqrt(
'R',
'C',n,m,n,l,nb,af(np1,1),m2,t,ldt,df,n,
268 $ df(1,np1),n,work,info)
273 CALL
zgemm(
'N',
'C', n, m2, m2, -one, d, n, q, m2, one, df, n )
274 resid = zlange(
'1', n, m2, df, n, rwork )
275 IF( cnorm.GT.zero )
THEN
276 result( 6 ) = resid / (eps*max(1,m2)*dnorm)
283 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zgemqrt(SIDE, TRANS, M, N, K, NB, V, LDV, T, LDT, C, LDC, WORK, INFO)
ZGEMQRT
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine zqrt05(M, N, L, NB, RESULT)
ZQRT05
subroutine zlarnv(IDIST, ISEED, N, X)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine ztpqrt(M, N, L, NB, A, LDA, B, LDB, T, LDT, WORK, INFO)
ZTPQRT
subroutine ztpmqrt(SIDE, TRANS, M, N, K, L, NB, V, LDV, T, LDT, A, LDA, B, LDB, WORK, INFO)
ZTPMQRT