TESTING/LIN/sgtt01.f

   1 *> \brief \b SGTT01
   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 SGTT01( N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK,
  12 *                          LDWORK, RWORK, RESID )
  13 *
  14 *       .. Scalar Arguments ..
  15 *       INTEGER            LDWORK, N
  16 *       REAL               RESID
  17 *       ..
  18 *       .. Array Arguments ..
  19 *       INTEGER            IPIV( * )
  20 *       REAL               D( * ), DF( * ), DL( * ), DLF( * ), DU( * ),
  21 *      $                   DU2( * ), DUF( * ), RWORK( * ),
  22 *      $                   WORK( LDWORK, * )
  23 *       ..
  24 *
  25 *
  26 *> \par Purpose:
  27 *  =============
  28 *>
  29 *> \verbatim
  30 *>
  31 *> SGTT01 reconstructs a tridiagonal matrix A from its LU factorization
  32 *> and computes the residual
  33 *>    norm(L*U - A) / ( norm(A) * EPS ),
  34 *> where EPS is the machine epsilon.
  35 *> \endverbatim
  36 *
  37 *  Arguments:
  38 *  ==========
  39 *
  40 *> \param[in] N
  41 *> \verbatim
  42 *>          N is INTEGTER
  43 *>          The order of the matrix A.  N >= 0.
  44 *> \endverbatim
  45 *>
  46 *> \param[in] DL
  47 *> \verbatim
  48 *>          DL is REAL array, dimension (N-1)
  49 *>          The (n-1) sub-diagonal elements of A.
  50 *> \endverbatim
  51 *>
  52 *> \param[in] D
  53 *> \verbatim
  54 *>          D is REAL array, dimension (N)
  55 *>          The diagonal elements of A.
  56 *> \endverbatim
  57 *>
  58 *> \param[in] DU
  59 *> \verbatim
  60 *>          DU is REAL array, dimension (N-1)
  61 *>          The (n-1) super-diagonal elements of A.
  62 *> \endverbatim
  63 *>
  64 *> \param[in] DLF
  65 *> \verbatim
  66 *>          DLF is REAL array, dimension (N-1)
  67 *>          The (n-1) multipliers that define the matrix L from the
  68 *>          LU factorization of A.
  69 *> \endverbatim
  70 *>
  71 *> \param[in] DF
  72 *> \verbatim
  73 *>          DF is REAL array, dimension (N)
  74 *>          The n diagonal elements of the upper triangular matrix U from
  75 *>          the LU factorization of A.
  76 *> \endverbatim
  77 *>
  78 *> \param[in] DUF
  79 *> \verbatim
  80 *>          DUF is REAL array, dimension (N-1)
  81 *>          The (n-1) elements of the first super-diagonal of U.
  82 *> \endverbatim
  83 *>
  84 *> \param[in] DU2
  85 *> \verbatim
  86 *>          DU2 is REAL array, dimension (N-2)
  87 *>          The (n-2) elements of the second super-diagonal of U.
  88 *> \endverbatim
  89 *>
  90 *> \param[in] IPIV
  91 *> \verbatim
  92 *>          IPIV is INTEGER array, dimension (N)
  93 *>          The pivot indices; for 1 <= i <= n, row i of the matrix was
  94 *>          interchanged with row IPIV(i).  IPIV(i) will always be either
  95 *>          i or i+1; IPIV(i) = i indicates a row interchange was not
  96 *>          required.
  97 *> \endverbatim
  98 *>
  99 *> \param[out] WORK
 100 *> \verbatim
 101 *>          WORK is REAL array, dimension (LDWORK,N)
 102 *> \endverbatim
 103 *>
 104 *> \param[in] LDWORK
 105 *> \verbatim
 106 *>          LDWORK is INTEGER
 107 *>          The leading dimension of the array WORK.  LDWORK >= max(1,N).
 108 *> \endverbatim
 109 *>
 110 *> \param[out] RWORK
 111 *> \verbatim
 112 *>          RWORK is REAL array, dimension (N)
 113 *> \endverbatim
 114 *>
 115 *> \param[out] RESID
 116 *> \verbatim
 117 *>          RESID is REAL
 118 *>          The scaled residual:  norm(L*U - A) / (norm(A) * EPS)
 119 *> \endverbatim
 120 *
 121 *  Authors:
 122 *  ========
 123 *
 124 *> \author Univ. of Tennessee
 125 *> \author Univ. of California Berkeley
 126 *> \author Univ. of Colorado Denver
 127 *> \author NAG Ltd.
 128 *
 129 *> \date November 2011
 130 *
 131 *> \ingroup single_lin
 132 *
 133 *  =====================================================================
 134       SUBROUTINE SGTT01( N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK,
 135      $                   LDWORK, RWORK, RESID )
 136 *
 137 *  -- LAPACK test routine (version 3.4.0) --
 138 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 139 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 140 *     November 2011
 141 *
 142 *     .. Scalar Arguments ..
 143       INTEGER            LDWORK, N
 144       REAL               RESID
 145 *     ..
 146 *     .. Array Arguments ..
 147       INTEGER            IPIV( * )
 148       REAL               D( * ), DF( * ), DL( * ), DLF( * ), DU( * ),
 149      $                   DU2( * ), DUF( * ), RWORK( * ),
 150      $                   WORK( LDWORK, * )
 151 *     ..
 152 *
 153 *  =====================================================================
 154 *
 155 *     .. Parameters ..
 156       REAL               ONE, ZERO
 157       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
 158 *     ..
 159 *     .. Local Scalars ..
 160       INTEGER            I, IP, J, LASTJ
 161       REAL               ANORM, EPS, LI
 162 *     ..
 163 *     .. External Functions ..
 164       REAL               SLAMCH, SLANGT, SLANHS
 165       EXTERNAL           SLAMCH, SLANGT, SLANHS
 166 *     ..
 167 *     .. Intrinsic Functions ..
 168       INTRINSIC          MIN
 169 *     ..
 170 *     .. External Subroutines ..
 171       EXTERNAL           SAXPY, SSWAP
 172 *     ..
 173 *     .. Executable Statements ..
 174 *
 175 *     Quick return if possible
 176 *
 177       IF( N.LE.0 ) THEN
 178          RESID = ZERO
 179          RETURN
 180       END IF
 181 *
 182       EPS = SLAMCH( 'Epsilon' )
 183 *
 184 *     Copy the matrix U to WORK.
 185 *
 186       DO 20 J = 1, N
 187          DO 10 I = 1, N
 188             WORK( I, J ) = ZERO
 189    10    CONTINUE
 190    20 CONTINUE
 191       DO 30 I = 1, N
 192          IF( I.EQ.1 ) THEN
 193             WORK( I, I ) = DF( I )
 194             IF( N.GE.2 )
 195      $         WORK( I, I+1 ) = DUF( I )
 196             IF( N.GE.3 )
 197      $         WORK( I, I+2 ) = DU2( I )
 198          ELSE IF( I.EQ.N ) THEN
 199             WORK( I, I ) = DF( I )
 200          ELSE
 201             WORK( I, I ) = DF( I )
 202             WORK( I, I+1 ) = DUF( I )
 203             IF( I.LT.N-1 )
 204      $         WORK( I, I+2 ) = DU2( I )
 205          END IF
 206    30 CONTINUE
 207 *
 208 *     Multiply on the left by L.
 209 *
 210       LASTJ = N
 211       DO 40 I = N - 1, 1, -1
 212          LI = DLF( I )
 213          CALL SAXPY( LASTJ-I+1, LI, WORK( I, I ), LDWORK,
 214      $               WORK( I+1, I ), LDWORK )
 215          IP = IPIV( I )
 216          IF( IP.EQ.I ) THEN
 217             LASTJ = MIN( I+2, N )
 218          ELSE
 219             CALL SSWAP( LASTJ-I+1, WORK( I, I ), LDWORK, WORK( I+1, I ),
 220      $                  LDWORK )
 221          END IF
 222    40 CONTINUE
 223 *
 224 *     Subtract the matrix A.
 225 *
 226       WORK( 1, 1 ) = WORK( 1, 1 ) - D( 1 )
 227       IF( N.GT.1 ) THEN
 228          WORK( 1, 2 ) = WORK( 1, 2 ) - DU( 1 )
 229          WORK( N, N-1 ) = WORK( N, N-1 ) - DL( N-1 )
 230          WORK( N, N ) = WORK( N, N ) - D( N )
 231          DO 50 I = 2, N - 1
 232             WORK( I, I-1 ) = WORK( I, I-1 ) - DL( I-1 )
 233             WORK( I, I ) = WORK( I, I ) - D( I )
 234             WORK( I, I+1 ) = WORK( I, I+1 ) - DU( I )
 235    50    CONTINUE
 236       END IF
 237 *
 238 *     Compute the 1-norm of the tridiagonal matrix A.
 239 *
 240       ANORM = SLANGT( '1', N, DL, D, DU )
 241 *
 242 *     Compute the 1-norm of WORK, which is only guaranteed to be
 243 *     upper Hessenberg.
 244 *
 245       RESID = SLANHS( '1', N, WORK, LDWORK, RWORK )
 246 *
 247 *     Compute norm(L*U - A) / (norm(A) * EPS)
 248 *
 249       IF( ANORM.LE.ZERO ) THEN
 250          IF( RESID.NE.ZERO )
 251      $      RESID = ONE / EPS
 252       ELSE
 253          RESID = ( RESID / ANORM ) / EPS
 254       END IF
 255 *
 256       RETURN
 257 *
 258 *     End of SGTT01
 259 *
 260       END