1*> \brief \b DQRT05 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 DQRT05(M,N,L,NB,RESULT) 12* 13* .. Scalar Arguments .. 14* INTEGER LWORK, M, N, L, NB, LDT 15* .. Return values .. 16* DOUBLE PRECISION RESULT(6) 17* 18* 19*> \par Purpose: 20* ============= 21*> 22*> \verbatim 23*> 24*> DQRT05 tests DTPQRT and DTPMQRT. 25*> \endverbatim 26* 27* Arguments: 28* ========== 29* 30*> \param[in] M 31*> \verbatim 32*> M is INTEGER 33*> Number of rows in lower part of the test matrix. 34*> \endverbatim 35*> 36*> \param[in] N 37*> \verbatim 38*> N is INTEGER 39*> Number of columns in test matrix. 40*> \endverbatim 41*> 42*> \param[in] L 43*> \verbatim 44*> L is INTEGER 45*> The number of rows of the upper trapezoidal part the 46*> lower test matrix. 0 <= L <= M. 47*> \endverbatim 48*> 49*> \param[in] NB 50*> \verbatim 51*> NB is INTEGER 52*> Block size of test matrix. NB <= N. 53*> \endverbatim 54*> 55*> \param[out] RESULT 56*> \verbatim 57*> RESULT is DOUBLE PRECISION array, dimension (6) 58*> Results of each of the six tests below. 59*> 60*> RESULT(1) = | A - Q R | 61*> RESULT(2) = | I - Q^H Q | 62*> RESULT(3) = | Q C - Q C | 63*> RESULT(4) = | Q^H C - Q^H C | 64*> RESULT(5) = | C Q - C Q | 65*> RESULT(6) = | C Q^H - C Q^H | 66*> \endverbatim 67* 68* Authors: 69* ======== 70* 71*> \author Univ. of Tennessee 72*> \author Univ. of California Berkeley 73*> \author Univ. of Colorado Denver 74*> \author NAG Ltd. 75* 76*> \date April 2012 77* 78*> \ingroup double_lin 79* 80* ===================================================================== 81 SUBROUTINE DQRT05(M,N,L,NB,RESULT) 82 IMPLICIT NONE 83* 84* -- LAPACK test routine (version 3.4.1) -- 85* -- LAPACK is a software package provided by Univ. of Tennessee, -- 86* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 87* April 2012 88* 89* .. Scalar Arguments .. 90 INTEGER LWORK, M, N, L, NB, LDT 91* .. Return values .. 92 DOUBLE PRECISION RESULT(6) 93* 94* ===================================================================== 95* 96* .. 97* .. Local allocatable arrays 98 DOUBLE PRECISION, ALLOCATABLE :: AF(:,:), Q(:,:), 99 $ R(:,:), RWORK(:), WORK( : ), T(:,:), 100 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:) 101* 102* .. Parameters .. 103 DOUBLE PRECISION ONE, ZERO 104 PARAMETER( ZERO = 0.0, ONE = 1.0 ) 105* .. 106* .. Local Scalars .. 107 INTEGER INFO, J, K, M2, NP1 108 DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM 109* .. 110* .. Local Arrays .. 111 INTEGER ISEED( 4 ) 112* .. 113* .. External Functions .. 114 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY 115 LOGICAL LSAME 116 EXTERNAL DLAMCH, DLANGE, DLANSY, LSAME 117* .. 118* .. Data statements .. 119 DATA ISEED / 1988, 1989, 1990, 1991 / 120* 121 EPS = DLAMCH( 'Epsilon' ) 122 K = N 123 M2 = M+N 124 IF( M.GT.0 ) THEN 125 NP1 = N+1 126 ELSE 127 NP1 = 1 128 END IF 129 LWORK = M2*M2*NB 130* 131* Dynamically allocate all arrays 132* 133 ALLOCATE(A(M2,N),AF(M2,N),Q(M2,M2),R(M2,M2),RWORK(M2), 134 $ WORK(LWORK),T(NB,N),C(M2,N),CF(M2,N), 135 $ D(N,M2),DF(N,M2) ) 136* 137* Put random stuff into A 138* 139 LDT=NB 140 CALL DLASET( 'Full', M2, N, ZERO, ZERO, A, M2 ) 141 CALL DLASET( 'Full', NB, N, ZERO, ZERO, T, NB ) 142 DO J=1,N 143 CALL DLARNV( 2, ISEED, J, A( 1, J ) ) 144 END DO 145 IF( M.GT.0 ) THEN 146 DO J=1,N 147 CALL DLARNV( 2, ISEED, M-L, A( MIN(N+M,N+1), J ) ) 148 END DO 149 END IF 150 IF( L.GT.0 ) THEN 151 DO J=1,N 152 CALL DLARNV( 2, ISEED, MIN(J,L), A( MIN(N+M,N+M-L+1), J ) ) 153 END DO 154 END IF 155* 156* Copy the matrix A to the array AF. 157* 158 CALL DLACPY( 'Full', M2, N, A, M2, AF, M2 ) 159* 160* Factor the matrix A in the array AF. 161* 162 CALL DTPQRT( M,N,L,NB,AF,M2,AF(NP1,1),M2,T,LDT,WORK,INFO) 163* 164* Generate the (M+N)-by-(M+N) matrix Q by applying H to I 165* 166 CALL DLASET( 'Full', M2, M2, ZERO, ONE, Q, M2 ) 167 CALL DGEMQRT( 'R', 'N', M2, M2, K, NB, AF, M2, T, LDT, Q, M2, 168 $ WORK, INFO ) 169* 170* Copy R 171* 172 CALL DLASET( 'Full', M2, N, ZERO, ZERO, R, M2 ) 173 CALL DLACPY( 'Upper', M2, N, AF, M2, R, M2 ) 174* 175* Compute |R - Q'*A| / |A| and store in RESULT(1) 176* 177 CALL DGEMM( 'T', 'N', M2, N, M2, -ONE, Q, M2, A, M2, ONE, R, M2 ) 178 ANORM = DLANGE( '1', M2, N, A, M2, RWORK ) 179 RESID = DLANGE( '1', M2, N, R, M2, RWORK ) 180 IF( ANORM.GT.ZERO ) THEN 181 RESULT( 1 ) = RESID / (EPS*ANORM*MAX(1,M2)) 182 ELSE 183 RESULT( 1 ) = ZERO 184 END IF 185* 186* Compute |I - Q'*Q| and store in RESULT(2) 187* 188 CALL DLASET( 'Full', M2, M2, ZERO, ONE, R, M2 ) 189 CALL DSYRK( 'U', 'C', M2, M2, -ONE, Q, M2, ONE, R, M2 ) 190 RESID = DLANSY( '1', 'Upper', M2, R, M2, RWORK ) 191 RESULT( 2 ) = RESID / (EPS*MAX(1,M2)) 192* 193* Generate random m-by-n matrix C and a copy CF 194* 195 DO J=1,N 196 CALL DLARNV( 2, ISEED, M2, C( 1, J ) ) 197 END DO 198 CNORM = DLANGE( '1', M2, N, C, M2, RWORK) 199 CALL DLACPY( 'Full', M2, N, C, M2, CF, M2 ) 200* 201* Apply Q to C as Q*C 202* 203 CALL DTPMQRT( 'L','N', M,N,K,L,NB,AF(NP1,1),M2,T,LDT,CF,M2, 204 $ CF(NP1,1),M2,WORK,INFO) 205* 206* Compute |Q*C - Q*C| / |C| 207* 208 CALL DGEMM( 'N', 'N', M2, N, M2, -ONE, Q, M2, C, M2, ONE, CF, M2 ) 209 RESID = DLANGE( '1', M2, N, CF, M2, RWORK ) 210 IF( CNORM.GT.ZERO ) THEN 211 RESULT( 3 ) = RESID / (EPS*MAX(1,M2)*CNORM) 212 ELSE 213 RESULT( 3 ) = ZERO 214 END IF 215* 216* Copy C into CF again 217* 218 CALL DLACPY( 'Full', M2, N, C, M2, CF, M2 ) 219* 220* Apply Q to C as QT*C 221* 222 CALL DTPMQRT( 'L','T',M,N,K,L,NB,AF(NP1,1),M2,T,LDT,CF,M2, 223 $ CF(NP1,1),M2,WORK,INFO) 224* 225* Compute |QT*C - QT*C| / |C| 226* 227 CALL DGEMM('T','N',M2,N,M2,-ONE,Q,M2,C,M2,ONE,CF,M2) 228 RESID = DLANGE( '1', M2, N, CF, M2, RWORK ) 229 IF( CNORM.GT.ZERO ) THEN 230 RESULT( 4 ) = RESID / (EPS*MAX(1,M2)*CNORM) 231 ELSE 232 RESULT( 4 ) = ZERO 233 END IF 234* 235* Generate random n-by-m matrix D and a copy DF 236* 237 DO J=1,M2 238 CALL DLARNV( 2, ISEED, N, D( 1, J ) ) 239 END DO 240 DNORM = DLANGE( '1', N, M2, D, N, RWORK) 241 CALL DLACPY( 'Full', N, M2, D, N, DF, N ) 242* 243* Apply Q to D as D*Q 244* 245 CALL DTPMQRT('R','N',N,M,N,L,NB,AF(NP1,1),M2,T,LDT,DF,N, 246 $ DF(1,NP1),N,WORK,INFO) 247* 248* Compute |D*Q - D*Q| / |D| 249* 250 CALL DGEMM('N','N',N,M2,M2,-ONE,D,N,Q,M2,ONE,DF,N) 251 RESID = DLANGE('1',N, M2,DF,N,RWORK ) 252 IF( CNORM.GT.ZERO ) THEN 253 RESULT( 5 ) = RESID / (EPS*MAX(1,M2)*DNORM) 254 ELSE 255 RESULT( 5 ) = ZERO 256 END IF 257* 258* Copy D into DF again 259* 260 CALL DLACPY('Full',N,M2,D,N,DF,N ) 261* 262* Apply Q to D as D*QT 263* 264 CALL DTPMQRT('R','T',N,M,N,L,NB,AF(NP1,1),M2,T,LDT,DF,N, 265 $ DF(1,NP1),N,WORK,INFO) 266 267* 268* Compute |D*QT - D*QT| / |D| 269* 270 CALL DGEMM( 'N', 'T', N, M2, M2, -ONE, D, N, Q, M2, ONE, DF, N ) 271 RESID = DLANGE( '1', N, M2, DF, N, RWORK ) 272 IF( CNORM.GT.ZERO ) THEN 273 RESULT( 6 ) = RESID / (EPS*MAX(1,M2)*DNORM) 274 ELSE 275 RESULT( 6 ) = ZERO 276 END IF 277* 278* Deallocate all arrays 279* 280 DEALLOCATE ( A, AF, Q, R, RWORK, WORK, T, C, D, CF, DF) 281 RETURN 282 END 283 284