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, n2, np1,i
109 DOUBLE PRECISION anorm, eps, resid, cnorm, dnorm
121 DATA iseed / 1988, 1989, 1990, 1991 /
135 ALLOCATE(a(m,n2),af(m,n2),q(n2,n2),r(n2,n2),rwork(n2),
136 $ work(lwork),t(nb,m),c(n2,m),cf(n2,m),
142 CALL zlaset(
'Full', m, n2, czero, czero, a, m )
143 CALL zlaset(
'Full', nb, m, czero, czero, t, nb )
145 CALL zlarnv( 2, iseed, m-j+1, a( j, j ) )
149 CALL zlarnv( 2, iseed, m, a( 1, min(n+m,m+1) + j - 1 ) )
154 CALL zlarnv( 2, iseed, m-j+1, a( j, min(n+m,n+m-l+1)
161 CALL zlacpy(
'Full', m, n2, a, m, af, m )
165 CALL ztplqt( m,n,l,nb,af,m,af(1,np1),m,t,ldt,work,info)
169 CALL zlaset(
'Full', n2, n2, czero, one, q, n2 )
170 CALL zgemlqt(
'L',
'N', n2, n2, k, nb, af, m, t, ldt, q, n2,
175 CALL zlaset(
'Full', n2, n2, czero, czero, r, n2 )
176 CALL zlacpy(
'Lower', m, n2, af, m, r, n2 )
180 CALL zgemm(
'N',
'C', m, n2, n2, -one, a, m, q, n2, one, r, n2)
181 anorm =
zlange(
'1', m, n2, a, m, rwork )
182 resid =
zlange(
'1', m, n2, r, n2, rwork )
183 IF( anorm.GT.zero )
THEN
184 result( 1 ) = resid / (eps*anorm*max(1,n2))
191 CALL zlaset(
'Full', n2, n2, czero, one, r, n2 )
192 CALL zherk(
'U',
'N', n2, n2, dreal(-one), q, n2, dreal(one),
194 resid =
zlansy(
'1',
'Upper', n2, r, n2, rwork )
195 result( 2 ) = resid / (eps*max(1,n2))
199 CALL zlaset(
'Full', n2, m, czero, one, c, n2 )
201 CALL zlarnv( 2, iseed, n2, c( 1, j ) )
203 cnorm =
zlange(
'1', n2, m, c, n2, rwork)
204 CALL zlacpy(
'Full', n2, m, c, n2, cf, n2 )
208 CALL ztpmlqt(
'L',
'N', n,m,k,l,nb,af(1, np1),m,t,ldt,cf,n2,
209 $ cf(np1,1),n2,work,info)
213 CALL zgemm(
'N',
'N', n2, m, n2, -one, q, n2, c, n2, one, cf, n2 )
214 resid =
zlange(
'1', n2, m, cf, n2, rwork )
215 IF( cnorm.GT.zero )
THEN
216 result( 3 ) = resid / (eps*max(1,n2)*cnorm)
224 CALL zlacpy(
'Full', n2, m, c, n2, cf, n2 )
228 CALL ztpmlqt(
'L',
'C',n,m,k,l,nb,af(1,np1),m,t,ldt,cf,n2,
229 $ cf(np1,1),n2,work,info)
233 CALL zgemm(
'C',
'N',n2,m,n2,-one,q,n2,c,n2,one,cf,n2)
234 resid =
zlange(
'1', n2, m, cf, n2, rwork )
236 IF( cnorm.GT.zero )
THEN
237 result( 4 ) = resid / (eps*max(1,n2)*cnorm)
245 CALL zlarnv( 2, iseed, m, d( 1, j ) )
247 dnorm =
zlange(
'1', m, n2, d, m, rwork)
248 CALL zlacpy(
'Full', m, n2, d, m, df, m )
252 CALL ztpmlqt(
'R',
'N',m,n,k,l,nb,af(1,np1),m,t,ldt,df,m,
253 $ df(1,np1),m,work,info)
257 CALL zgemm(
'N',
'N',m,n2,n2,-one,d,m,q,n2,one,df,m)
258 resid =
zlange(
'1',m, n2,df,m,rwork )
259 IF( cnorm.GT.zero )
THEN
260 result( 5 ) = resid / (eps*max(1,n2)*dnorm)
267 CALL zlacpy(
'Full',m,n2,d,m,df,m )
271 CALL ztpmlqt(
'R',
'C',m,n,k,l,nb,af(1,np1),m,t,ldt,df,m,
272 $ df(1,np1),m,work,info)
277 CALL zgemm(
'N',
'C', m, n2, n2, -one, d, m, q, n2, one, df, m )
278 resid =
zlange(
'1', m, n2, df, m, rwork )
279 IF( cnorm.GT.zero )
THEN
280 result( 6 ) = resid / (eps*max(1,n2)*dnorm)
287 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
double precision function dlamch(CMACH)
DLAMCH
subroutine zgemlqt(SIDE, TRANS, M, N, K, MB, V, LDV, T, LDT, C, LDC, WORK, INFO)
ZGEMLQT
subroutine ztplqt(M, N, L, MB, A, LDA, B, LDB, T, LDT, WORK, INFO)
ZTPLQT
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine ztpmlqt(SIDE, TRANS, M, N, K, L, MB, V, LDV, T, LDT, A, LDA, B, LDB, WORK, INFO)
ZTPMLQT
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
double precision function zlansy(NORM, UPLO, N, A, LDA, WORK)
ZLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex symmetric matrix.
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...
logical function lsame(CA, CB)
LSAME
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
subroutine zlarnv(IDIST, ISEED, N, X)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.