LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cunhr_col01.f
Go to the documentation of this file.
1*> \brief \b CUNHR_COL01
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 CUNHR_COL01( M, N, MB1, NB1, NB2, RESULT )
12*
13* .. Scalar Arguments ..
14* INTEGER M, N, MB1, NB1, NB2
15* .. Return values ..
16* DOUBLE PRECISION RESULT(6)
17*
18*
19*> \par Purpose:
20* =============
21*>
22*> \verbatim
23*>
24*> CUNHR_COL01 tests CUNGTSQR and CUNHR_COL using CLATSQR, CGEMQRT.
25*> Therefore, CLATSQR (part of CGEQR), CGEMQRT (part of CGEMQR)
26*> have to be tested before this test.
27*>
28*> \endverbatim
29*
30* Arguments:
31* ==========
32*
33*> \param[in] M
34*> \verbatim
35*> M is INTEGER
36*> Number of rows in test matrix.
37*> \endverbatim
38*> \param[in] N
39*> \verbatim
40*> N is INTEGER
41*> Number of columns in test matrix.
42*> \endverbatim
43*> \param[in] MB1
44*> \verbatim
45*> MB1 is INTEGER
46*> Number of row in row block in an input test matrix.
47*> \endverbatim
48*>
49*> \param[in] NB1
50*> \verbatim
51*> NB1 is INTEGER
52*> Number of columns in column block an input test matrix.
53*> \endverbatim
54*>
55*> \param[in] NB2
56*> \verbatim
57*> NB2 is INTEGER
58*> Number of columns in column block in an output test matrix.
59*> \endverbatim
60*>
61*> \param[out] RESULT
62*> \verbatim
63*> RESULT is REAL array, dimension (6)
64*> Results of each of the six tests below.
65*>
66*> A is a m-by-n test input matrix to be factored.
67*> so that A = Q_gr * ( R )
68*> ( 0 ),
69*>
70*> Q_qr is an implicit m-by-m unitary Q matrix, the result
71*> of factorization in blocked WY-representation,
72*> stored in CGEQRT output format.
73*>
74*> R is a n-by-n upper-triangular matrix,
75*>
76*> 0 is a (m-n)-by-n zero matrix,
77*>
78*> Q is an explicit m-by-m unitary matrix Q = Q_gr * I
79*>
80*> C is an m-by-n random matrix,
81*>
82*> D is an n-by-m random matrix.
83*>
84*> The six tests are:
85*>
86*> RESULT(1) = |R - (Q**H) * A| / ( eps * m * |A| )
87*> is equivalent to test for | A - Q * R | / (eps * m * |A|),
88*>
89*> RESULT(2) = |I - (Q**H) * Q| / ( eps * m ),
90*>
91*> RESULT(3) = | Q_qr * C - Q * C | / (eps * m * |C|),
92*>
93*> RESULT(4) = | (Q_gr**H) * C - (Q**H) * C | / (eps * m * |C|)
94*>
95*> RESULT(5) = | D * Q_qr - D * Q | / (eps * m * |D|)
96*>
97*> RESULT(6) = | D * (Q_qr**H) - D * (Q**H) | / (eps * m * |D|),
98*>
99*> where:
100*> Q_qr * C, (Q_gr**H) * C, D * Q_qr, D * (Q_qr**H) are
101*> computed using CGEMQRT,
102*>
103*> Q * C, (Q**H) * C, D * Q, D * (Q**H) are
104*> computed using CGEMM.
105*> \endverbatim
106*
107* Authors:
108* ========
109*
110*> \author Univ. of Tennessee
111*> \author Univ. of California Berkeley
112*> \author Univ. of Colorado Denver
113*> \author NAG Ltd.
114*
115*> \ingroup complex_lin
116*
117* =====================================================================
118 SUBROUTINE cunhr_col01( M, N, MB1, NB1, NB2, RESULT )
119 IMPLICIT NONE
120*
121* -- LAPACK test routine --
122* -- LAPACK is a software package provided by Univ. of Tennessee, --
123* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
124*
125* .. Scalar Arguments ..
126 INTEGER M, N, MB1, NB1, NB2
127* .. Return values ..
128 REAL RESULT(6)
129*
130* =====================================================================
131*
132* ..
133* .. Local allocatable arrays
134 COMPLEX , ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:),
135 $ WORK( : ), T1(:,:), T2(:,:), DIAG(:),
136 $ C(:,:), CF(:,:), D(:,:), DF(:,:)
137 REAL , ALLOCATABLE :: RWORK(:)
138*
139* .. Parameters ..
140 REAL ZERO
141 parameter( zero = 0.0e+0 )
142 COMPLEX CONE, CZERO
143 parameter( cone = ( 1.0e+0, 0.0e+0 ),
144 $ czero = ( 0.0e+0, 0.0e+0 ) )
145* ..
146* .. Local Scalars ..
147 LOGICAL TESTZEROS
148 INTEGER INFO, I, J, K, L, LWORK, NB1_UB, NB2_UB, NRB
149 REAL ANORM, EPS, RESID, CNORM, DNORM
150* ..
151* .. Local Arrays ..
152 INTEGER ISEED( 4 )
153 COMPLEX WORKQUERY( 1 )
154* ..
155* .. External Functions ..
156 REAL SLAMCH, CLANGE, CLANSY
157 EXTERNAL slamch, clange, clansy
158* ..
159* .. External Subroutines ..
160 EXTERNAL clacpy, clarnv, claset, clatsqr, cunhr_col,
162* ..
163* .. Intrinsic Functions ..
164 INTRINSIC ceiling, real, max, min
165* ..
166* .. Scalars in Common ..
167 CHARACTER(LEN=32) SRNAMT
168* ..
169* .. Common blocks ..
170 COMMON / srmnamc / srnamt
171* ..
172* .. Data statements ..
173 DATA iseed / 1988, 1989, 1990, 1991 /
174*
175* TEST MATRICES WITH HALF OF MATRIX BEING ZEROS
176*
177 testzeros = .false.
178*
179 eps = slamch( 'Epsilon' )
180 k = min( m, n )
181 l = max( m, n, 1)
182*
183* Dynamically allocate local arrays
184*
185 ALLOCATE ( a(m,n), af(m,n), q(l,l), r(m,l), rwork(l),
186 $ c(m,n), cf(m,n),
187 $ d(n,m), df(n,m) )
188*
189* Put random numbers into A and copy to AF
190*
191 DO j = 1, n
192 CALL clarnv( 2, iseed, m, a( 1, j ) )
193 END DO
194 IF( testzeros ) THEN
195 IF( m.GE.4 ) THEN
196 DO j = 1, n
197 CALL clarnv( 2, iseed, m/2, a( m/4, j ) )
198 END DO
199 END IF
200 END IF
201 CALL clacpy( 'Full', m, n, a, m, af, m )
202*
203* Number of row blocks in CLATSQR
204*
205 nrb = max( 1, ceiling( real( m - n ) / real( mb1 - n ) ) )
206*
207 ALLOCATE ( t1( nb1, n * nrb ) )
208 ALLOCATE ( t2( nb2, n ) )
209 ALLOCATE ( diag( n ) )
210*
211* Begin determine LWORK for the array WORK and allocate memory.
212*
213* CLATSQR requires NB1 to be bounded by N.
214*
215 nb1_ub = min( nb1, n)
216*
217* CGEMQRT requires NB2 to be bounded by N.
218*
219 nb2_ub = min( nb2, n)
220*
221 CALL clatsqr( m, n, mb1, nb1_ub, af, m, t1, nb1,
222 $ workquery, -1, info )
223 lwork = int( workquery( 1 ) )
224 CALL cungtsqr( m, n, mb1, nb1, af, m, t1, nb1, workquery, -1,
225 $ info )
226
227 lwork = max( lwork, int( workquery( 1 ) ) )
228*
229* In CGEMQRT, WORK is N*NB2_UB if SIDE = 'L',
230* or M*NB2_UB if SIDE = 'R'.
231*
232 lwork = max( lwork, nb2_ub * n, nb2_ub * m )
233*
234 ALLOCATE ( work( lwork ) )
235*
236* End allocate memory for WORK.
237*
238*
239* Begin Householder reconstruction routines
240*
241* Factor the matrix A in the array AF.
242*
243 srnamt = 'CLATSQR'
244 CALL clatsqr( m, n, mb1, nb1_ub, af, m, t1, nb1, work, lwork,
245 $ info )
246*
247* Copy the factor R into the array R.
248*
249 srnamt = 'CLACPY'
250 CALL clacpy( 'U', n, n, af, m, r, m )
251*
252* Reconstruct the orthogonal matrix Q.
253*
254 srnamt = 'CUNGTSQR'
255 CALL cungtsqr( m, n, mb1, nb1, af, m, t1, nb1, work, lwork,
256 $ info )
257*
258* Perform the Householder reconstruction, the result is stored
259* the arrays AF and T2.
260*
261 srnamt = 'CUNHR_COL'
262 CALL cunhr_col( m, n, nb2, af, m, t2, nb2, diag, info )
263*
264* Compute the factor R_hr corresponding to the Householder
265* reconstructed Q_hr and place it in the upper triangle of AF to
266* match the Q storage format in CGEQRT. R_hr = R_tsqr * S,
267* this means changing the sign of I-th row of the matrix R_tsqr
268* according to sign of of I-th diagonal element DIAG(I) of the
269* matrix S.
270*
271 srnamt = 'CLACPY'
272 CALL clacpy( 'U', n, n, r, m, af, m )
273*
274 DO i = 1, n
275 IF( diag( i ).EQ.-cone ) THEN
276 CALL cscal( n+1-i, -cone, af( i, i ), m )
277 END IF
278 END DO
279*
280* End Householder reconstruction routines.
281*
282*
283* Generate the m-by-m matrix Q
284*
285 CALL claset( 'Full', m, m, czero, cone, q, m )
286*
287 srnamt = 'CGEMQRT'
288 CALL cgemqrt( 'L', 'N', m, m, k, nb2_ub, af, m, t2, nb2, q, m,
289 $ work, info )
290*
291* Copy R
292*
293 CALL claset( 'Full', m, n, czero, czero, r, m )
294*
295 CALL clacpy( 'Upper', m, n, af, m, r, m )
296*
297* TEST 1
298* Compute |R - (Q**H)*A| / ( eps * m * |A| ) and store in RESULT(1)
299*
300 CALL cgemm( 'C', 'N', m, n, m, -cone, q, m, a, m, cone, r, m )
301*
302 anorm = clange( '1', m, n, a, m, rwork )
303 resid = clange( '1', m, n, r, m, rwork )
304 IF( anorm.GT.zero ) THEN
305 result( 1 ) = resid / ( eps * max( 1, m ) * anorm )
306 ELSE
307 result( 1 ) = zero
308 END IF
309*
310* TEST 2
311* Compute |I - (Q**H)*Q| / ( eps * m ) and store in RESULT(2)
312*
313 CALL claset( 'Full', m, m, czero, cone, r, m )
314 CALL cherk( 'U', 'C', m, m, real(-cone), q, m, real(cone), r, m )
315 resid = clansy( '1', 'Upper', m, r, m, rwork )
316 result( 2 ) = resid / ( eps * max( 1, m ) )
317*
318* Generate random m-by-n matrix C
319*
320 DO j = 1, n
321 CALL clarnv( 2, iseed, m, c( 1, j ) )
322 END DO
323 cnorm = clange( '1', m, n, c, m, rwork )
324 CALL clacpy( 'Full', m, n, c, m, cf, m )
325*
326* Apply Q to C as Q*C = CF
327*
328 srnamt = 'CGEMQRT'
329 CALL cgemqrt( 'L', 'N', m, n, k, nb2_ub, af, m, t2, nb2, cf, m,
330 $ work, info )
331*
332* TEST 3
333* Compute |CF - Q*C| / ( eps * m * |C| )
334*
335 CALL cgemm( 'N', 'N', m, n, m, -cone, q, m, c, m, cone, cf, m )
336 resid = clange( '1', m, n, cf, m, rwork )
337 IF( cnorm.GT.zero ) THEN
338 result( 3 ) = resid / ( eps * max( 1, m ) * cnorm )
339 ELSE
340 result( 3 ) = zero
341 END IF
342*
343* Copy C into CF again
344*
345 CALL clacpy( 'Full', m, n, c, m, cf, m )
346*
347* Apply Q to C as (Q**H)*C = CF
348*
349 srnamt = 'CGEMQRT'
350 CALL cgemqrt( 'L', 'C', m, n, k, nb2_ub, af, m, t2, nb2, cf, m,
351 $ work, info )
352*
353* TEST 4
354* Compute |CF - (Q**H)*C| / ( eps * m * |C|)
355*
356 CALL cgemm( 'C', 'N', m, n, m, -cone, q, m, c, m, cone, cf, m )
357 resid = clange( '1', m, n, cf, m, rwork )
358 IF( cnorm.GT.zero ) THEN
359 result( 4 ) = resid / ( eps * max( 1, m ) * cnorm )
360 ELSE
361 result( 4 ) = zero
362 END IF
363*
364* Generate random n-by-m matrix D and a copy DF
365*
366 DO j = 1, m
367 CALL clarnv( 2, iseed, n, d( 1, j ) )
368 END DO
369 dnorm = clange( '1', n, m, d, n, rwork )
370 CALL clacpy( 'Full', n, m, d, n, df, n )
371*
372* Apply Q to D as D*Q = DF
373*
374 srnamt = 'CGEMQRT'
375 CALL cgemqrt( 'R', 'N', n, m, k, nb2_ub, af, m, t2, nb2, df, n,
376 $ work, info )
377*
378* TEST 5
379* Compute |DF - D*Q| / ( eps * m * |D| )
380*
381 CALL cgemm( 'N', 'N', n, m, m, -cone, d, n, q, m, cone, df, n )
382 resid = clange( '1', n, m, df, n, rwork )
383 IF( dnorm.GT.zero ) THEN
384 result( 5 ) = resid / ( eps * max( 1, m ) * dnorm )
385 ELSE
386 result( 5 ) = zero
387 END IF
388*
389* Copy D into DF again
390*
391 CALL clacpy( 'Full', n, m, d, n, df, n )
392*
393* Apply Q to D as D*QT = DF
394*
395 srnamt = 'CGEMQRT'
396 CALL cgemqrt( 'R', 'C', n, m, k, nb2_ub, af, m, t2, nb2, df, n,
397 $ work, info )
398*
399* TEST 6
400* Compute |DF - D*(Q**H)| / ( eps * m * |D| )
401*
402 CALL cgemm( 'N', 'C', n, m, m, -cone, d, n, q, m, cone, df, n )
403 resid = clange( '1', n, m, df, n, rwork )
404 IF( dnorm.GT.zero ) THEN
405 result( 6 ) = resid / ( eps * max( 1, m ) * dnorm )
406 ELSE
407 result( 6 ) = zero
408 END IF
409*
410* Deallocate all arrays
411*
412 DEALLOCATE ( a, af, q, r, rwork, work, t1, t2, diag,
413 $ c, d, cf, df )
414*
415 RETURN
416*
417* End of CUNHR_COL01
418*
419 END
subroutine cunhr_col01(m, n, mb1, nb1, nb2, result)
CUNHR_COL01
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
subroutine cgemqrt(side, trans, m, n, k, nb, v, ldv, t, ldt, c, ldc, work, info)
CGEMQRT
Definition cgemqrt.f:168
subroutine cherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
CHERK
Definition cherk.f:173
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
subroutine clarnv(idist, iseed, n, x)
CLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition clarnv.f:99
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:106
subroutine clatsqr(m, n, mb, nb, a, lda, t, ldt, work, lwork, info)
CLATSQR
Definition clatsqr.f:169
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
subroutine cungtsqr(m, n, mb, nb, a, lda, t, ldt, work, lwork, info)
CUNGTSQR
Definition cungtsqr.f:176
subroutine cunhr_col(m, n, nb, a, lda, t, ldt, d, info)
CUNHR_COL
Definition cunhr_col.f:259