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 $ L(:,:), 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, LL, LWORK
98 REAL ANORM, EPS, RESID, CNORM, DNORM
99
100
101 INTEGER ISEED( 4 )
102
103
104 REAL SLAMCH, SLANGE, SLANSY
105 LOGICAL LSAME
107
108
109 INTRINSIC max, min
110
111
112 DATA iseed / 1988, 1989, 1990, 1991 /
113
115 k = min(m,n)
116 ll = max(m,n)
117 lwork = max(2,ll)*max(2,ll)*nb
118
119
120
121 ALLOCATE ( a(m,n), af(m,n), q(n,n), l(ll,n), rwork(ll),
122 $ work(lwork), t(nb,n), c(m,n), cf(m,n),
123 $ d(n,m), df(n,m) )
124
125
126
127 ldt=nb
128 DO j=1,n
129 CALL slarnv( 2, iseed, m, a( 1, j ) )
130 END DO
131 CALL slacpy(
'Full', m, n, a, m, af, m )
132
133
134
135 CALL sgelqt( m, n, nb, af, m, t, ldt, work, info )
136
137
138
139 CALL slaset(
'Full', n, n, zero, one, q, n )
140 CALL sgemlqt(
'R',
'N', n, n, k, nb, af, m, t, ldt, q, n,
141 $ work, info )
142
143
144
145 CALL slaset(
'Full', m, n, zero, zero, l, ll )
146 CALL slacpy(
'Lower', m, n, af, m, l, ll )
147
148
149
150 CALL sgemm(
'N',
'T', m, n, n, -one, a, m, q, n, one, l, ll )
151 anorm =
slange(
'1', m, n, a, m, rwork )
152 resid =
slange(
'1', m, n, l, ll, rwork )
153 IF( anorm.GT.zero ) THEN
154 result( 1 ) = resid / (eps*max(1,m)*anorm)
155 ELSE
156 result( 1 ) = zero
157 END IF
158
159
160
161 CALL slaset(
'Full', n, n, zero, one, l, ll )
162 CALL ssyrk(
'U',
'C', n, n, -one, q, n, one, l, ll )
163 resid =
slansy(
'1',
'Upper', n, l, ll, rwork )
164 result( 2 ) = resid / (eps*max(1,n))
165
166
167
168 DO j=1,m
169 CALL slarnv( 2, iseed, n, d( 1, j ) )
170 END DO
171 dnorm =
slange(
'1', n, m, d, n, rwork)
172 CALL slacpy(
'Full', n, m, d, n, df, n )
173
174
175
176 CALL sgemlqt(
'L',
'N', n, m, k, nb, af, m, t, nb, df, n,
177 $ work, info)
178
179
180
181 CALL sgemm(
'N',
'N', n, m, n, -one, q, n, d, n, one, df, n )
182 resid =
slange(
'1', n, m, df, n, rwork )
183 IF( dnorm.GT.zero ) THEN
184 result( 3 ) = resid / (eps*max(1,m)*dnorm)
185 ELSE
186 result( 3 ) = zero
187 END IF
188
189
190
191 CALL slacpy(
'Full', n, m, d, n, df, n )
192
193
194
195 CALL sgemlqt(
'L',
'T', n, m, k, nb, af, m, t, nb, df, n,
196 $ work, info)
197
198
199
200 CALL sgemm(
'T',
'N', n, m, n, -one, q, n, d, n, one, df, n )
201 resid =
slange(
'1', n, m, df, n, rwork )
202 IF( dnorm.GT.zero ) THEN
203 result( 4 ) = resid / (eps*max(1,m)*dnorm)
204 ELSE
205 result( 4 ) = zero
206 END IF
207
208
209
210 DO j=1,n
211 CALL slarnv( 2, iseed, m, c( 1, j ) )
212 END DO
213 cnorm =
slange(
'1', m, n, c, m, rwork)
214 CALL slacpy(
'Full', m, n, c, m, cf, m )
215
216
217
218 CALL sgemlqt(
'R',
'N', m, n, k, nb, af, m, t, nb, cf, m,
219 $ work, info)
220
221
222
223 CALL sgemm(
'N',
'N', m, n, n, -one, c, m, q, n, one, cf, m )
224 resid =
slange(
'1', n, m, df, n, rwork )
225 IF( cnorm.GT.zero ) THEN
226 result( 5 ) = resid / (eps*max(1,m)*dnorm)
227 ELSE
228 result( 5 ) = zero
229 END IF
230
231
232
233 CALL slacpy(
'Full', m, n, c, m, cf, m )
234
235
236
237 CALL sgemlqt(
'R',
'T', m, n, k, nb, af, m, t, nb, cf, m,
238 $ work, info)
239
240
241
242 CALL sgemm(
'N',
'T', m, n, n, -one, c, m, q, n, one, cf, m )
243 resid =
slange(
'1', m, n, cf, m, rwork )
244 IF( cnorm.GT.zero ) THEN
245 result( 6 ) = resid / (eps*max(1,m)*dnorm)
246 ELSE
247 result( 6 ) = zero
248 END IF
249
250
251
252 DEALLOCATE ( a, af, q, l, rwork, work, t, c, d, cf, df)
253
254 RETURN
subroutine sgelqt(m, n, mb, a, lda, t, ldt, work, info)
SGELQT
subroutine sgemlqt(side, trans, m, n, k, mb, v, ldv, t, ldt, c, ldc, work, info)
SGEMLQT
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
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