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
121 DATA iseed / 1988, 1989, 1990, 1991 /
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
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 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
logical function lsame(CA, CB)
LSAME
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 zlarnv(IDIST, ISEED, N, X)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
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 dlamch(CMACH)
DLAMCH
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