LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
ctsqr01.f
Go to the documentation of this file.
1 *> \brief \b CTSQR01
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 CTSQR01(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 *> \date April 2012
81 *
82 * =====================================================================
83  SUBROUTINE ctsqr01(TSSW, M, N, MB, NB, RESULT)
84  IMPLICIT NONE
85 *
86 * -- LAPACK test routine (version 3.7.0) --
87 * -- LAPACK is a software package provided by Univ. of Tennessee, --
88 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
89 * April 2012
90 *
91 * .. Scalar Arguments ..
92  CHARACTER TSSW
93  INTEGER M, N, MB, NB
94 * .. Return values ..
95  REAL RESULT(6)
96 *
97 * =====================================================================
98 *
99 * ..
100 * .. Local allocatable arrays
101  COMPLEX, ALLOCATABLE :: AF(:,:), Q(:,:),
102  $ R(:,:), RWORK(:), WORK( : ), T(:),
103  $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:), LQ(:,:)
104 *
105 * .. Parameters ..
106  REAL ZERO
107  COMPLEX ONE, CZERO
108  parameter( zero = 0.0, one = (1.0,0.0), czero=(0.0,0.0) )
109 * ..
110 * .. Local Scalars ..
111  LOGICAL TESTZEROS, TS
112  INTEGER INFO, J, K, L, LWORK, TSIZE, MNB
113  REAL ANORM, EPS, RESID, CNORM, DNORM
114 * ..
115 * .. Local Arrays ..
116  INTEGER ISEED( 4 )
117  COMPLEX TQUERY( 5 ), WORKQUERY( 1 )
118 * ..
119 * .. External Functions ..
120  REAL SLAMCH, CLANGE, CLANSY
121  LOGICAL LSAME
122  INTEGER ILAENV
123  EXTERNAL slamch, clange, clansy, lsame, ilaenv
124 * ..
125 * .. Intrinsic Functions ..
126  INTRINSIC max, min
127 * .. Scalars in Common ..
128  CHARACTER*32 srnamt
129 * ..
130 * .. Common blocks ..
131  COMMON / srnamc / srnamt
132 * ..
133 * .. Data statements ..
134  DATA iseed / 1988, 1989, 1990, 1991 /
135 *
136 * TEST TALL SKINNY OR SHORT WIDE
137 *
138  ts = lsame(tssw, 'TS')
139 *
140 * TEST MATRICES WITH HALF OF MATRIX BEING ZEROS
141 *
142  testzeros = .false.
143 *
144  eps = slamch( 'Epsilon' )
145  k = min(m,n)
146  l = max(m,n,1)
147  mnb = max( mb, nb)
148  lwork = max(3,l)*mnb
149 *
150 * Dynamically allocate local arrays
151 *
152  ALLOCATE ( a(m,n), af(m,n), q(l,l), r(m,l), rwork(l),
153  $ c(m,n), cf(m,n),
154  $ d(n,m), df(n,m), lq(l,n) )
155 *
156 * Put random numbers into A and copy to AF
157 *
158  DO j=1,n
159  CALL clarnv( 2, iseed, m, a( 1, j ) )
160  END DO
161  IF (testzeros) THEN
162  IF (m.GE.4) THEN
163  DO j=1,n
164  CALL clarnv( 2, iseed, m/2, a( m/4, j ) )
165  END DO
166  END IF
167  END IF
168  CALL clacpy( 'Full', m, n, a, m, af, m )
169 *
170  IF (ts) THEN
171 *
172 * Factor the matrix A in the array AF.
173 *
174  CALL cgeqr( m, n, af, m, tquery, -1, workquery, -1, info )
175  tsize = int( tquery( 1 ) )
176  lwork = int( workquery( 1 ) )
177  CALL cgemqr( 'L', 'N', m, m, k, af, m, tquery, tsize, cf, m,
178  $ workquery, -1, info)
179  lwork = max( lwork, int( workquery( 1 ) ) )
180  CALL cgemqr( 'L', 'N', m, n, k, af, m, tquery, tsize, cf, m,
181  $ workquery, -1, info)
182  lwork = max( lwork, int( workquery( 1 ) ) )
183  CALL cgemqr( 'L', 'C', m, n, k, af, m, tquery, tsize, cf, m,
184  $ workquery, -1, info)
185  lwork = max( lwork, int( workquery( 1 ) ) )
186  CALL cgemqr( 'R', 'N', n, m, k, af, m, tquery, tsize, df, n,
187  $ workquery, -1, info)
188  lwork = max( lwork, int( workquery( 1 ) ) )
189  CALL cgemqr( 'R', 'C', n, m, k, af, m, tquery, tsize, df, n,
190  $ workquery, -1, info)
191  lwork = max( lwork, int( workquery( 1 ) ) )
192  ALLOCATE ( t( tsize ) )
193  ALLOCATE ( work( lwork ) )
194  srnamt = 'CGEQR'
195  CALL cgeqr( m, n, af, m, t, tsize, work, lwork, info )
196 *
197 * Generate the m-by-m matrix Q
198 *
199  CALL claset( 'Full', m, m, czero, one, q, m )
200  srnamt = 'CGEMQR'
201  CALL cgemqr( 'L', 'N', m, m, k, af, m, t, tsize, q, m,
202  $ work, lwork, info )
203 *
204 * Copy R
205 *
206  CALL claset( 'Full', m, n, czero, czero, r, m )
207  CALL clacpy( 'Upper', m, n, af, m, r, m )
208 *
209 * Compute |R - Q'*A| / |A| and store in RESULT(1)
210 *
211  CALL cgemm( 'C', 'N', m, n, m, -one, q, m, a, m, one, r, m )
212  anorm = clange( '1', m, n, a, m, rwork )
213  resid = clange( '1', m, n, r, m, rwork )
214  IF( anorm.GT.zero ) THEN
215  result( 1 ) = resid / (eps*max(1,m)*anorm)
216  ELSE
217  result( 1 ) = zero
218  END IF
219 *
220 * Compute |I - Q'*Q| and store in RESULT(2)
221 *
222  CALL claset( 'Full', m, m, czero, one, r, m )
223  CALL cherk( 'U', 'C', m, m, real(-one), q, m, real(one), r, m )
224  resid = clansy( '1', 'Upper', m, r, m, rwork )
225  result( 2 ) = resid / (eps*max(1,m))
226 *
227 * Generate random m-by-n matrix C and a copy CF
228 *
229  DO j=1,n
230  CALL clarnv( 2, iseed, m, c( 1, j ) )
231  END DO
232  cnorm = clange( '1', m, n, c, m, rwork)
233  CALL clacpy( 'Full', m, n, c, m, cf, m )
234 *
235 * Apply Q to C as Q*C
236 *
237  srnamt = 'CGEMQR'
238  CALL cgemqr( 'L', 'N', m, n, k, af, m, t, tsize, cf, m,
239  $ work, lwork, info)
240 *
241 * Compute |Q*C - Q*C| / |C|
242 *
243  CALL cgemm( 'N', 'N', m, n, m, -one, q, m, c, m, one, cf, m )
244  resid = clange( '1', m, n, cf, m, rwork )
245  IF( cnorm.GT.zero ) THEN
246  result( 3 ) = resid / (eps*max(1,m)*cnorm)
247  ELSE
248  result( 3 ) = zero
249  END IF
250 *
251 * Copy C into CF again
252 *
253  CALL clacpy( 'Full', m, n, c, m, cf, m )
254 *
255 * Apply Q to C as QT*C
256 *
257  srnamt = 'CGEMQR'
258  CALL cgemqr( 'L', 'C', m, n, k, af, m, t, tsize, cf, m,
259  $ work, lwork, info)
260 *
261 * Compute |QT*C - QT*C| / |C|
262 *
263  CALL cgemm( 'C', 'N', m, n, m, -one, q, m, c, m, one, cf, m )
264  resid = clange( '1', m, n, cf, m, rwork )
265  IF( cnorm.GT.zero ) THEN
266  result( 4 ) = resid / (eps*max(1,m)*cnorm)
267  ELSE
268  result( 4 ) = zero
269  END IF
270 *
271 * Generate random n-by-m matrix D and a copy DF
272 *
273  DO j=1,m
274  CALL clarnv( 2, iseed, n, d( 1, j ) )
275  END DO
276  dnorm = clange( '1', n, m, d, n, rwork)
277  CALL clacpy( 'Full', n, m, d, n, df, n )
278 *
279 * Apply Q to D as D*Q
280 *
281  srnamt = 'CGEMQR'
282  CALL cgemqr( 'R', 'N', n, m, k, af, m, t, tsize, df, n,
283  $ work, lwork, info)
284 *
285 * Compute |D*Q - D*Q| / |D|
286 *
287  CALL cgemm( 'N', 'N', n, m, m, -one, d, n, q, m, one, df, n )
288  resid = clange( '1', n, m, df, n, rwork )
289  IF( dnorm.GT.zero ) THEN
290  result( 5 ) = resid / (eps*max(1,m)*dnorm)
291  ELSE
292  result( 5 ) = zero
293  END IF
294 *
295 * Copy D into DF again
296 *
297  CALL clacpy( 'Full', n, m, d, n, df, n )
298 *
299 * Apply Q to D as D*QT
300 *
301  CALL cgemqr( 'R', 'C', n, m, k, af, m, t, tsize, df, n,
302  $ work, lwork, info)
303 *
304 * Compute |D*QT - D*QT| / |D|
305 *
306  CALL cgemm( 'N', 'C', n, m, m, -one, d, n, q, m, one, df, n )
307  resid = clange( '1', n, m, df, n, rwork )
308  IF( cnorm.GT.zero ) THEN
309  result( 6 ) = resid / (eps*max(1,m)*dnorm)
310  ELSE
311  result( 6 ) = zero
312  END IF
313 *
314 * Short and wide
315 *
316  ELSE
317  CALL cgelq( m, n, af, m, tquery, -1, workquery, -1, info )
318  tsize = int( tquery( 1 ) )
319  lwork = int( workquery( 1 ) )
320  CALL cgemlq( 'R', 'N', n, n, k, af, m, tquery, tsize, q, n,
321  $ workquery, -1, info )
322  lwork = max( lwork, int( workquery( 1 ) ) )
323  CALL cgemlq( 'L', 'N', n, m, k, af, m, tquery, tsize, df, n,
324  $ workquery, -1, info)
325  lwork = max( lwork, int( workquery( 1 ) ) )
326  CALL cgemlq( 'L', 'C', n, m, k, af, m, tquery, tsize, df, n,
327  $ workquery, -1, info)
328  lwork = max( lwork, int( workquery( 1 ) ) )
329  CALL cgemlq( 'R', 'N', m, n, k, af, m, tquery, tsize, cf, m,
330  $ workquery, -1, info)
331  lwork = max( lwork, int( workquery( 1 ) ) )
332  CALL cgemlq( 'R', 'C', m, n, k, af, m, tquery, tsize, cf, m,
333  $ workquery, -1, info)
334  lwork = max( lwork, int( workquery( 1 ) ) )
335  ALLOCATE ( t( tsize ) )
336  ALLOCATE ( work( lwork ) )
337  srnamt = 'CGELQ'
338  CALL cgelq( m, n, af, m, t, tsize, work, lwork, info )
339 *
340 *
341 * Generate the n-by-n matrix Q
342 *
343  CALL claset( 'Full', n, n, czero, one, q, n )
344  srnamt = 'CGEMLQ'
345  CALL cgemlq( 'R', 'N', n, n, k, af, m, t, tsize, q, n,
346  $ work, lwork, info )
347 *
348 * Copy R
349 *
350  CALL claset( 'Full', m, n, czero, czero, lq, l )
351  CALL clacpy( 'Lower', m, n, af, m, lq, l )
352 *
353 * Compute |L - A*Q'| / |A| and store in RESULT(1)
354 *
355  CALL cgemm( 'N', 'C', m, n, n, -one, a, m, q, n, one, lq, l )
356  anorm = clange( '1', m, n, a, m, rwork )
357  resid = clange( '1', m, n, lq, l, rwork )
358  IF( anorm.GT.zero ) THEN
359  result( 1 ) = resid / (eps*max(1,n)*anorm)
360  ELSE
361  result( 1 ) = zero
362  END IF
363 *
364 * Compute |I - Q'*Q| and store in RESULT(2)
365 *
366  CALL claset( 'Full', n, n, czero, one, lq, l )
367  CALL cherk( 'U', 'C', n, n, real(-one), q, n, real(one), lq, l)
368  resid = clansy( '1', 'Upper', n, lq, l, rwork )
369  result( 2 ) = resid / (eps*max(1,n))
370 *
371 * Generate random m-by-n matrix C and a copy CF
372 *
373  DO j=1,m
374  CALL clarnv( 2, iseed, n, d( 1, j ) )
375  END DO
376  dnorm = clange( '1', n, m, d, n, rwork)
377  CALL clacpy( 'Full', n, m, d, n, df, n )
378 *
379 * Apply Q to C as Q*C
380 *
381  CALL cgemlq( 'L', 'N', n, m, k, af, m, t, tsize, df, n,
382  $ work, lwork, info)
383 *
384 * Compute |Q*D - Q*D| / |D|
385 *
386  CALL cgemm( 'N', 'N', n, m, n, -one, q, n, d, n, one, df, n )
387  resid = clange( '1', n, m, df, n, rwork )
388  IF( dnorm.GT.zero ) THEN
389  result( 3 ) = resid / (eps*max(1,n)*dnorm)
390  ELSE
391  result( 3 ) = zero
392  END IF
393 *
394 * Copy D into DF again
395 *
396  CALL clacpy( 'Full', n, m, d, n, df, n )
397 *
398 * Apply Q to D as QT*D
399 *
400  CALL cgemlq( 'L', 'C', n, m, k, af, m, t, tsize, df, n,
401  $ work, lwork, info)
402 *
403 * Compute |QT*D - QT*D| / |D|
404 *
405  CALL cgemm( 'C', 'N', n, m, n, -one, q, n, d, n, one, df, n )
406  resid = clange( '1', n, m, df, n, rwork )
407  IF( dnorm.GT.zero ) THEN
408  result( 4 ) = resid / (eps*max(1,n)*dnorm)
409  ELSE
410  result( 4 ) = zero
411  END IF
412 *
413 * Generate random n-by-m matrix D and a copy DF
414 *
415  DO j=1,n
416  CALL clarnv( 2, iseed, m, c( 1, j ) )
417  END DO
418  cnorm = clange( '1', m, n, c, m, rwork)
419  CALL clacpy( 'Full', m, n, c, m, cf, m )
420 *
421 * Apply Q to C as C*Q
422 *
423  CALL cgemlq( 'R', 'N', m, n, k, af, m, t, tsize, cf, m,
424  $ work, lwork, info)
425 *
426 * Compute |C*Q - C*Q| / |C|
427 *
428  CALL cgemm( 'N', 'N', m, n, n, -one, c, m, q, n, one, cf, m )
429  resid = clange( '1', n, m, df, n, rwork )
430  IF( cnorm.GT.zero ) THEN
431  result( 5 ) = resid / (eps*max(1,n)*cnorm)
432  ELSE
433  result( 5 ) = zero
434  END IF
435 *
436 * Copy C into CF again
437 *
438  CALL clacpy( 'Full', m, n, c, m, cf, m )
439 *
440 * Apply Q to D as D*QT
441 *
442  CALL cgemlq( 'R', 'C', m, n, k, af, m, t, tsize, cf, m,
443  $ work, lwork, info)
444 *
445 * Compute |C*QT - C*QT| / |C|
446 *
447  CALL cgemm( 'N', 'C', m, n, n, -one, c, m, q, n, one, cf, m )
448  resid = clange( '1', m, n, cf, m, rwork )
449  IF( cnorm.GT.zero ) THEN
450  result( 6 ) = resid / (eps*max(1,n)*cnorm)
451  ELSE
452  result( 6 ) = zero
453  END IF
454 *
455  END IF
456 *
457 * Deallocate all arrays
458 *
459  DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
460 *
461  RETURN
462  END
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
cgelq
subroutine cgelq(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
CGELQ
Definition: cgelq.f:172
clacpy
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
cherk
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
Definition: cherk.f:175
cgeqr
subroutine cgeqr(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
CGEQR
Definition: cgeqr.f:174
cgemlq
subroutine cgemlq(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
CGEMLQ
Definition: cgemlq.f:170
claset
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:108
ctsqr01
subroutine ctsqr01(TSSW, M, N, MB, NB, RESULT)
CTSQR01
Definition: ctsqr01.f:84
cgemqr
subroutine cgemqr(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
CGEMQR
Definition: cgemqr.f:172
clarnv
subroutine clarnv(IDIST, ISEED, N, X)
CLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: clarnv.f:101