90 INTEGER lwork, m, n, l, nb, ldt
98 COMPLEX,
ALLOCATABLE :: af(:,:), q(:,:),
99 $ r(:,:), rwork(:), work( : ), t(:,:),
100 $ cf(:,:), df(:,:), a(:,:), c(:,:), d(:,:)
105 parameter( zero = 0.0, one = (1.0,0.0), czero=(0.0,0.0) )
108 INTEGER info, j, k, m2, np1
109 REAL 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 claset(
'Full', m2, n, czero, czero, a, m2 )
143 CALL claset(
'Full', nb, n, czero, czero, t, nb )
145 CALL clarnv( 2, iseed, j, a( 1, j ) )
149 CALL clarnv( 2, iseed, m-l, a( min(n+m,n+1), j ) )
154 CALL clarnv( 2, iseed, min(j,l), a( min(n+m,n+m-l+1), j ) )
160 CALL clacpy(
'Full', m2, n, a, m2, af, m2 )
164 CALL ctpqrt( m,n,l,nb,af,m2,af(np1,1),m2,t,ldt,work,info)
168 CALL claset(
'Full', m2, m2, czero, one, q, m2 )
169 CALL cgemqrt(
'R',
'N', m2, m2, k, nb, af, m2, t, ldt, q, m2,
174 CALL claset(
'Full', m2, n, czero, czero, r, m2 )
175 CALL clacpy(
'Upper', m2, n, af, m2, r, m2 )
179 CALL cgemm(
'C',
'N', m2, n, m2, -one, q, m2, a, m2, one, r, m2 )
180 anorm =
clange(
'1', m2, n, a, m2, rwork )
181 resid =
clange(
'1', m2, n, r, m2, rwork )
182 IF( anorm.GT.zero )
THEN
183 result( 1 ) = resid / (eps*anorm*max(1,m2))
190 CALL claset(
'Full', m2, m2, czero, one, r, m2 )
191 CALL cherk(
'U',
'C', m2, m2,
REAL(-ONE), q, m2,
REAL(ONE),
193 resid =
clansy(
'1',
'Upper', m2, r, m2, rwork )
194 result( 2 ) = resid / (eps*max(1,m2))
199 CALL clarnv( 2, iseed, m2, c( 1, j ) )
201 cnorm =
clange(
'1', m2, n, c, m2, rwork)
202 CALL clacpy(
'Full', m2, n, c, m2, cf, m2 )
206 CALL ctpmqrt(
'L',
'N', m,n,k,l,nb,af(np1,1),m2,t,ldt,cf,m2,
207 $ cf(np1,1),m2,work,info)
211 CALL cgemm(
'N',
'N', m2, n, m2, -one, q, m2, c, m2, one, cf, m2 )
212 resid =
clange(
'1', m2, n, cf, m2, rwork )
213 IF( cnorm.GT.zero )
THEN
214 result( 3 ) = resid / (eps*max(1,m2)*cnorm)
221 CALL clacpy(
'Full', m2, n, c, m2, cf, m2 )
225 CALL ctpmqrt(
'L',
'C',m,n,k,l,nb,af(np1,1),m2,t,ldt,cf,m2,
226 $ cf(np1,1),m2,work,info)
230 CALL cgemm(
'C',
'N',m2,n,m2,-one,q,m2,c,m2,one,cf,m2)
231 resid =
clange(
'1', m2, n, cf, m2, rwork )
232 IF( cnorm.GT.zero )
THEN
233 result( 4 ) = resid / (eps*max(1,m2)*cnorm)
241 CALL clarnv( 2, iseed, n, d( 1, j ) )
243 dnorm =
clange(
'1', n, m2, d, n, rwork)
244 CALL clacpy(
'Full', n, m2, d, n, df, n )
248 CALL ctpmqrt(
'R',
'N',n,m,n,l,nb,af(np1,1),m2,t,ldt,df,n,
249 $ df(1,np1),n,work,info)
253 CALL cgemm(
'N',
'N',n,m2,m2,-one,d,n,q,m2,one,df,n)
254 resid =
clange(
'1',n, m2,df,n,rwork )
255 IF( cnorm.GT.zero )
THEN
256 result( 5 ) = resid / (eps*max(1,m2)*dnorm)
263 CALL clacpy(
'Full',n,m2,d,n,df,n )
267 CALL ctpmqrt(
'R',
'C',n,m,n,l,nb,af(np1,1),m2,t,ldt,df,n,
268 $ df(1,np1),n,work,info)
273 CALL cgemm(
'N',
'C', n, m2, m2, -one, d, n, q, m2, one, df, n )
274 resid =
clange(
'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 claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
logical function lsame(CA, CB)
LSAME
subroutine ctpqrt(M, N, L, NB, A, LDA, B, LDB, T, LDT, WORK, INFO)
CTPQRT
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
real function slamch(CMACH)
SLAMCH
subroutine clarnv(IDIST, ISEED, N, X)
CLARNV returns a vector of random numbers from a uniform or normal distribution.
real function clansy(NORM, UPLO, N, A, LDA, WORK)
CLANSY 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 cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
subroutine ctpmqrt(SIDE, TRANS, M, N, K, L, NB, V, LDV, T, LDT, A, LDA, B, LDB, WORK, INFO)
CTPMQRT
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
subroutine cgemqrt(SIDE, TRANS, M, N, K, NB, V, LDV, T, LDT, C, LDC, WORK, INFO)
CGEMQRT