90 INTEGER lwork, m, n, l, nb, ldt
98 REAL,
ALLOCATABLE :: af(:,:), q(:,:),
99 $ r(:,:), rwork(:), work( : ), t(:,:),
100 $ cf(:,:), df(:,:), a(:,:), c(:,:), d(:,:)
104 parameter( zero = 0.0, one = 1.0 )
107 INTEGER info,
j, k, m2, np1
108 REAL anorm, eps, resid, cnorm, dnorm
120 DATA iseed / 1988, 1989, 1990, 1991 /
134 ALLOCATE(a(m2,n),af(m2,n),q(m2,m2),r(m2,m2),rwork(m2),
135 $ work(lwork),t(nb,n),c(m2,n),cf(m2,n),
141 CALL
slaset(
'Full', m2, n, zero, zero, a, m2 )
142 CALL
slaset(
'Full', nb, n, zero, zero, t, nb )
144 CALL
slarnv( 2, iseed,
j, a( 1,
j ) )
148 CALL
slarnv( 2, iseed, m-l, a( n+1,
j ) )
153 CALL
slarnv( 2, iseed, min(
j,l), a( n+m-l+1,
j ) )
159 CALL
slacpy(
'Full', m2, n, a, m2, af, m2 )
163 CALL
stpqrt( m,n,l,nb,af,m2,af(np1,1),m2,t,ldt,work,info)
167 CALL
slaset(
'Full', m2, m2, zero, one, q, m2 )
168 CALL
sgemqrt(
'R',
'N', m2, m2, k, nb, af, m2, t, ldt, q, m2,
173 CALL
slaset(
'Full', m2, n, zero, zero, r, m2 )
174 CALL
slacpy(
'Upper', m2, n, af, m2, r, m2 )
178 CALL
sgemm(
'T',
'N', m2, n, m2, -one, q, m2, a, m2, one, r, m2 )
179 anorm =
slange(
'1', m2, n, a, m2, rwork )
180 resid =
slange(
'1', m2, n, r, m2, rwork )
181 IF( anorm.GT.zero )
THEN
182 result( 1 ) = resid / (eps*anorm*max(1,m2))
189 CALL
slaset(
'Full', m2, m2, zero, one, r, m2 )
190 CALL
ssyrk(
'U',
'C', m2, m2, -one, q, m2, one,
192 resid =
slansy(
'1',
'Upper', m2, r, m2, rwork )
193 result( 2 ) = resid / (eps*max(1,m2))
198 CALL
slarnv( 2, iseed, m2, c( 1,
j ) )
200 cnorm =
slange(
'1', m2, n, c, m2, rwork)
201 CALL
slacpy(
'Full', m2, n, c, m2, cf, m2 )
205 CALL
stpmqrt(
'L',
'N', m,n,k,l,nb,af(np1,1),m2,t,ldt,cf,
206 $ m2,cf(np1,1),m2,work,info)
210 CALL
sgemm(
'N',
'N', m2, n, m2, -one, q,m2,c,m2,one,cf,m2)
211 resid =
slange(
'1', m2, n, cf, m2, rwork )
212 IF( cnorm.GT.zero )
THEN
213 result( 3 ) = resid / (eps*max(1,m2)*cnorm)
220 CALL
slacpy(
'Full', m2, n, c, m2, cf, m2 )
224 CALL
stpmqrt(
'L',
'T',m,n,k,l,nb,af(np1,1),m2,t,ldt,cf,m2,
225 $ cf(np1,1),m2,work,info)
229 CALL
sgemm(
'T',
'N',m2,n,m2,-one,q,m2,c,m2,one,cf,m2)
230 resid =
slange(
'1', m2, n, cf, m2, rwork )
231 IF( cnorm.GT.zero )
THEN
232 result( 4 ) = resid / (eps*max(1,m2)*cnorm)
240 CALL
slarnv( 2, iseed, n, d( 1,
j ) )
242 dnorm =
slange(
'1', n, m2, d, n, rwork)
243 CALL
slacpy(
'Full', n, m2, d, n, df, n )
247 CALL
stpmqrt(
'R',
'N',n,m,n,l,nb,af(np1,1),m2,t,ldt,df,n,
248 $ df(1,np1),n,work,info)
252 CALL
sgemm(
'N',
'N',n,m2,m2,-one,d,n,q,m2,one,df,n)
253 resid =
slange(
'1',n, m2,df,n,rwork )
254 IF( cnorm.GT.zero )
THEN
255 result( 5 ) = resid / (eps*max(1,m2)*dnorm)
262 CALL
slacpy(
'Full',n,m2,d,n,df,n )
266 CALL
stpmqrt(
'R',
'T',n,m,n,l,nb,af(np1,1),m2,t,ldt,df,n,
267 $ df(1,np1),n,work,info)
272 CALL
sgemm(
'N',
'T', n, m2, m2, -one, d, n, q, m2, one, df, n )
273 resid =
slange(
'1', n, m2, df, n, rwork )
274 IF( cnorm.GT.zero )
THEN
275 result( 6 ) = resid / (eps*max(1,m2)*dnorm)
282 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
subroutine sgemqrt(SIDE, TRANS, M, N, K, NB, V, LDV, T, LDT, C, LDC, WORK, INFO)
SGEMQRT
subroutine sqrt05(M, N, L, NB, RESULT)
SQRT05
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
REAL function slansy(NORM, UPLO, N, A, LDA, WORK)
SLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric matrix.
LOGICAL function lsame(CA, CB)
LSAME
REAL function slamch(CMACH)
SLAMCH
subroutine stpqrt(M, N, L, NB, A, LDA, B, LDB, T, LDT, WORK, INFO)
STPQRT
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
subroutine stpmqrt(SIDE, TRANS, M, N, K, L, NB, V, LDV, T, LDT, A, LDA, B, LDB, WORK, INFO)
STPMQRT
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
REAL function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine slarnv(IDIST, ISEED, N, X)
SLARNV 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