TESTING/EIG/sgrqts.f

   1 *> \brief \b SGRQTS
   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 SGRQTS( M, P, N, A, AF, Q, R, LDA, TAUA, B, BF, Z, T,
  12 *                          BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT )
  13 *
  14 *       .. Scalar Arguments ..
  15 *       INTEGER            LDA, LDB, LWORK, M, P, N
  16 *       ..
  17 *       .. Array Arguments ..
  18 *       REAL               A( LDA, * ), AF( LDA, * ), R( LDA, * ),
  19 *      $                   Q( LDA, * ),
  20 *      $                   B( LDB, * ), BF( LDB, * ), T( LDB, * ),
  21 *      $                   Z( LDB, * ), BWK( LDB, * ),
  22 *      $                   TAUA( * ), TAUB( * ),
  23 *      $                   RESULT( 4 ), RWORK( * ), WORK( LWORK )
  24 *       ..
  25 *
  26 *
  27 *> \par Purpose:
  28 *  =============
  29 *>
  30 *> \verbatim
  31 *>
  32 *> SGRQTS tests SGGRQF, which computes the GRQ factorization of an
  33 *> M-by-N matrix A and a P-by-N matrix B: A = R*Q and B = Z*T*Q.
  34 *> \endverbatim
  35 *
  36 *  Arguments:
  37 *  ==========
  38 *
  39 *> \param[in] M
  40 *> \verbatim
  41 *>          M is INTEGER
  42 *>          The number of rows of the matrix A.  M >= 0.
  43 *> \endverbatim
  44 *>
  45 *> \param[in] P
  46 *> \verbatim
  47 *>          P is INTEGER
  48 *>          The number of rows of the matrix B.  P >= 0.
  49 *> \endverbatim
  50 *>
  51 *> \param[in] N
  52 *> \verbatim
  53 *>          N is INTEGER
  54 *>          The number of columns of the matrices A and B.  N >= 0.
  55 *> \endverbatim
  56 *>
  57 *> \param[in] A
  58 *> \verbatim
  59 *>          A is REAL array, dimension (LDA,N)
  60 *>          The M-by-N matrix A.
  61 *> \endverbatim
  62 *>
  63 *> \param[out] AF
  64 *> \verbatim
  65 *>          AF is REAL array, dimension (LDA,N)
  66 *>          Details of the GRQ factorization of A and B, as returned
  67 *>          by SGGRQF, see SGGRQF for further details.
  68 *> \endverbatim
  69 *>
  70 *> \param[out] Q
  71 *> \verbatim
  72 *>          Q is REAL array, dimension (LDA,N)
  73 *>          The N-by-N orthogonal matrix Q.
  74 *> \endverbatim
  75 *>
  76 *> \param[out] R
  77 *> \verbatim
  78 *>          R is REAL array, dimension (LDA,MAX(M,N))
  79 *> \endverbatim
  80 *>
  81 *> \param[in] LDA
  82 *> \verbatim
  83 *>          LDA is INTEGER
  84 *>          The leading dimension of the arrays A, AF, R and Q.
  85 *>          LDA >= max(M,N).
  86 *> \endverbatim
  87 *>
  88 *> \param[out] TAUA
  89 *> \verbatim
  90 *>          TAUA is REAL array, dimension (min(M,N))
  91 *>          The scalar factors of the elementary reflectors, as returned
  92 *>          by SGGQRC.
  93 *> \endverbatim
  94 *>
  95 *> \param[in] B
  96 *> \verbatim
  97 *>          B is REAL array, dimension (LDB,N)
  98 *>          On entry, the P-by-N matrix A.
  99 *> \endverbatim
 100 *>
 101 *> \param[out] BF
 102 *> \verbatim
 103 *>          BF is REAL array, dimension (LDB,N)
 104 *>          Details of the GQR factorization of A and B, as returned
 105 *>          by SGGRQF, see SGGRQF for further details.
 106 *> \endverbatim
 107 *>
 108 *> \param[out] Z
 109 *> \verbatim
 110 *>          Z is REAL array, dimension (LDB,P)
 111 *>          The P-by-P orthogonal matrix Z.
 112 *> \endverbatim
 113 *>
 114 *> \param[out] T
 115 *> \verbatim
 116 *>          T is REAL array, dimension (LDB,max(P,N))
 117 *> \endverbatim
 118 *>
 119 *> \param[out] BWK
 120 *> \verbatim
 121 *>          BWK is REAL array, dimension (LDB,N)
 122 *> \endverbatim
 123 *>
 124 *> \param[in] LDB
 125 *> \verbatim
 126 *>          LDB is INTEGER
 127 *>          The leading dimension of the arrays B, BF, Z and T.
 128 *>          LDB >= max(P,N).
 129 *> \endverbatim
 130 *>
 131 *> \param[out] TAUB
 132 *> \verbatim
 133 *>          TAUB is REAL array, dimension (min(P,N))
 134 *>          The scalar factors of the elementary reflectors, as returned
 135 *>          by SGGRQF.
 136 *> \endverbatim
 137 *>
 138 *> \param[out] WORK
 139 *> \verbatim
 140 *>          WORK is REAL array, dimension (LWORK)
 141 *> \endverbatim
 142 *>
 143 *> \param[in] LWORK
 144 *> \verbatim
 145 *>          LWORK is INTEGER
 146 *>          The dimension of the array WORK, LWORK >= max(M,P,N)**2.
 147 *> \endverbatim
 148 *>
 149 *> \param[out] RWORK
 150 *> \verbatim
 151 *>          RWORK is REAL array, dimension (M)
 152 *> \endverbatim
 153 *>
 154 *> \param[out] RESULT
 155 *> \verbatim
 156 *>          RESULT is REAL array, dimension (4)
 157 *>          The test ratios:
 158 *>            RESULT(1) = norm( R - A*Q' ) / ( MAX(M,N)*norm(A)*ULP)
 159 *>            RESULT(2) = norm( T*Q - Z'*B ) / (MAX(P,N)*norm(B)*ULP)
 160 *>            RESULT(3) = norm( I - Q'*Q ) / ( N*ULP )
 161 *>            RESULT(4) = norm( I - Z'*Z ) / ( P*ULP )
 162 *> \endverbatim
 163 *
 164 *  Authors:
 165 *  ========
 166 *
 167 *> \author Univ. of Tennessee
 168 *> \author Univ. of California Berkeley
 169 *> \author Univ. of Colorado Denver
 170 *> \author NAG Ltd.
 171 *
 172 *> \date November 2011
 173 *
 174 *> \ingroup single_eig
 175 *
 176 *  =====================================================================
 177       SUBROUTINE SGRQTS( M, P, N, A, AF, Q, R, LDA, TAUA, B, BF, Z, T,
 178      $                   BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT )
 179 *
 180 *  -- LAPACK test routine (version 3.4.0) --
 181 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 182 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 183 *     November 2011
 184 *
 185 *     .. Scalar Arguments ..
 186       INTEGER            LDA, LDB, LWORK, M, P, N
 187 *     ..
 188 *     .. Array Arguments ..
 189       REAL               A( LDA, * ), AF( LDA, * ), R( LDA, * ),
 190      $                   Q( LDA, * ),
 191      $                   B( LDB, * ), BF( LDB, * ), T( LDB, * ),
 192      $                   Z( LDB, * ), BWK( LDB, * ),
 193      $                   TAUA( * ), TAUB( * ),
 194      $                   RESULT( 4 ), RWORK( * ), WORK( LWORK )
 195 *     ..
 196 *
 197 *  =====================================================================
 198 *
 199 *     .. Parameters ..
 200       REAL               ZERO, ONE
 201       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 202       REAL               ROGUE
 203       PARAMETER          ( ROGUE = -1.0E+10 )
 204 *     ..
 205 *     .. Local Scalars ..
 206       INTEGER            INFO
 207       REAL               ANORM, BNORM, ULP, UNFL, RESID
 208 *     ..
 209 *     .. External Functions ..
 210       REAL               SLAMCH, SLANGE, SLANSY
 211       EXTERNAL           SLAMCH, SLANGE, SLANSY
 212 *     ..
 213 *     .. External Subroutines ..
 214       EXTERNAL           SGEMM, SGGRQF, SLACPY, SLASET, SORGQR,
 215      $                   SORGRQ, SSYRK
 216 *     ..
 217 *     .. Intrinsic Functions ..
 218       INTRINSIC          MAX, MIN, REAL
 219 *     ..
 220 *     .. Executable Statements ..
 221 *
 222       ULP = SLAMCH( 'Precision' )
 223       UNFL = SLAMCH( 'Safe minimum' )
 224 *
 225 *     Copy the matrix A to the array AF.
 226 *
 227       CALL SLACPY( 'Full', M, N, A, LDA, AF, LDA )
 228       CALL SLACPY( 'Full', P, N, B, LDB, BF, LDB )
 229 *
 230       ANORM = MAX( SLANGE( '1', M, N, A, LDA, RWORK ), UNFL )
 231       BNORM = MAX( SLANGE( '1', P, N, B, LDB, RWORK ), UNFL )
 232 *
 233 *     Factorize the matrices A and B in the arrays AF and BF.
 234 *
 235       CALL SGGRQF( M, P, N, AF, LDA, TAUA, BF, LDB, TAUB, WORK,
 236      $             LWORK, INFO )
 237 *
 238 *     Generate the N-by-N matrix Q
 239 *
 240       CALL SLASET( 'Full', N, N, ROGUE, ROGUE, Q, LDA )
 241       IF( M.LE.N ) THEN
 242          IF( M.GT.0 .AND. M.LT.N )
 243      $      CALL SLACPY( 'Full', M, N-M, AF, LDA, Q( N-M+1, 1 ), LDA )
 244          IF( M.GT.1 )
 245      $      CALL SLACPY( 'Lower', M-1, M-1, AF( 2, N-M+1 ), LDA,
 246      $                   Q( N-M+2, N-M+1 ), LDA )
 247       ELSE
 248          IF( N.GT.1 )
 249      $      CALL SLACPY( 'Lower', N-1, N-1, AF( M-N+2, 1 ), LDA,
 250      $                   Q( 2, 1 ), LDA )
 251       END IF
 252       CALL SORGRQ( N, N, MIN( M, N ), Q, LDA, TAUA, WORK, LWORK, INFO )
 253 *
 254 *     Generate the P-by-P matrix Z
 255 *
 256       CALL SLASET( 'Full', P, P, ROGUE, ROGUE, Z, LDB )
 257       IF( P.GT.1 )
 258      $   CALL SLACPY( 'Lower', P-1, N, BF( 2,1 ), LDB, Z( 2,1 ), LDB )
 259       CALL SORGQR( P, P, MIN( P,N ), Z, LDB, TAUB, WORK, LWORK, INFO )
 260 *
 261 *     Copy R
 262 *
 263       CALL SLASET( 'Full', M, N, ZERO, ZERO, R, LDA )
 264       IF( M.LE.N )THEN
 265          CALL SLACPY( 'Upper', M, M, AF( 1, N-M+1 ), LDA, R( 1, N-M+1 ),
 266      $                LDA )
 267       ELSE
 268          CALL SLACPY( 'Full', M-N, N, AF, LDA, R, LDA )
 269          CALL SLACPY( 'Upper', N, N, AF( M-N+1, 1 ), LDA, R( M-N+1, 1 ),
 270      $                LDA )
 271       END IF
 272 *
 273 *     Copy T
 274 *
 275       CALL SLASET( 'Full', P, N, ZERO, ZERO, T, LDB )
 276       CALL SLACPY( 'Upper', P, N, BF, LDB, T, LDB )
 277 *
 278 *     Compute R - A*Q'
 279 *
 280       CALL SGEMM( 'No transpose', 'Transpose', M, N, N, -ONE, A, LDA, Q,
 281      $            LDA, ONE, R, LDA )
 282 *
 283 *     Compute norm( R - A*Q' ) / ( MAX(M,N)*norm(A)*ULP ) .
 284 *
 285       RESID = SLANGE( '1', M, N, R, LDA, RWORK )
 286       IF( ANORM.GT.ZERO ) THEN
 287          RESULT( 1 ) = ( ( RESID / REAL(MAX(1,M,N) ) ) / ANORM ) / ULP
 288       ELSE
 289          RESULT( 1 ) = ZERO
 290       END IF
 291 *
 292 *     Compute T*Q - Z'*B
 293 *
 294       CALL SGEMM( 'Transpose', 'No transpose', P, N, P, ONE, Z, LDB, B,
 295      $            LDB, ZERO, BWK, LDB )
 296       CALL SGEMM( 'No transpose', 'No transpose', P, N, N, ONE, T, LDB,
 297      $            Q, LDA, -ONE, BWK, LDB )
 298 *
 299 *     Compute norm( T*Q - Z'*B ) / ( MAX(P,N)*norm(A)*ULP ) .
 300 *
 301       RESID = SLANGE( '1', P, N, BWK, LDB, RWORK )
 302       IF( BNORM.GT.ZERO ) THEN
 303          RESULT( 2 ) = ( ( RESID / REAL( MAX( 1,P,M ) ) )/BNORM ) / ULP
 304       ELSE
 305          RESULT( 2 ) = ZERO
 306       END IF
 307 *
 308 *     Compute I - Q*Q'
 309 *
 310       CALL SLASET( 'Full', N, N, ZERO, ONE, R, LDA )
 311       CALL SSYRK( 'Upper', 'No Transpose', N, N, -ONE, Q, LDA, ONE, R,
 312      $            LDA )
 313 *
 314 *     Compute norm( I - Q'*Q ) / ( N * ULP ) .
 315 *
 316       RESID = SLANSY( '1', 'Upper', N, R, LDA, RWORK )
 317       RESULT( 3 ) = ( RESID / REAL( MAX( 1,N ) ) ) / ULP
 318 *
 319 *     Compute I - Z'*Z
 320 *
 321       CALL SLASET( 'Full', P, P, ZERO, ONE, T, LDB )
 322       CALL SSYRK( 'Upper', 'Transpose', P, P, -ONE, Z, LDB, ONE, T,
 323      $            LDB )
 324 *
 325 *     Compute norm( I - Z'*Z ) / ( P*ULP ) .
 326 *
 327       RESID = SLANSY( '1', 'Upper', P, T, LDB, RWORK )
 328       RESULT( 4 ) = ( RESID / REAL( MAX( 1,P ) ) ) / ULP
 329 *
 330       RETURN
 331 *
 332 *     End of SGRQTS
 333 *
 334       END