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...
LOGICAL function lsame(CA, CB)
LSAME
REAL function slamch(CMACH)
SLAMCH
subroutine ctpqrt(M, N, L, NB, A, LDA, B, LDB, T, LDT, WORK, INFO)
CTPQRT
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 ...
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine cqrt05(M, N, L, NB, RESULT)
CQRT05
subroutine clarnv(IDIST, ISEED, N, X)
CLARNV returns a vector of random numbers from a uniform or normal distribution.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
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 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