73 IMPLICIT NONE
74
75
76
77
78
79
80 INTEGER M, N, NB, LDT
81
82 REAL RESULT(6)
83
84
85
86
87
88 REAL, ALLOCATABLE :: AF(:,:), Q(:,:),
89 $ R(:,:), RWORK(:), WORK( : ), T(:,:),
90 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:)
91
92
93 REAL ONE, ZERO
94 parameter( zero = 0.0, one = 1.0 )
95
96
97 INTEGER INFO, J, K, L, LWORK
98 REAL ANORM, EPS, RESID, CNORM, DNORM
99
100
101 INTEGER ISEED( 4 )
102
103
105
106
107 REAL SLAMCH
108 REAL SLANGE, SLANSY
109 LOGICAL LSAME
111
112
113 INTRINSIC max, min
114
115
116 DATA iseed / 1988, 1989, 1990, 1991 /
117
119 k = min(m,n)
120 l = max(m,n)
121 lwork = max(2,l)*max(2,l)*nb
122
123
124
125 ALLOCATE ( a(m,n), af(m,n), q(m,m), r(m,l), rwork(l),
126 $ work(lwork), t(nb,n), c(m,n), cf(m,n),
127 $ d(n,m), df(n,m) )
128
129
130
131 ldt=nb
132 DO j=1,n
133 CALL slarnv( 2, iseed, m, a( 1, j ) )
134 END DO
135 CALL slacpy(
'Full', m, n, a, m, af, m )
136
137
138
139 CALL sgeqrt( m, n, nb, af, m, t, ldt, work, info )
140
141
142
143 CALL slaset(
'Full', m, m, zero, one, q, m )
144 CALL sgemqrt(
'R',
'N', m, m, k, nb, af, m, t, ldt, q, m,
145 $ work, info )
146
147
148
149 CALL slaset(
'Full', m, n, zero, zero, r, m )
150 CALL slacpy(
'Upper', m, n, af, m, r, m )
151
152
153
154 CALL sgemm(
'T',
'N', m, n, m, -one, q, m, a, m, one, r, m )
155 anorm =
slange(
'1', m, n, a, m, rwork )
156 resid =
slange(
'1', m, n, r, m, rwork )
157 IF( anorm.GT.zero ) THEN
158 result( 1 ) = resid / (eps*max(1,m)*anorm)
159 ELSE
160 result( 1 ) = zero
161 END IF
162
163
164
165 CALL slaset(
'Full', m, m, zero, one, r, m )
166 CALL ssyrk(
'U',
'C', m, m, -one, q, m, one, r, m )
167 resid =
slansy(
'1',
'Upper', m, r, m, rwork )
168 result( 2 ) = resid / (eps*max(1,m))
169
170
171
172 DO j=1,n
173 CALL slarnv( 2, iseed, m, c( 1, j ) )
174 END DO
175 cnorm =
slange(
'1', m, n, c, m, rwork)
176 CALL slacpy(
'Full', m, n, c, m, cf, m )
177
178
179
180 CALL sgemqrt(
'L',
'N', m, n, k, nb, af, m, t, nb, cf, m,
181 $ work, info)
182
183
184
185 CALL sgemm(
'N',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
186 resid =
slange(
'1', m, n, cf, m, rwork )
187 IF( cnorm.GT.zero ) THEN
188 result( 3 ) = resid / (eps*max(1,m)*cnorm)
189 ELSE
190 result( 3 ) = zero
191 END IF
192
193
194
195 CALL slacpy(
'Full', m, n, c, m, cf, m )
196
197
198
199 CALL sgemqrt(
'L',
'T', m, n, k, nb, af, m, t, nb, cf, m,
200 $ work, info)
201
202
203
204 CALL sgemm(
'T',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
205 resid =
slange(
'1', m, n, cf, m, rwork )
206 IF( cnorm.GT.zero ) THEN
207 result( 4 ) = resid / (eps*max(1,m)*cnorm)
208 ELSE
209 result( 4 ) = zero
210 END IF
211
212
213
214 DO j=1,m
215 CALL slarnv( 2, iseed, n, d( 1, j ) )
216 END DO
217 dnorm =
slange(
'1', n, m, d, n, rwork)
218 CALL slacpy(
'Full', n, m, d, n, df, n )
219
220
221
222 CALL sgemqrt(
'R',
'N', n, m, k, nb, af, m, t, nb, df, n,
223 $ work, info)
224
225
226
227 CALL sgemm(
'N',
'N', n, m, m, -one, d, n, q, m, one, df, n )
228 resid =
slange(
'1', n, m, df, n, rwork )
229 IF( cnorm.GT.zero ) THEN
230 result( 5 ) = resid / (eps*max(1,m)*dnorm)
231 ELSE
232 result( 5 ) = zero
233 END IF
234
235
236
237 CALL slacpy(
'Full', n, m, d, n, df, n )
238
239
240
241 CALL sgemqrt(
'R',
'T', n, m, k, nb, af, m, t, nb, df, n,
242 $ work, info)
243
244
245
246 CALL sgemm(
'N',
'T', n, m, m, -one, d, n, q, m, one, df, n )
247 resid =
slange(
'1', n, m, df, n, rwork )
248 IF( cnorm.GT.zero ) THEN
249 result( 6 ) = resid / (eps*max(1,m)*dnorm)
250 ELSE
251 result( 6 ) = zero
252 END IF
253
254
255
256 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
257
258 RETURN
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
subroutine sgemqrt(side, trans, m, n, k, nb, v, ldv, t, ldt, c, ldc, work, info)
SGEMQRT
subroutine sgeqrt(m, n, nb, a, lda, t, ldt, work, info)
SGEQRT
subroutine ssyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
SSYRK
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
real function slamch(cmach)
SLAMCH
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 ...
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,...
subroutine slarnv(idist, iseed, n, x)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
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