1*> \brief \b DORHR_COL02 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 DORHR_COL02( 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*> DORHR_COL02 tests DORGTSQR_ROW and DORHR_COL inside DGETSQRHRT 25*> (which calls DLATSQR, DORGTSQR_ROW and DORHR_COL) using DGEMQRT. 26*> Therefore, DLATSQR (part of DGEQR), DGEMQRT (part of DGEMQR) 27*> have to be tested before this test. 28*> 29*> \endverbatim 30* 31* Arguments: 32* ========== 33* 34*> \param[in] M 35*> \verbatim 36*> M is INTEGER 37*> Number of rows in test matrix. 38*> \endverbatim 39*> \param[in] N 40*> \verbatim 41*> N is INTEGER 42*> Number of columns in test matrix. 43*> \endverbatim 44*> \param[in] MB1 45*> \verbatim 46*> MB1 is INTEGER 47*> Number of row in row block in an input test matrix. 48*> \endverbatim 49*> 50*> \param[in] NB1 51*> \verbatim 52*> NB1 is INTEGER 53*> Number of columns in column block an input test matrix. 54*> \endverbatim 55*> 56*> \param[in] NB2 57*> \verbatim 58*> NB2 is INTEGER 59*> Number of columns in column block in an output test matrix. 60*> \endverbatim 61*> 62*> \param[out] RESULT 63*> \verbatim 64*> RESULT is DOUBLE PRECISION array, dimension (6) 65*> Results of each of the six tests below. 66*> 67*> A is a m-by-n test input matrix to be factored. 68*> so that A = Q_gr * ( R ) 69*> ( 0 ), 70*> 71*> Q_qr is an implicit m-by-m orthogonal Q matrix, the result 72*> of factorization in blocked WY-representation, 73*> stored in ZGEQRT output format. 74*> 75*> R is a n-by-n upper-triangular matrix, 76*> 77*> 0 is a (m-n)-by-n zero matrix, 78*> 79*> Q is an explicit m-by-m orthogonal matrix Q = Q_gr * I 80*> 81*> C is an m-by-n random matrix, 82*> 83*> D is an n-by-m random matrix. 84*> 85*> The six tests are: 86*> 87*> RESULT(1) = |R - (Q**H) * A| / ( eps * m * |A| ) 88*> is equivalent to test for | A - Q * R | / (eps * m * |A|), 89*> 90*> RESULT(2) = |I - (Q**H) * Q| / ( eps * m ), 91*> 92*> RESULT(3) = | Q_qr * C - Q * C | / (eps * m * |C|), 93*> 94*> RESULT(4) = | (Q_gr**H) * C - (Q**H) * C | / (eps * m * |C|) 95*> 96*> RESULT(5) = | D * Q_qr - D * Q | / (eps * m * |D|) 97*> 98*> RESULT(6) = | D * (Q_qr**H) - D * (Q**H) | / (eps * m * |D|), 99*> 100*> where: 101*> Q_qr * C, (Q_gr**H) * C, D * Q_qr, D * (Q_qr**H) are 102*> computed using DGEMQRT, 103*> 104*> Q * C, (Q**H) * C, D * Q, D * (Q**H) are 105*> computed using DGEMM. 106*> \endverbatim 107* 108* Authors: 109* ======== 110* 111*> \author Univ. of Tennessee 112*> \author Univ. of California Berkeley 113*> \author Univ. of Colorado Denver 114*> \author NAG Ltd. 115* 116*> \ingroup double_lin 117* 118* ===================================================================== 119 SUBROUTINE DORHR_COL02( M, N, MB1, NB1, NB2, RESULT ) 120 IMPLICIT NONE 121* 122* -- LAPACK test routine -- 123* -- LAPACK is a software package provided by Univ. of Tennessee, -- 124* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 125* 126* .. Scalar Arguments .. 127 INTEGER M, N, MB1, NB1, NB2 128* .. Return values .. 129 DOUBLE PRECISION RESULT(6) 130* 131* ===================================================================== 132* 133* .. 134* .. Local allocatable arrays 135 DOUBLE PRECISION, ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:), 136 $ RWORK(:), WORK( : ), T1(:,:), T2(:,:), DIAG(:), 137 $ C(:,:), CF(:,:), D(:,:), DF(:,:) 138* 139* .. Parameters .. 140 DOUBLE PRECISION ONE, ZERO 141 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 142* .. 143* .. Local Scalars .. 144 LOGICAL TESTZEROS 145 INTEGER INFO, J, K, L, LWORK, NB2_UB, NRB 146 DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM 147* .. 148* .. Local Arrays .. 149 INTEGER ISEED( 4 ) 150 DOUBLE PRECISION WORKQUERY( 1 ) 151* .. 152* .. External Functions .. 153 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY 154 EXTERNAL DLAMCH, DLANGE, DLANSY 155* .. 156* .. External Subroutines .. 157 EXTERNAL DLACPY, DLARNV, DLASET, DGETSQRHRT, 158 $ DSCAL, DGEMM, DGEMQRT, DSYRK 159* .. 160* .. Intrinsic Functions .. 161 INTRINSIC CEILING, DBLE, MAX, MIN 162* .. 163* .. Scalars in Common .. 164 CHARACTER(LEN=32) SRNAMT 165* .. 166* .. Common blocks .. 167 COMMON / SRMNAMC / SRNAMT 168* .. 169* .. Data statements .. 170 DATA ISEED / 1988, 1989, 1990, 1991 / 171* 172* TEST MATRICES WITH HALF OF MATRIX BEING ZEROS 173* 174 TESTZEROS = .FALSE. 175* 176 EPS = DLAMCH( 'Epsilon' ) 177 K = MIN( M, N ) 178 L = MAX( M, N, 1) 179* 180* Dynamically allocate local arrays 181* 182 ALLOCATE ( A(M,N), AF(M,N), Q(L,L), R(M,L), RWORK(L), 183 $ C(M,N), CF(M,N), 184 $ D(N,M), DF(N,M) ) 185* 186* Put random numbers into A and copy to AF 187* 188 DO J = 1, N 189 CALL DLARNV( 2, ISEED, M, A( 1, J ) ) 190 END DO 191 IF( TESTZEROS ) THEN 192 IF( M.GE.4 ) THEN 193 DO J = 1, N 194 CALL DLARNV( 2, ISEED, M/2, A( M/4, J ) ) 195 END DO 196 END IF 197 END IF 198 CALL DLACPY( 'Full', M, N, A, M, AF, M ) 199* 200* Number of row blocks in DLATSQR 201* 202 NRB = MAX( 1, CEILING( DBLE( M - N ) / DBLE( MB1 - N ) ) ) 203* 204 ALLOCATE ( T1( NB1, N * NRB ) ) 205 ALLOCATE ( T2( NB2, N ) ) 206 ALLOCATE ( DIAG( N ) ) 207* 208* Begin determine LWORK for the array WORK and allocate memory. 209* 210* DGEMQRT requires NB2 to be bounded by N. 211* 212 NB2_UB = MIN( NB2, N) 213* 214* 215 CALL DGETSQRHRT( M, N, MB1, NB1, NB2, AF, M, T2, NB2, 216 $ WORKQUERY, -1, INFO ) 217* 218 LWORK = INT( WORKQUERY( 1 ) ) 219* 220* In DGEMQRT, WORK is N*NB2_UB if SIDE = 'L', 221* or M*NB2_UB if SIDE = 'R'. 222* 223 LWORK = MAX( LWORK, NB2_UB * N, NB2_UB * M ) 224* 225 ALLOCATE ( WORK( LWORK ) ) 226* 227* End allocate memory for WORK. 228* 229* 230* Begin Householder reconstruction routines 231* 232* Factor the matrix A in the array AF. 233* 234 SRNAMT = 'DGETSQRHRT' 235 CALL DGETSQRHRT( M, N, MB1, NB1, NB2, AF, M, T2, NB2, 236 $ WORK, LWORK, INFO ) 237* 238* End Householder reconstruction routines. 239* 240* 241* Generate the m-by-m matrix Q 242* 243 CALL DLASET( 'Full', M, M, ZERO, ONE, Q, M ) 244* 245 SRNAMT = 'DGEMQRT' 246 CALL DGEMQRT( 'L', 'N', M, M, K, NB2_UB, AF, M, T2, NB2, Q, M, 247 $ WORK, INFO ) 248* 249* Copy R 250* 251 CALL DLASET( 'Full', M, N, ZERO, ZERO, R, M ) 252* 253 CALL DLACPY( 'Upper', M, N, AF, M, R, M ) 254* 255* TEST 1 256* Compute |R - (Q**T)*A| / ( eps * m * |A| ) and store in RESULT(1) 257* 258 CALL DGEMM( 'T', 'N', M, N, M, -ONE, Q, M, A, M, ONE, R, M ) 259* 260 ANORM = DLANGE( '1', M, N, A, M, RWORK ) 261 RESID = DLANGE( '1', M, N, R, M, RWORK ) 262 IF( ANORM.GT.ZERO ) THEN 263 RESULT( 1 ) = RESID / ( EPS * MAX( 1, M ) * ANORM ) 264 ELSE 265 RESULT( 1 ) = ZERO 266 END IF 267* 268* TEST 2 269* Compute |I - (Q**T)*Q| / ( eps * m ) and store in RESULT(2) 270* 271 CALL DLASET( 'Full', M, M, ZERO, ONE, R, M ) 272 CALL DSYRK( 'U', 'T', M, M, -ONE, Q, M, ONE, R, M ) 273 RESID = DLANSY( '1', 'Upper', M, R, M, RWORK ) 274 RESULT( 2 ) = RESID / ( EPS * MAX( 1, M ) ) 275* 276* Generate random m-by-n matrix C 277* 278 DO J = 1, N 279 CALL DLARNV( 2, ISEED, M, C( 1, J ) ) 280 END DO 281 CNORM = DLANGE( '1', M, N, C, M, RWORK ) 282 CALL DLACPY( 'Full', M, N, C, M, CF, M ) 283* 284* Apply Q to C as Q*C = CF 285* 286 SRNAMT = 'DGEMQRT' 287 CALL DGEMQRT( 'L', 'N', M, N, K, NB2_UB, AF, M, T2, NB2, CF, M, 288 $ WORK, INFO ) 289* 290* TEST 3 291* Compute |CF - Q*C| / ( eps * m * |C| ) 292* 293 CALL DGEMM( 'N', 'N', M, N, M, -ONE, Q, M, C, M, ONE, CF, M ) 294 RESID = DLANGE( '1', M, N, CF, M, RWORK ) 295 IF( CNORM.GT.ZERO ) THEN 296 RESULT( 3 ) = RESID / ( EPS * MAX( 1, M ) * CNORM ) 297 ELSE 298 RESULT( 3 ) = ZERO 299 END IF 300* 301* Copy C into CF again 302* 303 CALL DLACPY( 'Full', M, N, C, M, CF, M ) 304* 305* Apply Q to C as (Q**T)*C = CF 306* 307 SRNAMT = 'DGEMQRT' 308 CALL DGEMQRT( 'L', 'T', M, N, K, NB2_UB, AF, M, T2, NB2, CF, M, 309 $ WORK, INFO ) 310* 311* TEST 4 312* Compute |CF - (Q**T)*C| / ( eps * m * |C|) 313* 314 CALL DGEMM( 'T', 'N', M, N, M, -ONE, Q, M, C, M, ONE, CF, M ) 315 RESID = DLANGE( '1', M, N, CF, M, RWORK ) 316 IF( CNORM.GT.ZERO ) THEN 317 RESULT( 4 ) = RESID / ( EPS * MAX( 1, M ) * CNORM ) 318 ELSE 319 RESULT( 4 ) = ZERO 320 END IF 321* 322* Generate random n-by-m matrix D and a copy DF 323* 324 DO J = 1, M 325 CALL DLARNV( 2, ISEED, N, D( 1, J ) ) 326 END DO 327 DNORM = DLANGE( '1', N, M, D, N, RWORK ) 328 CALL DLACPY( 'Full', N, M, D, N, DF, N ) 329* 330* Apply Q to D as D*Q = DF 331* 332 SRNAMT = 'DGEMQRT' 333 CALL DGEMQRT( 'R', 'N', N, M, K, NB2_UB, AF, M, T2, NB2, DF, N, 334 $ WORK, INFO ) 335* 336* TEST 5 337* Compute |DF - D*Q| / ( eps * m * |D| ) 338* 339 CALL DGEMM( 'N', 'N', N, M, M, -ONE, D, N, Q, M, ONE, DF, N ) 340 RESID = DLANGE( '1', N, M, DF, N, RWORK ) 341 IF( DNORM.GT.ZERO ) THEN 342 RESULT( 5 ) = RESID / ( EPS * MAX( 1, M ) * DNORM ) 343 ELSE 344 RESULT( 5 ) = ZERO 345 END IF 346* 347* Copy D into DF again 348* 349 CALL DLACPY( 'Full', N, M, D, N, DF, N ) 350* 351* Apply Q to D as D*QT = DF 352* 353 SRNAMT = 'DGEMQRT' 354 CALL DGEMQRT( 'R', 'T', N, M, K, NB2_UB, AF, M, T2, NB2, DF, N, 355 $ WORK, INFO ) 356* 357* TEST 6 358* Compute |DF - D*(Q**T)| / ( eps * m * |D| ) 359* 360 CALL DGEMM( 'N', 'T', N, M, M, -ONE, D, N, Q, M, ONE, DF, N ) 361 RESID = DLANGE( '1', N, M, DF, N, RWORK ) 362 IF( DNORM.GT.ZERO ) THEN 363 RESULT( 6 ) = RESID / ( EPS * MAX( 1, M ) * DNORM ) 364 ELSE 365 RESULT( 6 ) = ZERO 366 END IF 367* 368* Deallocate all arrays 369* 370 DEALLOCATE ( A, AF, Q, R, RWORK, WORK, T1, T2, DIAG, 371 $ C, D, CF, DF ) 372* 373 RETURN 374* 375* End of DORHR_COL02 376* 377 END 378