79 IMPLICIT NONE
80
81
82
83
84
85
86 INTEGER LWORK, M, N, L, NB, LDT
87
88 REAL RESULT(6)
89
90
91
92
93
94 REAL, ALLOCATABLE :: AF(:,:), Q(:,:),
95 $ R(:,:), RWORK(:), WORK( : ), T(:,:),
96 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:)
97
98
99 REAL ONE, ZERO
100 parameter( zero = 0.0, one = 1.0 )
101
102
103 INTEGER INFO, J, K, N2, NP1,i
104 REAL ANORM, EPS, RESID, CNORM, DNORM
105
106
107 INTEGER ISEED( 4 )
108
109
110 REAL SLAMCH, SLANGE, SLANSY
111 LOGICAL LSAME
113
114
115 DATA iseed / 1988, 1989, 1990, 1991 /
116
118 k = m
119 n2 = m+n
120 IF( n.GT.0 ) THEN
121 np1 = m+1
122 ELSE
123 np1 = 1
124 END IF
125 lwork = n2*n2*nb
126
127
128
129 ALLOCATE(a(m,n2),af(m,n2),q(n2,n2),r(n2,n2),rwork(n2),
130 $ work(lwork),t(nb,m),c(n2,m),cf(n2,m),
131 $ d(m,n2),df(m,n2) )
132
133
134
135 ldt=nb
136 CALL slaset(
'Full', m, n2, zero, zero, a, m )
137 CALL slaset(
'Full', nb, m, zero, zero, t, nb )
138 DO j=1,m
139 CALL slarnv( 2, iseed, m-j+1, a( j, j ) )
140 END DO
141 IF( n.GT.0 ) THEN
142 DO j=1,n-l
143 CALL slarnv( 2, iseed, m, a( 1, min(n+m,m+1) + j - 1 ) )
144 END DO
145 END IF
146 IF( l.GT.0 ) THEN
147 DO j=1,l
148 CALL slarnv( 2, iseed, m-j+1, a( j, min(n+m,n+m-l+1)
149 $ + j - 1 ) )
150 END DO
151 END IF
152
153
154
155 CALL slacpy(
'Full', m, n2, a, m, af, m )
156
157
158
159 CALL stplqt( m,n,l,nb,af,m,af(1,np1),m,t,ldt,work,info)
160
161
162
163 CALL slaset(
'Full', n2, n2, zero, one, q, n2 )
164 CALL sgemlqt(
'L',
'N', n2, n2, k, nb, af, m, t, ldt, q, n2,
165 $ work, info )
166
167
168
169 CALL slaset(
'Full', n2, n2, zero, zero, r, n2 )
170 CALL slacpy(
'Lower', m, n2, af, m, r, n2 )
171
172
173
174 CALL sgemm(
'N',
'T', m, n2, n2, -one, a, m, q, n2, one, r, n2)
175 anorm =
slange(
'1', m, n2, a, m, rwork )
176 resid =
slange(
'1', m, n2, r, n2, rwork )
177 IF( anorm.GT.zero ) THEN
178 result( 1 ) = resid / (eps*anorm*max(1,n2))
179 ELSE
180 result( 1 ) = zero
181 END IF
182
183
184
185 CALL slaset(
'Full', n2, n2, zero, one, r, n2 )
186 CALL ssyrk(
'U',
'N', n2, n2, -one, q, n2, one, r, n2 )
187 resid =
slansy(
'1',
'Upper', n2, r, n2, rwork )
188 result( 2 ) = resid / (eps*max(1,n2))
189
190
191
192 CALL slaset(
'Full', n2, m, zero, one, c, n2 )
193 DO j=1,m
194 CALL slarnv( 2, iseed, n2, c( 1, j ) )
195 END DO
196 cnorm =
slange(
'1', n2, m, c, n2, rwork)
197 CALL slacpy(
'Full', n2, m, c, n2, cf, n2 )
198
199
200
201 CALL stpmlqt(
'L',
'N', n,m,k,l,nb,af(1, np1),m,t,ldt,cf,n2,
202 $ cf(np1,1),n2,work,info)
203
204
205
206 CALL sgemm(
'N',
'N', n2, m, n2, -one, q, n2, c, n2, one, cf, n2 )
207 resid =
slange(
'1', n2, m, cf, n2, rwork )
208 IF( cnorm.GT.zero ) THEN
209 result( 3 ) = resid / (eps*max(1,n2)*cnorm)
210 ELSE
211 result( 3 ) = zero
212 END IF
213
214
215
216
217 CALL slacpy(
'Full', n2, m, c, n2, cf, n2 )
218
219
220
221 CALL stpmlqt(
'L',
'T',n,m,k,l,nb,af(1,np1),m,t,ldt,cf,n2,
222 $ cf(np1,1),n2,work,info)
223
224
225
226 CALL sgemm(
'T',
'N',n2,m,n2,-one,q,n2,c,n2,one,cf,n2)
227 resid =
slange(
'1', n2, m, cf, n2, rwork )
228
229 IF( cnorm.GT.zero ) THEN
230 result( 4 ) = resid / (eps*max(1,n2)*cnorm)
231 ELSE
232 result( 4 ) = zero
233 END IF
234
235
236
237 DO j=1,n2
238 CALL slarnv( 2, iseed, m, d( 1, j ) )
239 END DO
240 dnorm =
slange(
'1', m, n2, d, m, rwork)
241 CALL slacpy(
'Full', m, n2, d, m, df, m )
242
243
244
245 CALL stpmlqt(
'R',
'N',m,n,k,l,nb,af(1,np1),m,t,ldt,df,m,
246 $ df(1,np1),m,work,info)
247
248
249
250 CALL sgemm(
'N',
'N',m,n2,n2,-one,d,m,q,n2,one,df,m)
251 resid =
slange(
'1',m, n2,df,m,rwork )
252 IF( cnorm.GT.zero ) THEN
253 result( 5 ) = resid / (eps*max(1,n2)*dnorm)
254 ELSE
255 result( 5 ) = zero
256 END IF
257
258
259
260 CALL slacpy(
'Full',m,n2,d,m,df,m )
261
262
263
264 CALL stpmlqt(
'R',
'T',m,n,k,l,nb,af(1,np1),m,t,ldt,df,m,
265 $ df(1,np1),m,work,info)
266
267
268
269
270 CALL sgemm(
'N',
'T', m, n2, n2, -one, d, m, q, n2, one, df, m )
271 resid =
slange(
'1', m, n2, df, m, rwork )
272 IF( cnorm.GT.zero ) THEN
273 result( 6 ) = resid / (eps*max(1,n2)*dnorm)
274 ELSE
275 result( 6 ) = zero
276 END IF
277
278
279
280 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
281 RETURN
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
subroutine stplqt(m, n, l, mb, a, lda, b, ldb, t, ldt, work, info)
STPLQT
subroutine stpmlqt(side, trans, m, n, k, l, mb, v, ldv, t, ldt, a, lda, b, ldb, work, info)
STPMLQT