81 SUBROUTINE sqrt05(M,N,L,NB,RESULT)
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
117 EXTERNAL slamch, slange, slansy, lsame
120 DATA iseed / 1988, 1989, 1990, 1991 /
122 eps = slamch(
'Epsilon' )
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 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 sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine stpmqrt(SIDE, TRANS, M, N, K, L, NB, V, LDV, T, LDT, A, LDA, B, LDB, WORK, INFO)
STPMQRT
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...
subroutine slarnv(IDIST, ISEED, N, X)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine sgemqrt(SIDE, TRANS, M, N, K, NB, V, LDV, T, LDT, C, LDC, WORK, INFO)
SGEMQRT
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.