LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
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  IF( cnorm.GT.zero ) 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  IF( cnorm.GT.zero ) 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  IF( dnorm.GT.zero ) 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  IF( cnorm.GT.zero ) 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  IF( anorm.GT.zero ) 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  IF( dnorm.GT.zero ) 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  IF( dnorm.GT.zero ) 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  IF( cnorm.GT.zero ) 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  IF( cnorm.GT.zero ) 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
subroutine slarnv(IDIST, ISEED, N, X)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: slarnv.f:97
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
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
subroutine stsqr01(TSSW, M, N, MB, NB, RESULT)
STSQR01
Definition: stsqr01.f:84
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
Definition: ssyrk.f:169
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
subroutine sgelq(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
SGELQ
Definition: sgelq.f:172
subroutine sgemlq(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
SGEMLQ
Definition: sgemlq.f:170
subroutine sgemqr(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
SGEMQR
Definition: sgemqr.f:172
subroutine sgeqr(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
SGEQR
Definition: sgeqr.f:174