TESTING/LIN/zqrt05.f

   1 *> \brief \b ZQRT05
   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 ZQRT05(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 *> ZQRT05 tests ZTPQRT and ZTPMQRT.
  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 complex16_lin
  79 *
  80 *  =====================================================================
  81       SUBROUTINE ZQRT05(M,N,L,NB,RESULT)
  82       IMPLICIT NONE
  83 *
  84 *  -- LAPACK test routine (version 3.6.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       COMPLEX*16, ALLOCATABLE :: AF(:,:), Q(:,:),
  99      $  R(:,:), WORK( : ), T(:,:),
 100      $  CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:)
 101       DOUBLE PRECISION, ALLOCATABLE :: RWORK(:)
 102 *
 103 *     .. Parameters ..
 104       DOUBLE PRECISION ZERO
 105       COMPLEX*16 ONE, CZERO
 106       PARAMETER( ZERO = 0.0, ONE = (1.0,0.0), CZERO=(0.0,0.0) )
 107 *     ..
 108 *     .. Local Scalars ..
 109       INTEGER INFO, J, K, M2, NP1
 110       DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
 111 *     ..
 112 *     .. Local Arrays ..
 113       INTEGER            ISEED( 4 )
 114 *     ..
 115 *     .. External Functions ..
 116       DOUBLE PRECISION DLAMCH
 117       DOUBLE PRECISION ZLANGE, ZLANSY
 118       LOGICAL  LSAME
 119       EXTERNAL DLAMCH, ZLANGE, ZLANSY, LSAME
 120 *     ..
 121 *     .. Data statements ..
 122       DATA ISEED / 1988, 1989, 1990, 1991 /
 123 *
 124       EPS = DLAMCH( 'Epsilon' )
 125       K = N
 126       M2 = M+N
 127       IF( M.GT.0 ) THEN
 128          NP1 = N+1
 129       ELSE
 130          NP1 = 1
 131       END IF
 132       LWORK = M2*M2*NB
 133 *
 134 *     Dynamically allocate all arrays
 135 *
 136       ALLOCATE(A(M2,N),AF(M2,N),Q(M2,M2),R(M2,M2),RWORK(M2),
 137      $           WORK(LWORK),T(NB,N),C(M2,N),CF(M2,N),
 138      $           D(N,M2),DF(N,M2) )
 139 *
 140 *     Put random stuff into A
 141 *
 142       LDT=NB
 143       CALL ZLASET( 'Full', M2, N, CZERO, CZERO, A, M2 )
 144       CALL ZLASET( 'Full', NB, N, CZERO, CZERO, T, NB )
 145       DO J=1,N
 146          CALL ZLARNV( 2, ISEED, J, A( 1, J ) )
 147       END DO
 148       IF( M.GT.0 ) THEN
 149          DO J=1,N
 150             CALL ZLARNV( 2, ISEED, M-L, A( MIN(N+M,N+1), J ) )
 151          END DO
 152       END IF
 153       IF( L.GT.0 ) THEN
 154          DO J=1,N
 155             CALL ZLARNV( 2, ISEED, MIN(J,L), A( MIN(N+M,N+M-L+1), J ) )
 156          END DO
 157       END IF
 158 *
 159 *     Copy the matrix A to the array AF.
 160 *
 161       CALL ZLACPY( 'Full', M2, N, A, M2, AF, M2 )
 162 *
 163 *     Factor the matrix A in the array AF.
 164 *
 165       CALL ZTPQRT( M,N,L,NB,AF,M2,AF(NP1,1),M2,T,LDT,WORK,INFO)
 166 *
 167 *     Generate the (M+N)-by-(M+N) matrix Q by applying H to I
 168 *
 169       CALL ZLASET( 'Full', M2, M2, CZERO, ONE, Q, M2 )
 170       CALL ZGEMQRT( 'R', 'N', M2, M2, K, NB, AF, M2, T, LDT, Q, M2,
 171      $              WORK, INFO )
 172 *
 173 *     Copy R
 174 *
 175       CALL ZLASET( 'Full', M2, N, CZERO, CZERO, R, M2 )
 176       CALL ZLACPY( 'Upper', M2, N, AF, M2, R, M2 )
 177 *
 178 *     Compute |R - Q'*A| / |A| and store in RESULT(1)
 179 *
 180       CALL ZGEMM( 'C', 'N', M2, N, M2, -ONE, Q, M2, A, M2, ONE, R, M2 )
 181       ANORM = ZLANGE( '1', M2, N, A, M2, RWORK )
 182       RESID = ZLANGE( '1', M2, N, R, M2, RWORK )
 183       IF( ANORM.GT.ZERO ) THEN
 184          RESULT( 1 ) = RESID / (EPS*ANORM*MAX(1,M2))
 185       ELSE
 186          RESULT( 1 ) = ZERO
 187       END IF
 188 *
 189 *     Compute |I - Q'*Q| and store in RESULT(2)
 190 *
 191       CALL ZLASET( 'Full', M2, M2, CZERO, ONE, R, M2 )
 192       CALL ZHERK( 'U', 'C', M2, M2, DREAL(-ONE), Q, M2, DREAL(ONE),
 193      $            R, M2 )
 194       RESID = ZLANSY( '1', 'Upper', M2, R, M2, RWORK )
 195       RESULT( 2 ) = RESID / (EPS*MAX(1,M2))
 196 *
 197 *     Generate random m-by-n matrix C and a copy CF
 198 *
 199       DO J=1,N
 200          CALL ZLARNV( 2, ISEED, M2, C( 1, J ) )
 201       END DO
 202       CNORM = ZLANGE( '1', M2, N, C, M2, RWORK)
 203       CALL ZLACPY( 'Full', M2, N, C, M2, CF, M2 )
 204 *
 205 *     Apply Q to C as Q*C
 206 *
 207       CALL ZTPMQRT( 'L','N', M,N,K,L,NB,AF(NP1,1),M2,T,LDT,CF,M2,
 208      $               CF(NP1,1),M2,WORK,INFO)
 209 *
 210 *     Compute |Q*C - Q*C| / |C|
 211 *
 212       CALL ZGEMM( 'N', 'N', M2, N, M2, -ONE, Q, M2, C, M2, ONE, CF, M2 )
 213       RESID = ZLANGE( '1', M2, N, CF, M2, RWORK )
 214       IF( CNORM.GT.ZERO ) THEN
 215          RESULT( 3 ) = RESID / (EPS*MAX(1,M2)*CNORM)
 216       ELSE
 217          RESULT( 3 ) = ZERO
 218       END IF
 219 *
 220 *     Copy C into CF again
 221 *
 222       CALL ZLACPY( 'Full', M2, N, C, M2, CF, M2 )
 223 *
 224 *     Apply Q to C as QT*C
 225 *
 226       CALL ZTPMQRT( 'L','C',M,N,K,L,NB,AF(NP1,1),M2,T,LDT,CF,M2,
 227      $              CF(NP1,1),M2,WORK,INFO)
 228 *
 229 *     Compute |QT*C - QT*C| / |C|
 230 *
 231       CALL ZGEMM('C','N',M2,N,M2,-ONE,Q,M2,C,M2,ONE,CF,M2)
 232       RESID = ZLANGE( '1', M2, N, CF, M2, RWORK )
 233       IF( CNORM.GT.ZERO ) THEN
 234          RESULT( 4 ) = RESID / (EPS*MAX(1,M2)*CNORM)
 235       ELSE
 236          RESULT( 4 ) = ZERO
 237       END IF
 238 *
 239 *     Generate random n-by-m matrix D and a copy DF
 240 *
 241       DO J=1,M2
 242          CALL ZLARNV( 2, ISEED, N, D( 1, J ) )
 243       END DO
 244       DNORM = ZLANGE( '1', N, M2, D, N, RWORK)
 245       CALL ZLACPY( 'Full', N, M2, D, N, DF, N )
 246 *
 247 *     Apply Q to D as D*Q
 248 *
 249       CALL ZTPMQRT('R','N',N,M,N,L,NB,AF(NP1,1),M2,T,LDT,DF,N,
 250      $             DF(1,NP1),N,WORK,INFO)
 251 *
 252 *     Compute |D*Q - D*Q| / |D|
 253 *
 254       CALL ZGEMM('N','N',N,M2,M2,-ONE,D,N,Q,M2,ONE,DF,N)
 255       RESID = ZLANGE('1',N, M2,DF,N,RWORK )
 256       IF( CNORM.GT.ZERO ) THEN
 257          RESULT( 5 ) = RESID / (EPS*MAX(1,M2)*DNORM)
 258       ELSE
 259          RESULT( 5 ) = ZERO
 260       END IF
 261 *
 262 *     Copy D into DF again
 263 *
 264       CALL ZLACPY('Full',N,M2,D,N,DF,N )
 265 *
 266 *     Apply Q to D as D*QT
 267 *
 268       CALL ZTPMQRT('R','C',N,M,N,L,NB,AF(NP1,1),M2,T,LDT,DF,N,
 269      $             DF(1,NP1),N,WORK,INFO)
 270
 271 *
 272 *     Compute |D*QT - D*QT| / |D|
 273 *
 274       CALL ZGEMM( 'N', 'C', N, M2, M2, -ONE, D, N, Q, M2, ONE, DF, N )
 275       RESID = ZLANGE( '1', N, M2, DF, N, RWORK )
 276       IF( CNORM.GT.ZERO ) THEN
 277          RESULT( 6 ) = RESID / (EPS*MAX(1,M2)*DNORM)
 278       ELSE
 279          RESULT( 6 ) = ZERO
 280       END IF
 281 *
 282 *     Deallocate all arrays
 283 *
 284       DEALLOCATE ( A, AF, Q, R, RWORK, WORK, T, C, D, CF, DF)
 285       RETURN
 286       END
 287