OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
stsqr01.f
Go to the documentation of this file.
1*> \brief \b STSQR01
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE STSQR01(TSSW, M,N, MB, NB, RESULT)
12*
13* .. Scalar Arguments ..
14* INTEGER M, N, MB
15* .. Return values ..
16* REAL RESULT(6)
17*
18*
19*> \par Purpose:
20* =============
21*>
22*> \verbatim
23*>
24*> DTSQR01 tests DGEQR , DGELQ, DGEMLQ and DGEMQR.
25*> \endverbatim
26*
27* Arguments:
28* ==========
29*
30*> \param[in] TSSW
31*> \verbatim
32*> TSSW is CHARACTER
33*> 'TS' for testing tall skinny QR
34*> and anything else for testing short wide LQ
35*> \endverbatim
36*> \param[in] M
37*> \verbatim
38*> M is INTEGER
39*> Number of rows in test matrix.
40*> \endverbatim
41*>
42*> \param[in] N
43*> \verbatim
44*> N is INTEGER
45*> Number of columns in test matrix.
46*> \endverbatim
47*> \param[in] MB
48*> \verbatim
49*> MB is INTEGER
50*> Number of row in row block in test matrix.
51*> \endverbatim
52*>
53*> \param[in] NB
54*> \verbatim
55*> NB is INTEGER
56*> Number of columns in column block test matrix.
57*> \endverbatim
58*>
59*> \param[out] RESULT
60*> \verbatim
61*> RESULT is REAL array, dimension (6)
62*> Results of each of the six tests below.
63*>
64*> RESULT(1) = | A - Q R | or | A - L Q |
65*> RESULT(2) = | I - Q^H Q | or | I - Q Q^H |
66*> RESULT(3) = | Q C - Q C |
67*> RESULT(4) = | Q^H C - Q^H C |
68*> RESULT(5) = | C Q - C Q |
69*> RESULT(6) = | C Q^H - C Q^H |
70*> \endverbatim
71*
72* Authors:
73* ========
74*
75*> \author Univ. of Tennessee
76*> \author Univ. of California Berkeley
77*> \author Univ. of Colorado Denver
78*> \author NAG Ltd.
79*
80*> \ingroup double_lin
81*
82* =====================================================================
83 SUBROUTINE stsqr01(TSSW, M, N, MB, NB, RESULT)
84 IMPLICIT NONE
85*
86* -- LAPACK test routine --
87* -- LAPACK is a software package provided by Univ. of Tennessee, --
88* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
89*
90* .. Scalar Arguments ..
91 CHARACTER TSSW
92 INTEGER M, N, MB, NB
93* .. Return values ..
94 REAL RESULT(6)
95*
96* =====================================================================
97*
98* ..
99* .. Local allocatable arrays
100 REAL, ALLOCATABLE :: AF(:,:), Q(:,:),
101 $ R(:,:), RWORK(:), WORK( : ), T(:),
102 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:), LQ(:,:)
103*
104* .. Parameters ..
105 REAL ONE, ZERO
106 parameter( zero = 0.0, one = 1.0 )
107* ..
108* .. Local Scalars ..
109 LOGICAL TESTZEROS, TS
110 INTEGER INFO, J, K, L, LWORK, TSIZE, MNB
111 REAL ANORM, EPS, RESID, CNORM, DNORM
112* ..
113* .. Local Arrays ..
114 INTEGER ISEED( 4 )
115 REAL TQUERY( 5 ), WORKQUERY( 1 )
116* ..
117* .. External Functions ..
118 REAL SLAMCH, SLANGE, SLANSY
119 LOGICAL LSAME
120 INTEGER ILAENV
121 EXTERNAL slamch, slarnv, slange, slansy, lsame, ilaenv
122* ..
123* .. Intrinsic Functions ..
124 INTRINSIC max, min
125* .. Scalars in Common ..
126 CHARACTER*32 srnamt
127* ..
128* .. Common blocks ..
129 COMMON / srnamc / srnamt
130* ..
131* .. Data statements ..
132 DATA iseed / 1988, 1989, 1990, 1991 /
133*
134* TEST TALL SKINNY OR SHORT WIDE
135*
136 ts = lsame(tssw, 'TS')
137*
138* TEST MATRICES WITH HALF OF MATRIX BEING ZEROS
139*
140 testzeros = .false.
141*
142 eps = slamch( 'Epsilon' )
143 k = min(m,n)
144 l = max(m,n,1)
145 mnb = max( mb, nb)
146 lwork = max(3,l)*mnb
147*
148* Dynamically allocate local arrays
149*
150 ALLOCATE ( a(m,n), af(m,n), q(l,l), r(m,l), rwork(l),
151 $ c(m,n), cf(m,n),
152 $ d(n,m), df(n,m), lq(l,n) )
153*
154* Put random numbers into A and copy to AF
155*
156 DO j=1,n
157 CALL slarnv( 2, iseed, m, a( 1, j ) )
158 END DO
159 IF (testzeros) THEN
160 IF (m.GE.4) THEN
161 DO j=1,n
162 CALL slarnv( 2, iseed, m/2, a( m/4, j ) )
163 END DO
164 END IF
165 END IF
166 CALL slacpy( 'Full', m, n, a, m, af, m )
167*
168 IF (ts) THEN
169*
170* Factor the matrix A in the array AF.
171*
172 CALL sgeqr( m, n, af, m, tquery, -1, workquery, -1, info )
173 tsize = int( tquery( 1 ) )
174 lwork = int( workquery( 1 ) )
175 CALL sgemqr( 'L', 'N', m, m, k, af, m, tquery, tsize, cf, m,
176 $ workquery, -1, info)
177 lwork = max( lwork, int( workquery( 1 ) ) )
178 CALL sgemqr( 'L', 'N', m, n, k, af, m, tquery, tsize, cf, m,
179 $ workquery, -1, info)
180 lwork = max( lwork, int( workquery( 1 ) ) )
181 CALL sgemqr( 'L', 'T', m, n, k, af, m, tquery, tsize, cf, m,
182 $ workquery, -1, info)
183 lwork = max( lwork, int( workquery( 1 ) ) )
184 CALL sgemqr( 'R', 'N', n, m, k, af, m, tquery, tsize, df, n,
185 $ workquery, -1, info)
186 lwork = max( lwork, int( workquery( 1 ) ) )
187 CALL sgemqr( 'R', 'T', n, m, k, af, m, tquery, tsize, df, n,
188 $ workquery, -1, info)
189 lwork = max( lwork, int( workquery( 1 ) ) )
190 ALLOCATE ( t( tsize ) )
191 ALLOCATE ( work( lwork ) )
192 srnamt = 'SGEQR'
193 CALL sgeqr( m, n, af, m, t, tsize, work, lwork, info )
194*
195* Generate the m-by-m matrix Q
196*
197 CALL slaset( 'Full', m, m, zero, one, q, m )
198 srnamt = 'SGEMQR'
199 CALL sgemqr( 'L', 'N', m, m, k, af, m, t, tsize, q, m,
200 $ work, lwork, info )
201*
202* Copy R
203*
204 CALL slaset( 'Full', m, n, zero, zero, r, m )
205 CALL slacpy( 'Upper', m, n, af, m, r, m )
206*
207* Compute |R - Q'*A| / |A| and store in RESULT(1)
208*
209 CALL sgemm( 'T', 'N', m, n, m, -one, q, m, a, m, one, r, m )
210 anorm = slange( '1', m, n, a, m, rwork )
211 resid = slange( '1', m, n, r, m, rwork )
212 IF( anorm.GT.zero ) THEN
213 result( 1 ) = resid / (eps*max(1,m)*anorm)
214 ELSE
215 result( 1 ) = zero
216 END IF
217*
218* Compute |I - Q'*Q| and store in RESULT(2)
219*
220 CALL slaset( 'Full', m, m, zero, one, r, m )
221 CALL ssyrk( 'U', 'C', m, m, -one, q, m, one, r, m )
222 resid = slansy( '1', 'Upper', m, r, m, rwork )
223 result( 2 ) = resid / (eps*max(1,m))
224*
225* Generate random m-by-n matrix C and a copy CF
226*
227 DO j=1,n
228 CALL slarnv( 2, iseed, m, c( 1, j ) )
229 END DO
230 cnorm = slange( '1', m, n, c, m, rwork)
231 CALL slacpy( 'Full', m, n, c, m, cf, m )
232*
233* Apply Q to C as Q*C
234*
235 srnamt = 'DGEQR'
236 CALL sgemqr( 'l', 'n', M, N, K, AF, M, T, TSIZE, CF, M,
237 $ WORK, LWORK, INFO)
238*
239* Compute |Q*C - Q*C| / |C|
240*
241 CALL SGEMM( 'n', 'n', M, N, M, -ONE, Q, M, C, M, ONE, CF, M )
242 RESID = SLANGE( '1', M, N, CF, M, RWORK )
243.GT. IF( CNORMZERO ) THEN
244 RESULT( 3 ) = RESID / (EPS*MAX(1,M)*CNORM)
245 ELSE
246 RESULT( 3 ) = ZERO
247 END IF
248*
249* Copy C into CF again
250*
251 CALL SLACPY( 'full', M, N, C, M, CF, M )
252*
253* Apply Q to C as QT*C
254*
255 srnamt = 'dgeqr'
256 CALL SGEMQR( 'l', 't', M, N, K, AF, M, T, TSIZE, CF, M,
257 $ WORK, LWORK, INFO)
258*
259* Compute |QT*C - QT*C| / |C|
260*
261 CALL SGEMM( 't', 'n', M, N, M, -ONE, Q, M, C, M, ONE, CF, M )
262 RESID = SLANGE( '1', M, N, CF, M, RWORK )
263.GT. IF( CNORMZERO ) THEN
264 RESULT( 4 ) = RESID / (EPS*MAX(1,M)*CNORM)
265 ELSE
266 RESULT( 4 ) = ZERO
267 END IF
268*
269* Generate random n-by-m matrix D and a copy DF
270*
271 DO J=1,M
272 CALL SLARNV( 2, ISEED, N, D( 1, J ) )
273 END DO
274 DNORM = SLANGE( '1', N, M, D, N, RWORK)
275 CALL SLACPY( 'full', N, M, D, N, DF, N )
276*
277* Apply Q to D as D*Q
278*
279 srnamt = 'dgeqr'
280 CALL SGEMQR( 'r', 'n', N, M, K, AF, M, T, TSIZE, DF, N,
281 $ WORK, LWORK, INFO)
282*
283* Compute |D*Q - D*Q| / |D|
284*
285 CALL SGEMM( 'n', 'n', N, M, M, -ONE, D, N, Q, M, ONE, DF, N )
286 RESID = SLANGE( '1', N, M, DF, N, RWORK )
287.GT. IF( DNORMZERO ) THEN
288 RESULT( 5 ) = RESID / (EPS*MAX(1,M)*DNORM)
289 ELSE
290 RESULT( 5 ) = ZERO
291 END IF
292*
293* Copy D into DF again
294*
295 CALL SLACPY( 'full', N, M, D, N, DF, N )
296*
297* Apply Q to D as D*QT
298*
299 CALL SGEMQR( 'r', 't', N, M, K, AF, M, T, TSIZE, DF, N,
300 $ WORK, LWORK, INFO)
301*
302* Compute |D*QT - D*QT| / |D|
303*
304 CALL SGEMM( 'n', 't', N, M, M, -ONE, D, N, Q, M, ONE, DF, N )
305 RESID = SLANGE( '1', N, M, DF, N, RWORK )
306.GT. IF( CNORMZERO ) THEN
307 RESULT( 6 ) = RESID / (EPS*MAX(1,M)*DNORM)
308 ELSE
309 RESULT( 6 ) = ZERO
310 END IF
311*
312* Short and wide
313*
314 ELSE
315 CALL SGELQ( M, N, AF, M, TQUERY, -1, WORKQUERY, -1, INFO )
316 TSIZE = INT( TQUERY( 1 ) )
317 LWORK = INT( WORKQUERY( 1 ))
318 CALL SGEMLQ( 'r', 'n', N, N, K, AF, M, TQUERY, TSIZE, Q, N,
319 $ WORKQUERY, -1, INFO )
320 LWORK = MAX( LWORK, INT( WORKQUERY( 1 ) ) )
321 CALL SGEMLQ( 'l', 'n', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
322 $ WORKQUERY, -1, INFO)
323 LWORK = MAX( LWORK, INT( WORKQUERY( 1 ) ) )
324 CALL SGEMLQ( 'l', 't', N, M, K, AF, M, TQUERY, TSIZE, DF, N,
325 $ WORKQUERY, -1, INFO)
326 LWORK = MAX( LWORK, INT( WORKQUERY( 1 ) ) )
327 CALL SGEMLQ( 'r', 'n', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
328 $ WORKQUERY, -1, INFO)
329 LWORK = MAX( LWORK, INT( WORKQUERY( 1 ) ) )
330 CALL SGEMLQ( 'r', 't', M, N, K, AF, M, TQUERY, TSIZE, CF, M,
331 $ WORKQUERY, -1, INFO)
332 LWORK = MAX( LWORK, INT( WORKQUERY( 1 ) ) )
333 ALLOCATE ( T( TSIZE ) )
334 ALLOCATE ( WORK( LWORK ) )
335 srnamt = 'sgelq'
336 CALL SGELQ( M, N, AF, M, T, TSIZE, WORK, LWORK, INFO )
337*
338*
339* Generate the n-by-n matrix Q
340*
341 CALL SLASET( 'full', N, N, ZERO, ONE, Q, N )
342 srnamt = 'sgemlq'
343 CALL SGEMLQ( 'r', 'n', N, N, K, AF, M, T, TSIZE, Q, N,
344 $ WORK, LWORK, INFO )
345*
346* Copy R
347*
348 CALL SLASET( 'full', M, N, ZERO, ZERO, LQ, L )
349 CALL SLACPY( 'lower', M, N, AF, M, LQ, L )
350*
351* Compute |L - A*Q'| / |A| and store in RESULT(1)
352*
353 CALL SGEMM( 'n', 't', M, N, N, -ONE, A, M, Q, N, ONE, LQ, L )
354 ANORM = SLANGE( '1', M, N, A, M, RWORK )
355 RESID = SLANGE( '1', M, N, LQ, L, RWORK )
356.GT. IF( ANORMZERO ) THEN
357 RESULT( 1 ) = RESID / (EPS*MAX(1,N)*ANORM)
358 ELSE
359 RESULT( 1 ) = ZERO
360 END IF
361*
362* Compute |I - Q'*Q| and store in RESULT(2)
363*
364 CALL SLASET( 'full', N, N, ZERO, ONE, LQ, L )
365 CALL SSYRK( 'u', 'c', N, N, -ONE, Q, N, ONE, LQ, L )
366 RESID = SLANSY( '1', 'upper', N, LQ, L, RWORK )
367 RESULT( 2 ) = RESID / (EPS*MAX(1,N))
368*
369* Generate random m-by-n matrix C and a copy CF
370*
371 DO J=1,M
372 CALL SLARNV( 2, ISEED, N, D( 1, J ) )
373 END DO
374 DNORM = SLANGE( '1', N, M, D, N, RWORK)
375 CALL SLACPY( 'full', N, M, D, N, DF, N )
376*
377* Apply Q to C as Q*C
378*
379 CALL SGEMLQ( 'l', 'n', N, M, K, AF, M, T, TSIZE, DF, N,
380 $ WORK, LWORK, INFO)
381*
382* Compute |Q*D - Q*D| / |D|
383*
384 CALL SGEMM( 'n', 'n', N, M, N, -ONE, Q, N, D, N, ONE, DF, N )
385 RESID = SLANGE( '1', N, M, DF, N, RWORK )
386.GT. IF( DNORMZERO ) THEN
387 RESULT( 3 ) = RESID / (EPS*MAX(1,N)*DNORM)
388 ELSE
389 RESULT( 3 ) = ZERO
390 END IF
391*
392* Copy D into DF again
393*
394 CALL SLACPY( 'full', N, M, D, N, DF, N )
395*
396* Apply Q to D as QT*D
397*
398 CALL SGEMLQ( 'l', 't', N, M, K, AF, M, T, TSIZE, DF, N,
399 $ WORK, LWORK, INFO)
400*
401* Compute |QT*D - QT*D| / |D|
402*
403 CALL SGEMM( 't', 'n', N, M, N, -ONE, Q, N, D, N, ONE, DF, N )
404 RESID = SLANGE( '1', N, M, DF, N, RWORK )
405.GT. IF( DNORMZERO ) THEN
406 RESULT( 4 ) = RESID / (EPS*MAX(1,N)*DNORM)
407 ELSE
408 RESULT( 4 ) = ZERO
409 END IF
410*
411* Generate random n-by-m matrix D and a copy DF
412*
413 DO J=1,N
414 CALL SLARNV( 2, ISEED, M, C( 1, J ) )
415 END DO
416 CNORM = SLANGE( '1', M, N, C, M, RWORK)
417 CALL SLACPY( 'full', M, N, C, M, CF, M )
418*
419* Apply Q to C as C*Q
420*
421 CALL SGEMLQ( 'r', 'n', M, N, K, AF, M, T, TSIZE, CF, M,
422 $ WORK, LWORK, INFO)
423*
424* Compute |C*Q - C*Q| / |C|
425*
426 CALL SGEMM( 'n', 'n', M, N, N, -ONE, C, M, Q, N, ONE, CF, M )
427 RESID = SLANGE( '1', N, M, DF, N, RWORK )
428.GT. IF( CNORMZERO ) THEN
429 RESULT( 5 ) = RESID / (EPS*MAX(1,N)*CNORM)
430 ELSE
431 RESULT( 5 ) = ZERO
432 END IF
433*
434* Copy C into CF again
435*
436 CALL SLACPY( 'full', M, N, C, M, CF, M )
437*
438* Apply Q to D as D*QT
439*
440 CALL SGEMLQ( 'r', 't', M, N, K, AF, M, T, TSIZE, CF, M,
441 $ WORK, LWORK, INFO)
442*
443* Compute |C*QT - C*QT| / |C|
444*
445 CALL SGEMM( 'n', 't', M, N, N, -ONE, C, M, Q, N, ONE, CF, M )
446 RESID = SLANGE( '1', M, N, CF, M, RWORK )
447.GT. IF( CNORMZERO ) THEN
448 RESULT( 6 ) = RESID / (EPS*MAX(1,N)*CNORM)
449 ELSE
450 RESULT( 6 ) = ZERO
451 END IF
452*
453 END IF
454*
455* Deallocate all arrays
456*
457 DEALLOCATE ( A, AF, Q, R, RWORK, WORK, T, C, D, CF, DF)
458*
459 RETURN
460 END
dgeqr
subroutine dgeqr(m, n, a, lda, t, tsize, work, lwork, info)
DGEQR
Definition dgeqr.f:174
slarnv
subroutine slarnv(idist, iseed, n, x)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition slarnv.f:97
slaset
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.
Definition slaset.f:110
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
stsqr01
subroutine stsqr01(tssw, m, n, mb, nb, result)
STSQR01
Definition stsqr01.f:84
ssyrk
subroutine ssyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
SSYRK
Definition ssyrk.f:169
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:187
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sgelq
subroutine sgelq(m, n, a, lda, t, tsize, work, lwork, info)
SGELQ
Definition sgelq.f:172
sgemlq
subroutine sgemlq(side, trans, m, n, k, a, lda, t, tsize, c, ldc, work, lwork, info)
SGEMLQ
Definition sgemlq.f:170
sgemqr
subroutine sgemqr(side, trans, m, n, k, a, lda, t, tsize, c, ldc, work, lwork, info)
SGEMQR
Definition sgemqr.f:172
sgeqr
subroutine sgeqr(m, n, a, lda, t, tsize, work, lwork, info)
SGEQR
Definition sgeqr.f:174