LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zgtt01.f
Go to the documentation of this file.
00001 *> \brief \b ZGTT01
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE ZGTT01( N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK,
00012 *                          LDWORK, RWORK, RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       INTEGER            LDWORK, N
00016 *       DOUBLE PRECISION   RESID
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       INTEGER            IPIV( * )
00020 *       DOUBLE PRECISION   RWORK( * )
00021 *       COMPLEX*16         D( * ), DF( * ), DL( * ), DLF( * ), DU( * ),
00022 *      $                   DU2( * ), DUF( * ), WORK( LDWORK, * )
00023 *       ..
00024 *  
00025 *
00026 *> \par Purpose:
00027 *  =============
00028 *>
00029 *> \verbatim
00030 *>
00031 *> ZGTT01 reconstructs a tridiagonal matrix A from its LU factorization
00032 *> and computes the residual
00033 *>    norm(L*U - A) / ( norm(A) * EPS ),
00034 *> where EPS is the machine epsilon.
00035 *> \endverbatim
00036 *
00037 *  Arguments:
00038 *  ==========
00039 *
00040 *> \param[in] N
00041 *> \verbatim
00042 *>          N is INTEGTER
00043 *>          The order of the matrix A.  N >= 0.
00044 *> \endverbatim
00045 *>
00046 *> \param[in] DL
00047 *> \verbatim
00048 *>          DL is COMPLEX*16 array, dimension (N-1)
00049 *>          The (n-1) sub-diagonal elements of A.
00050 *> \endverbatim
00051 *>
00052 *> \param[in] D
00053 *> \verbatim
00054 *>          D is COMPLEX*16 array, dimension (N)
00055 *>          The diagonal elements of A.
00056 *> \endverbatim
00057 *>
00058 *> \param[in] DU
00059 *> \verbatim
00060 *>          DU is COMPLEX*16 array, dimension (N-1)
00061 *>          The (n-1) super-diagonal elements of A.
00062 *> \endverbatim
00063 *>
00064 *> \param[in] DLF
00065 *> \verbatim
00066 *>          DLF is COMPLEX*16 array, dimension (N-1)
00067 *>          The (n-1) multipliers that define the matrix L from the
00068 *>          LU factorization of A.
00069 *> \endverbatim
00070 *>
00071 *> \param[in] DF
00072 *> \verbatim
00073 *>          DF is COMPLEX*16 array, dimension (N)
00074 *>          The n diagonal elements of the upper triangular matrix U from
00075 *>          the LU factorization of A.
00076 *> \endverbatim
00077 *>
00078 *> \param[in] DUF
00079 *> \verbatim
00080 *>          DUF is COMPLEX*16 array, dimension (N-1)
00081 *>          The (n-1) elements of the first super-diagonal of U.
00082 *> \endverbatim
00083 *>
00084 *> \param[in] DU2
00085 *> \verbatim
00086 *>          DU2 is COMPLEX*16 array, dimension (N-2)
00087 *>          The (n-2) elements of the second super-diagonal of U.
00088 *> \endverbatim
00089 *>
00090 *> \param[in] IPIV
00091 *> \verbatim
00092 *>          IPIV is INTEGER array, dimension (N)
00093 *>          The pivot indices; for 1 <= i <= n, row i of the matrix was
00094 *>          interchanged with row IPIV(i).  IPIV(i) will always be either
00095 *>          i or i+1; IPIV(i) = i indicates a row interchange was not
00096 *>          required.
00097 *> \endverbatim
00098 *>
00099 *> \param[out] WORK
00100 *> \verbatim
00101 *>          WORK is COMPLEX*16 array, dimension (LDWORK,N)
00102 *> \endverbatim
00103 *>
00104 *> \param[in] LDWORK
00105 *> \verbatim
00106 *>          LDWORK is INTEGER
00107 *>          The leading dimension of the array WORK.  LDWORK >= max(1,N).
00108 *> \endverbatim
00109 *>
00110 *> \param[out] RWORK
00111 *> \verbatim
00112 *>          RWORK is DOUBLE PRECISION array, dimension (N)
00113 *> \endverbatim
00114 *>
00115 *> \param[out] RESID
00116 *> \verbatim
00117 *>          RESID is DOUBLE PRECISION
00118 *>          The scaled residual:  norm(L*U - A) / (norm(A) * EPS)
00119 *> \endverbatim
00120 *
00121 *  Authors:
00122 *  ========
00123 *
00124 *> \author Univ. of Tennessee 
00125 *> \author Univ. of California Berkeley 
00126 *> \author Univ. of Colorado Denver 
00127 *> \author NAG Ltd. 
00128 *
00129 *> \date November 2011
00130 *
00131 *> \ingroup complex16_lin
00132 *
00133 *  =====================================================================
00134       SUBROUTINE ZGTT01( N, DL, D, DU, DLF, DF, DUF, DU2, IPIV, WORK,
00135      $                   LDWORK, RWORK, RESID )
00136 *
00137 *  -- LAPACK test routine (version 3.4.0) --
00138 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00139 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00140 *     November 2011
00141 *
00142 *     .. Scalar Arguments ..
00143       INTEGER            LDWORK, N
00144       DOUBLE PRECISION   RESID
00145 *     ..
00146 *     .. Array Arguments ..
00147       INTEGER            IPIV( * )
00148       DOUBLE PRECISION   RWORK( * )
00149       COMPLEX*16         D( * ), DF( * ), DL( * ), DLF( * ), DU( * ),
00150      $                   DU2( * ), DUF( * ), WORK( LDWORK, * )
00151 *     ..
00152 *
00153 *  =====================================================================
00154 *
00155 *     .. Parameters ..
00156       DOUBLE PRECISION   ONE, ZERO
00157       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00158 *     ..
00159 *     .. Local Scalars ..
00160       INTEGER            I, IP, J, LASTJ
00161       DOUBLE PRECISION   ANORM, EPS
00162       COMPLEX*16         LI
00163 *     ..
00164 *     .. External Functions ..
00165       DOUBLE PRECISION   DLAMCH, ZLANGT, ZLANHS
00166       EXTERNAL           DLAMCH, ZLANGT, ZLANHS
00167 *     ..
00168 *     .. Intrinsic Functions ..
00169       INTRINSIC          MIN
00170 *     ..
00171 *     .. External Subroutines ..
00172       EXTERNAL           ZAXPY, ZSWAP
00173 *     ..
00174 *     .. Executable Statements ..
00175 *
00176 *     Quick return if possible
00177 *
00178       IF( N.LE.0 ) THEN
00179          RESID = ZERO
00180          RETURN
00181       END IF
00182 *
00183       EPS = DLAMCH( 'Epsilon' )
00184 *
00185 *     Copy the matrix U to WORK.
00186 *
00187       DO 20 J = 1, N
00188          DO 10 I = 1, N
00189             WORK( I, J ) = ZERO
00190    10    CONTINUE
00191    20 CONTINUE
00192       DO 30 I = 1, N
00193          IF( I.EQ.1 ) THEN
00194             WORK( I, I ) = DF( I )
00195             IF( N.GE.2 )
00196      $         WORK( I, I+1 ) = DUF( I )
00197             IF( N.GE.3 )
00198      $         WORK( I, I+2 ) = DU2( I )
00199          ELSE IF( I.EQ.N ) THEN
00200             WORK( I, I ) = DF( I )
00201          ELSE
00202             WORK( I, I ) = DF( I )
00203             WORK( I, I+1 ) = DUF( I )
00204             IF( I.LT.N-1 )
00205      $         WORK( I, I+2 ) = DU2( I )
00206          END IF
00207    30 CONTINUE
00208 *
00209 *     Multiply on the left by L.
00210 *
00211       LASTJ = N
00212       DO 40 I = N - 1, 1, -1
00213          LI = DLF( I )
00214          CALL ZAXPY( LASTJ-I+1, LI, WORK( I, I ), LDWORK,
00215      $               WORK( I+1, I ), LDWORK )
00216          IP = IPIV( I )
00217          IF( IP.EQ.I ) THEN
00218             LASTJ = MIN( I+2, N )
00219          ELSE
00220             CALL ZSWAP( LASTJ-I+1, WORK( I, I ), LDWORK, WORK( I+1, I ),
00221      $                  LDWORK )
00222          END IF
00223    40 CONTINUE
00224 *
00225 *     Subtract the matrix A.
00226 *
00227       WORK( 1, 1 ) = WORK( 1, 1 ) - D( 1 )
00228       IF( N.GT.1 ) THEN
00229          WORK( 1, 2 ) = WORK( 1, 2 ) - DU( 1 )
00230          WORK( N, N-1 ) = WORK( N, N-1 ) - DL( N-1 )
00231          WORK( N, N ) = WORK( N, N ) - D( N )
00232          DO 50 I = 2, N - 1
00233             WORK( I, I-1 ) = WORK( I, I-1 ) - DL( I-1 )
00234             WORK( I, I ) = WORK( I, I ) - D( I )
00235             WORK( I, I+1 ) = WORK( I, I+1 ) - DU( I )
00236    50    CONTINUE
00237       END IF
00238 *
00239 *     Compute the 1-norm of the tridiagonal matrix A.
00240 *
00241       ANORM = ZLANGT( '1', N, DL, D, DU )
00242 *
00243 *     Compute the 1-norm of WORK, which is only guaranteed to be
00244 *     upper Hessenberg.
00245 *
00246       RESID = ZLANHS( '1', N, WORK, LDWORK, RWORK )
00247 *
00248 *     Compute norm(L*U - A) / (norm(A) * EPS)
00249 *
00250       IF( ANORM.LE.ZERO ) THEN
00251          IF( RESID.NE.ZERO )
00252      $      RESID = ONE / EPS
00253       ELSE
00254          RESID = ( RESID / ANORM ) / EPS
00255       END IF
00256 *
00257       RETURN
00258 *
00259 *     End of ZGTT01
00260 *
00261       END
 All Files Functions