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