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,
161     $                   CUNGTSQR, CSCAL, CGEMM, CGEMQRT, CHERK
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, -CONE, Q, M, 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
420