85 SUBROUTINE stsqr01(TSSW, M, N, MB, NB, RESULT)
103 REAL,
ALLOCATABLE :: AF(:,:), Q(:,:),
104 $ r(:,:), rwork(:), work( : ), t(:),
105 $ cf(:,:), df(:,:), a(:,:), c(:,:), d(:,:), lq(:,:)
109 parameter( zero = 0.0, one = 1.0 )
112 LOGICAL TESTZEROS, TS
113 INTEGER INFO, J, K, L, LWORK, TSIZE, MNB
114 REAL ANORM, EPS, RESID, CNORM, DNORM
118 REAL TQUERY( 5 ), WORKQUERY
121 REAL SLAMCH, SLANGE, SLANSY
124 EXTERNAL slamch,
slarnv, slange, slansy, lsame, ilaenv
132 COMMON / srnamc / srnamt
135 DATA iseed / 1988, 1989, 1990, 1991 /
139 ts = lsame(tssw,
'TS')
145 eps = slamch(
'Epsilon' )
153 ALLOCATE ( a(m,n), af(m,n), q(l,l), r(m,l), rwork(l),
155 $ d(n,m), df(n,m), lq(l,n) )
160 CALL slarnv( 2, iseed, m, a( 1, j ) )
165 CALL slarnv( 2, iseed, m/2, a( m/4, j ) )
169 CALL slacpy(
'Full', m, n, a, m, af, m )
175 CALL sgeqr( m, n, af, m, tquery, -1, workquery, -1, info )
176 tsize = int( tquery( 1 ) )
177 lwork = int( workquery )
178 CALL sgemqr(
'L',
'N', m, m, k, af, m, tquery, tsize, cf, m,
179 $ workquery, -1, info)
180 lwork = max( lwork, int( workquery ) )
181 CALL sgemqr(
'L',
'N', m, n, k, af, m, tquery, tsize, cf, m,
182 $ workquery, -1, info)
183 lwork = max( lwork, int( workquery ) )
184 CALL sgemqr(
'L',
'T', m, n, k, af, m, tquery, tsize, cf, m,
185 $ workquery, -1, info)
186 lwork = max( lwork, int( workquery ) )
187 CALL sgemqr(
'R',
'N', n, m, k, af, m, tquery, tsize, df, n,
188 $ workquery, -1, info)
189 lwork = max( lwork, int( workquery ) )
190 CALL sgemqr(
'R',
'T', n, m, k, af, m, tquery, tsize, df, n,
191 $ workquery, -1, info)
192 lwork = max( lwork, int( workquery ) )
193 ALLOCATE ( t( tsize ) )
194 ALLOCATE ( work( lwork ) )
196 CALL sgeqr( m, n, af, m, t, tsize, work, lwork, info )
200 CALL slaset(
'Full', m, m, zero, one, q, m )
202 CALL sgemqr(
'L',
'N', m, m, k, af, m, t, tsize, q, m,
203 $ work, lwork, info )
207 CALL slaset(
'Full', m, n, zero, zero, r, m )
208 CALL slacpy(
'Upper', m, n, af, m, r, m )
212 CALL sgemm(
'T',
'N', m, n, m, -one, q, m, a, m, one, r, m )
213 anorm = slange(
'1', m, n, a, m, rwork )
214 resid = slange(
'1', m, n, r, m, rwork )
215 IF( anorm.GT.zero )
THEN
216 result( 1 ) = resid / (eps*max(1,m)*anorm)
223 CALL slaset(
'Full', m, m, zero, one, r, m )
224 CALL ssyrk(
'U',
'C', m, m, -one, q, m, one, r, m )
225 resid = slansy(
'1',
'Upper', m, r, m, rwork )
226 result( 2 ) = resid / (eps*max(1,m))
231 CALL slarnv( 2, iseed, m, c( 1, j ) )
233 cnorm = slange(
'1', m, n, c, m, rwork)
234 CALL slacpy(
'Full', m, n, c, m, cf, m )
239 CALL sgemqr(
'L',
'N', m, n, k, af, m, t, tsize, cf, m,
244 CALL sgemm(
'N',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
245 resid = slange(
'1', m, n, cf, m, rwork )
246 IF( cnorm.GT.zero )
THEN
247 result( 3 ) = resid / (eps*max(1,m)*cnorm)
254 CALL slacpy(
'Full', m, n, c, m, cf, m )
259 CALL sgemqr(
'L',
'T', m, n, k, af, m, t, tsize, cf, m,
264 CALL sgemm(
'T',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
265 resid = slange(
'1', m, n, cf, m, rwork )
266 IF( cnorm.GT.zero )
THEN
267 result( 4 ) = resid / (eps*max(1,m)*cnorm)
275 CALL slarnv( 2, iseed, n, d( 1, j ) )
277 dnorm = slange(
'1', n, m, d, n, rwork)
278 CALL slacpy(
'Full', n, m, d, n, df, n )
283 CALL sgemqr(
'R',
'N', n, m, k, af, m, t, tsize, df, n,
288 CALL sgemm(
'N',
'N', n, m, m, -one, d, n, q, m, one, df, n )
289 resid = slange(
'1', n, m, df, n, rwork )
290 IF( dnorm.GT.zero )
THEN
291 result( 5 ) = resid / (eps*max(1,m)*dnorm)
298 CALL slacpy(
'Full', n, m, d, n, df, n )
302 CALL sgemqr(
'R',
'T', n, m, k, af, m, t, tsize, df, n,
307 CALL sgemm(
'N',
'T', n, m, m, -one, d, n, q, m, one, df, n )
308 resid = slange(
'1', n, m, df, n, rwork )
309 IF( cnorm.GT.zero )
THEN
310 result( 6 ) = resid / (eps*max(1,m)*dnorm)
318 CALL sgelq( m, n, af, m, tquery, -1, workquery, -1, info )
319 tsize = int( tquery( 1 ) )
320 lwork = int( workquery )
321 CALL sgemlq(
'R',
'N', n, n, k, af, m, tquery, tsize, q, n,
322 $ workquery, -1, info )
323 lwork = max( lwork, int( workquery ) )
324 CALL sgemlq(
'L',
'N', n, m, k, af, m, tquery, tsize, df, n,
325 $ workquery, -1, info)
326 lwork = max( lwork, int( workquery ) )
327 CALL sgemlq(
'L',
'T', n, m, k, af, m, tquery, tsize, df, n,
328 $ workquery, -1, info)
329 lwork = max( lwork, int( workquery ) )
330 CALL sgemlq(
'R',
'N', m, n, k, af, m, tquery, tsize, cf, m,
331 $ workquery, -1, info)
332 lwork = max( lwork, int( workquery ) )
333 CALL sgemlq(
'R',
'T', m, n, k, af, m, tquery, tsize, cf, m,
334 $ workquery, -1, info)
335 lwork = max( lwork, int( workquery ) )
336 ALLOCATE ( t( tsize ) )
337 ALLOCATE ( work( lwork ) )
339 CALL sgelq( m, n, af, m, t, tsize, work, lwork, info )
344 CALL slaset(
'Full', n, n, zero, one, q, n )
346 CALL sgemlq(
'R',
'N', n, n, k, af, m, t, tsize, q, n,
347 $ work, lwork, info )
351 CALL slaset(
'Full', m, n, zero, zero, lq, l )
352 CALL slacpy(
'Lower', m, n, af, m, lq, l )
356 CALL sgemm(
'N',
'T', m, n, n, -one, a, m, q, n, one, lq, l )
357 anorm = slange(
'1', m, n, a, m, rwork )
358 resid = slange(
'1', m, n, lq, l, rwork )
359 IF( anorm.GT.zero )
THEN
360 result( 1 ) = resid / (eps*max(1,n)*anorm)
367 CALL slaset(
'Full', n, n, zero, one, lq, l )
368 CALL ssyrk(
'U',
'C', n, n, -one, q, n, one, lq, l )
369 resid = slansy(
'1',
'Upper', n, lq, l, rwork )
370 result( 2 ) = resid / (eps*max(1,n))
375 CALL slarnv( 2, iseed, n, d( 1, j ) )
377 dnorm = slange(
'1', n, m, d, n, rwork)
378 CALL slacpy(
'Full', n, m, d, n, df, n )
382 CALL sgemlq(
'L',
'N', n, m, k, af, m, t, tsize, df, n,
387 CALL sgemm(
'N',
'N', n, m, n, -one, q, n, d, n, one, df, n )
388 resid = slange(
'1', n, m, df, n, rwork )
389 IF( dnorm.GT.zero )
THEN
390 result( 3 ) = resid / (eps*max(1,n)*dnorm)
397 CALL slacpy(
'Full', n, m, d, n, df, n )
401 CALL sgemlq(
'L',
'T', n, m, k, af, m, t, tsize, df, n,
406 CALL sgemm(
'T',
'N', n, m, n, -one, q, n, d, n, one, df, n )
407 resid = slange(
'1', n, m, df, n, rwork )
408 IF( dnorm.GT.zero )
THEN
409 result( 4 ) = resid / (eps*max(1,n)*dnorm)
417 CALL slarnv( 2, iseed, m, c( 1, j ) )
419 cnorm = slange(
'1', m, n, c, m, rwork)
420 CALL slacpy(
'Full', m, n, c, m, cf, m )
424 CALL sgemlq(
'R',
'N', m, n, k, af, m, t, tsize, cf, m,
429 CALL sgemm(
'N',
'N', m, n, n, -one, c, m, q, n, one, cf, m )
430 resid = slange(
'1', n, m, df, n, rwork )
431 IF( cnorm.GT.zero )
THEN
432 result( 5 ) = resid / (eps*max(1,n)*cnorm)
439 CALL slacpy(
'Full', m, n, c, m, cf, m )
443 CALL sgemlq(
'R',
'T', m, n, k, af, m, t, tsize, cf, m,
448 CALL sgemm(
'N',
'T', m, n, n, -one, c, m, q, n, one, cf, m )
449 resid = slange(
'1', m, n, cf, m, rwork )
450 IF( cnorm.GT.zero )
THEN
451 result( 6 ) = resid / (eps*max(1,n)*cnorm)
460 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
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 sgeqr(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
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 sgemqr(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
subroutine sgemlq(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
subroutine sgelq(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine stsqr01(TSSW, M, N, MB, NB, RESULT)
STSQR01