TESTING/LIN/sqrt14.f

   1 *> \brief \b SQRT14
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       REAL             FUNCTION SQRT14( TRANS, M, N, NRHS, A, LDA, X,
  12 *                        LDX, WORK, LWORK )
  13 *
  14 *       .. Scalar Arguments ..
  15 *       CHARACTER          TRANS
  16 *       INTEGER            LDA, LDX, LWORK, M, N, NRHS
  17 *       ..
  18 *       .. Array Arguments ..
  19 *       REAL               A( LDA, * ), WORK( LWORK ), X( LDX, * )
  20 *       ..
  21 *
  22 *
  23 *> \par Purpose:
  24 *  =============
  25 *>
  26 *> \verbatim
  27 *>
  28 *> SQRT14 checks whether X is in the row space of A or A'.  It does so
  29 *> by scaling both X and A such that their norms are in the range
  30 *> [sqrt(eps), 1/sqrt(eps)], then computing a QR factorization of [A,X]
  31 *> (if TRANS = 'T') or an LQ factorization of [A',X]' (if TRANS = 'N'),
  32 *> and returning the norm of the trailing triangle, scaled by
  33 *> MAX(M,N,NRHS)*eps.
  34 *> \endverbatim
  35 *
  36 *  Arguments:
  37 *  ==========
  38 *
  39 *> \param[in] TRANS
  40 *> \verbatim
  41 *>          TRANS is CHARACTER*1
  42 *>          = 'N':  No transpose, check for X in the row space of A
  43 *>          = 'T':  Transpose, check for X in the row space of A'.
  44 *> \endverbatim
  45 *>
  46 *> \param[in] M
  47 *> \verbatim
  48 *>          M is INTEGER
  49 *>          The number of rows of the matrix A.
  50 *> \endverbatim
  51 *>
  52 *> \param[in] N
  53 *> \verbatim
  54 *>          N is INTEGER
  55 *>          The number of columns of the matrix A.
  56 *> \endverbatim
  57 *>
  58 *> \param[in] NRHS
  59 *> \verbatim
  60 *>          NRHS is INTEGER
  61 *>          The number of right hand sides, i.e., the number of columns
  62 *>          of X.
  63 *> \endverbatim
  64 *>
  65 *> \param[in] A
  66 *> \verbatim
  67 *>          A is REAL array, dimension (LDA,N)
  68 *>          The M-by-N matrix A.
  69 *> \endverbatim
  70 *>
  71 *> \param[in] LDA
  72 *> \verbatim
  73 *>          LDA is INTEGER
  74 *>          The leading dimension of the array A.
  75 *> \endverbatim
  76 *>
  77 *> \param[in] X
  78 *> \verbatim
  79 *>          X is REAL array, dimension (LDX,NRHS)
  80 *>          If TRANS = 'N', the N-by-NRHS matrix X.
  81 *>          IF TRANS = 'T', the M-by-NRHS matrix X.
  82 *> \endverbatim
  83 *>
  84 *> \param[in] LDX
  85 *> \verbatim
  86 *>          LDX is INTEGER
  87 *>          The leading dimension of the array X.
  88 *> \endverbatim
  89 *>
  90 *> \param[out] WORK
  91 *> \verbatim
  92 *>          WORK is REAL array dimension (LWORK)
  93 *> \endverbatim
  94 *>
  95 *> \param[in] LWORK
  96 *> \verbatim
  97 *>          LWORK is INTEGER
  98 *>          length of workspace array required
  99 *>          If TRANS = 'N', LWORK >= (M+NRHS)*(N+2);
 100 *>          if TRANS = 'T', LWORK >= (N+NRHS)*(M+2).
 101 *> \endverbatim
 102 *
 103 *  Authors:
 104 *  ========
 105 *
 106 *> \author Univ. of Tennessee
 107 *> \author Univ. of California Berkeley
 108 *> \author Univ. of Colorado Denver
 109 *> \author NAG Ltd.
 110 *
 111 *> \date November 2011
 112 *
 113 *> \ingroup single_lin
 114 *
 115 *  =====================================================================
 116       REAL             FUNCTION SQRT14( TRANS, M, N, NRHS, A, LDA, X,
 117      $                 LDX, WORK, LWORK )
 118 *
 119 *  -- LAPACK test routine (version 3.4.0) --
 120 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 121 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 122 *     November 2011
 123 *
 124 *     .. Scalar Arguments ..
 125       CHARACTER          TRANS
 126       INTEGER            LDA, LDX, LWORK, M, N, NRHS
 127 *     ..
 128 *     .. Array Arguments ..
 129       REAL               A( LDA, * ), WORK( LWORK ), X( LDX, * )
 130 *     ..
 131 *
 132 *  =====================================================================
 133 *
 134 *     .. Parameters ..
 135       REAL               ZERO, ONE
 136       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
 137 *     ..
 138 *     .. Local Scalars ..
 139       LOGICAL            TPSD
 140       INTEGER            I, INFO, J, LDWORK
 141       REAL               ANRM, ERR, XNRM
 142 *     ..
 143 *     .. Local Arrays ..
 144       REAL               RWORK( 1 )
 145 *     ..
 146 *     .. External Functions ..
 147       LOGICAL            LSAME
 148       REAL               SLAMCH, SLANGE
 149       EXTERNAL           LSAME, SLAMCH, SLANGE
 150 *     ..
 151 *     .. External Subroutines ..
 152       EXTERNAL           SGELQ2, SGEQR2, SLACPY, SLASCL, XERBLA
 153 *     ..
 154 *     .. Intrinsic Functions ..
 155       INTRINSIC          ABS, MAX, MIN, REAL
 156 *     ..
 157 *     .. Executable Statements ..
 158 *
 159       SQRT14 = ZERO
 160       IF( LSAME( TRANS, 'N' ) ) THEN
 161          LDWORK = M + NRHS
 162          TPSD = .FALSE.
 163          IF( LWORK.LT.( M+NRHS )*( N+2 ) ) THEN
 164             CALL XERBLA( 'SQRT14', 10 )
 165             RETURN
 166          ELSE IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
 167             RETURN
 168          END IF
 169       ELSE IF( LSAME( TRANS, 'T' ) ) THEN
 170          LDWORK = M
 171          TPSD = .TRUE.
 172          IF( LWORK.LT.( N+NRHS )*( M+2 ) ) THEN
 173             CALL XERBLA( 'SQRT14', 10 )
 174             RETURN
 175          ELSE IF( M.LE.0 .OR. NRHS.LE.0 ) THEN
 176             RETURN
 177          END IF
 178       ELSE
 179          CALL XERBLA( 'SQRT14', 1 )
 180          RETURN
 181       END IF
 182 *
 183 *     Copy and scale A
 184 *
 185       CALL SLACPY( 'All', M, N, A, LDA, WORK, LDWORK )
 186       ANRM = SLANGE( 'M', M, N, WORK, LDWORK, RWORK )
 187       IF( ANRM.NE.ZERO )
 188      $   CALL SLASCL( 'G', 0, 0, ANRM, ONE, M, N, WORK, LDWORK, INFO )
 189 *
 190 *     Copy X or X' into the right place and scale it
 191 *
 192       IF( TPSD ) THEN
 193 *
 194 *        Copy X into columns n+1:n+nrhs of work
 195 *
 196          CALL SLACPY( 'All', M, NRHS, X, LDX, WORK( N*LDWORK+1 ),
 197      $                LDWORK )
 198          XNRM = SLANGE( 'M', M, NRHS, WORK( N*LDWORK+1 ), LDWORK,
 199      $          RWORK )
 200          IF( XNRM.NE.ZERO )
 201      $      CALL SLASCL( 'G', 0, 0, XNRM, ONE, M, NRHS,
 202      $                   WORK( N*LDWORK+1 ), LDWORK, INFO )
 203          ANRM = SLANGE( 'One-norm', M, N+NRHS, WORK, LDWORK, RWORK )
 204 *
 205 *        Compute QR factorization of X
 206 *
 207          CALL SGEQR2( M, N+NRHS, WORK, LDWORK,
 208      $                WORK( LDWORK*( N+NRHS )+1 ),
 209      $                WORK( LDWORK*( N+NRHS )+MIN( M, N+NRHS )+1 ),
 210      $                INFO )
 211 *
 212 *        Compute largest entry in upper triangle of
 213 *        work(n+1:m,n+1:n+nrhs)
 214 *
 215          ERR = ZERO
 216          DO 20 J = N + 1, N + NRHS
 217             DO 10 I = N + 1, MIN( M, J )
 218                ERR = MAX( ERR, ABS( WORK( I+( J-1 )*M ) ) )
 219    10       CONTINUE
 220    20    CONTINUE
 221 *
 222       ELSE
 223 *
 224 *        Copy X' into rows m+1:m+nrhs of work
 225 *
 226          DO 40 I = 1, N
 227             DO 30 J = 1, NRHS
 228                WORK( M+J+( I-1 )*LDWORK ) = X( I, J )
 229    30       CONTINUE
 230    40    CONTINUE
 231 *
 232          XNRM = SLANGE( 'M', NRHS, N, WORK( M+1 ), LDWORK, RWORK )
 233          IF( XNRM.NE.ZERO )
 234      $      CALL SLASCL( 'G', 0, 0, XNRM, ONE, NRHS, N, WORK( M+1 ),
 235      $                   LDWORK, INFO )
 236 *
 237 *        Compute LQ factorization of work
 238 *
 239          CALL SGELQ2( LDWORK, N, WORK, LDWORK, WORK( LDWORK*N+1 ),
 240      $                WORK( LDWORK*( N+1 )+1 ), INFO )
 241 *
 242 *        Compute largest entry in lower triangle in
 243 *        work(m+1:m+nrhs,m+1:n)
 244 *
 245          ERR = ZERO
 246          DO 60 J = M + 1, N
 247             DO 50 I = J, LDWORK
 248                ERR = MAX( ERR, ABS( WORK( I+( J-1 )*LDWORK ) ) )
 249    50       CONTINUE
 250    60    CONTINUE
 251 *
 252       END IF
 253 *
 254       SQRT14 = ERR / ( REAL( MAX( M, N, NRHS ) )*SLAMCH( 'Epsilon' ) )
 255 *
 256       RETURN
 257 *
 258 *     End of SQRT14
 259 *
 260       END