91 REAL,
ALLOCATABLE :: af(:,:), q(:,:),
92 $ l(:,:), rwork(:), work( : ), t(:,:),
93 $ cf(:,:), df(:,:), a(:,:), c(:,:), d(:,:)
97 parameter( zero = 0.0, one = 1.0 )
100 INTEGER info, j, k, ll, lwork
101 REAL anorm, eps, resid, cnorm, dnorm
115 DATA iseed / 1988, 1989, 1990, 1991 /
120 lwork = max(2,ll)*max(2,ll)*nb
124 ALLOCATE ( a(m,n), af(m,n), q(n,n), l(ll,n), rwork(ll),
125 $ work(lwork), t(nb,n), c(m,n), cf(m,n),
132 CALL slarnv( 2, iseed, m, a( 1, j ) )
134 CALL slacpy(
'Full', m, n, a, m, af, m )
138 CALL sgelqt( m, n, nb, af, m, t, ldt, work, info )
142 CALL slaset(
'Full', n, n, zero, one, q, n )
143 CALL sgemlqt(
'R',
'N', n, n, k, nb, af, m, t, ldt, q, n,
148 CALL slaset(
'Full', m, n, zero, zero, l, ll )
149 CALL slacpy(
'Lower', m, n, af, m, l, ll )
153 CALL sgemm(
'N',
'T', m, n, n, -one, a, m, q, n, one, l, ll )
154 anorm =
slange(
'1', m, n, a, m, rwork )
155 resid =
slange(
'1', m, n, l, ll, rwork )
156 IF( anorm.GT.zero )
THEN
157 result( 1 ) = resid / (eps*max(1,m)*anorm)
164 CALL slaset(
'Full', n, n, zero, one, l, ll )
165 CALL ssyrk(
'U',
'C', n, n, -one, q, n, one, l, ll )
166 resid =
slansy(
'1',
'Upper', n, l, ll, rwork )
167 result( 2 ) = resid / (eps*max(1,n))
172 CALL slarnv( 2, iseed, n, d( 1, j ) )
174 dnorm =
slange(
'1', n, m, d, n, rwork)
175 CALL slacpy(
'Full', n, m, d, n, df, n )
179 CALL sgemlqt(
'L',
'N', n, m, k, nb, af, m, t, nb, df, n,
184 CALL sgemm(
'N',
'N', n, m, n, -one, q, n, d, n, one, df, n )
185 resid =
slange(
'1', n, m, df, n, rwork )
186 IF( dnorm.GT.zero )
THEN
187 result( 3 ) = resid / (eps*max(1,m)*dnorm)
194 CALL slacpy(
'Full', n, m, d, n, df, n )
198 CALL sgemlqt(
'L',
'T', n, m, k, nb, af, m, t, nb, df, n,
203 CALL sgemm(
'T',
'N', n, m, n, -one, q, n, d, n, one, df, n )
204 resid =
slange(
'1', n, m, df, n, rwork )
205 IF( dnorm.GT.zero )
THEN
206 result( 4 ) = resid / (eps*max(1,m)*dnorm)
214 CALL slarnv( 2, iseed, m, c( 1, j ) )
216 cnorm =
slange(
'1', m, n, c, m, rwork)
217 CALL slacpy(
'Full', m, n, c, m, cf, m )
221 CALL sgemlqt(
'R',
'N', m, n, k, nb, af, m, t, nb, cf, m,
226 CALL sgemm(
'N',
'N', m, n, n, -one, c, m, q, n, one, cf, m )
227 resid =
slange(
'1', n, m, df, n, rwork )
228 IF( cnorm.GT.zero )
THEN
229 result( 5 ) = resid / (eps*max(1,m)*dnorm)
236 CALL slacpy(
'Full', m, n, c, m, cf, m )
240 CALL sgemlqt(
'R',
'T', m, n, k, nb, af, m, t, nb, cf, m,
245 CALL sgemm(
'N',
'T', m, n, n, -one, c, m, q, n, one, cf, m )
246 resid =
slange(
'1', m, n, cf, m, rwork )
247 IF( cnorm.GT.zero )
THEN
248 result( 6 ) = resid / (eps*max(1,m)*dnorm)
255 DEALLOCATE ( a, af, q, l, rwork, work, t, c, d, cf, df)
subroutine sgelqt(M, N, MB, A, LDA, T, LDT, WORK, INFO)
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
subroutine sgemlqt(SIDE, TRANS, M, N, K, MB, V, LDV, T, LDT, C, LDC, WORK, INFO)
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
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 slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET 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 slarnv(IDIST, ISEED, N, X)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
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.