LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dgbt01.f
Go to the documentation of this file.
00001 *> \brief \b DGBT01
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 DGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK,
00012 *                          RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       INTEGER            KL, KU, LDA, LDAFAC, M, N
00016 *       DOUBLE PRECISION   RESID
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       INTEGER            IPIV( * )
00020 *       DOUBLE PRECISION   A( LDA, * ), AFAC( LDAFAC, * ), WORK( * )
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *> DGBT01 reconstructs a band matrix  A  from its L*U factorization and
00030 *> computes the residual:
00031 *>    norm(L*U - A) / ( N * norm(A) * EPS ),
00032 *> where EPS is the machine epsilon.
00033 *>
00034 *> The expression L*U - A is computed one column at a time, so A and
00035 *> AFAC are not modified.
00036 *> \endverbatim
00037 *
00038 *  Arguments:
00039 *  ==========
00040 *
00041 *> \param[in] M
00042 *> \verbatim
00043 *>          M is INTEGER
00044 *>          The number of rows of the matrix A.  M >= 0.
00045 *> \endverbatim
00046 *>
00047 *> \param[in] N
00048 *> \verbatim
00049 *>          N is INTEGER
00050 *>          The number of columns of the matrix A.  N >= 0.
00051 *> \endverbatim
00052 *>
00053 *> \param[in] KL
00054 *> \verbatim
00055 *>          KL is INTEGER
00056 *>          The number of subdiagonals within the band of A.  KL >= 0.
00057 *> \endverbatim
00058 *>
00059 *> \param[in] KU
00060 *> \verbatim
00061 *>          KU is INTEGER
00062 *>          The number of superdiagonals within the band of A.  KU >= 0.
00063 *> \endverbatim
00064 *>
00065 *> \param[in,out] A
00066 *> \verbatim
00067 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
00068 *>          The original matrix A in band storage, stored in rows 1 to
00069 *>          KL+KU+1.
00070 *> \endverbatim
00071 *>
00072 *> \param[in] LDA
00073 *> \verbatim
00074 *>          LDA is INTEGER.
00075 *>          The leading dimension of the array A.  LDA >= max(1,KL+KU+1).
00076 *> \endverbatim
00077 *>
00078 *> \param[in] AFAC
00079 *> \verbatim
00080 *>          AFAC is DOUBLE PRECISION array, dimension (LDAFAC,N)
00081 *>          The factored form of the matrix A.  AFAC contains the banded
00082 *>          factors L and U from the L*U factorization, as computed by
00083 *>          DGBTRF.  U is stored as an upper triangular band matrix with
00084 *>          KL+KU superdiagonals in rows 1 to KL+KU+1, and the
00085 *>          multipliers used during the factorization are stored in rows
00086 *>          KL+KU+2 to 2*KL+KU+1.  See DGBTRF for further details.
00087 *> \endverbatim
00088 *>
00089 *> \param[in] LDAFAC
00090 *> \verbatim
00091 *>          LDAFAC is INTEGER
00092 *>          The leading dimension of the array AFAC.
00093 *>          LDAFAC >= max(1,2*KL*KU+1).
00094 *> \endverbatim
00095 *>
00096 *> \param[in] IPIV
00097 *> \verbatim
00098 *>          IPIV is INTEGER array, dimension (min(M,N))
00099 *>          The pivot indices from DGBTRF.
00100 *> \endverbatim
00101 *>
00102 *> \param[out] WORK
00103 *> \verbatim
00104 *>          WORK is DOUBLE PRECISION array, dimension (2*KL+KU+1)
00105 *> \endverbatim
00106 *>
00107 *> \param[out] RESID
00108 *> \verbatim
00109 *>          RESID is DOUBLE PRECISION
00110 *>          norm(L*U - A) / ( N * norm(A) * EPS )
00111 *> \endverbatim
00112 *
00113 *  Authors:
00114 *  ========
00115 *
00116 *> \author Univ. of Tennessee 
00117 *> \author Univ. of California Berkeley 
00118 *> \author Univ. of Colorado Denver 
00119 *> \author NAG Ltd. 
00120 *
00121 *> \date November 2011
00122 *
00123 *> \ingroup double_lin
00124 *
00125 *  =====================================================================
00126       SUBROUTINE DGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, IPIV, WORK,
00127      $                   RESID )
00128 *
00129 *  -- LAPACK test routine (version 3.4.0) --
00130 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00131 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00132 *     November 2011
00133 *
00134 *     .. Scalar Arguments ..
00135       INTEGER            KL, KU, LDA, LDAFAC, M, N
00136       DOUBLE PRECISION   RESID
00137 *     ..
00138 *     .. Array Arguments ..
00139       INTEGER            IPIV( * )
00140       DOUBLE PRECISION   A( LDA, * ), AFAC( LDAFAC, * ), WORK( * )
00141 *     ..
00142 *
00143 *  =====================================================================
00144 *
00145 *     .. Parameters ..
00146       DOUBLE PRECISION   ZERO, ONE
00147       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00148 *     ..
00149 *     .. Local Scalars ..
00150       INTEGER            I, I1, I2, IL, IP, IW, J, JL, JU, JUA, KD, LENJ
00151       DOUBLE PRECISION   ANORM, EPS, T
00152 *     ..
00153 *     .. External Functions ..
00154       DOUBLE PRECISION   DASUM, DLAMCH
00155       EXTERNAL           DASUM, DLAMCH
00156 *     ..
00157 *     .. External Subroutines ..
00158       EXTERNAL           DAXPY, DCOPY
00159 *     ..
00160 *     .. Intrinsic Functions ..
00161       INTRINSIC          DBLE, MAX, MIN
00162 *     ..
00163 *     .. Executable Statements ..
00164 *
00165 *     Quick exit if M = 0 or N = 0.
00166 *
00167       RESID = ZERO
00168       IF( M.LE.0 .OR. N.LE.0 )
00169      $   RETURN
00170 *
00171 *     Determine EPS and the norm of A.
00172 *
00173       EPS = DLAMCH( 'Epsilon' )
00174       KD = KU + 1
00175       ANORM = ZERO
00176       DO 10 J = 1, N
00177          I1 = MAX( KD+1-J, 1 )
00178          I2 = MIN( KD+M-J, KL+KD )
00179          IF( I2.GE.I1 )
00180      $      ANORM = MAX( ANORM, DASUM( I2-I1+1, A( I1, J ), 1 ) )
00181    10 CONTINUE
00182 *
00183 *     Compute one column at a time of L*U - A.
00184 *
00185       KD = KL + KU + 1
00186       DO 40 J = 1, N
00187 *
00188 *        Copy the J-th column of U to WORK.
00189 *
00190          JU = MIN( KL+KU, J-1 )
00191          JL = MIN( KL, M-J )
00192          LENJ = MIN( M, J ) - J + JU + 1
00193          IF( LENJ.GT.0 ) THEN
00194             CALL DCOPY( LENJ, AFAC( KD-JU, J ), 1, WORK, 1 )
00195             DO 20 I = LENJ + 1, JU + JL + 1
00196                WORK( I ) = ZERO
00197    20       CONTINUE
00198 *
00199 *           Multiply by the unit lower triangular matrix L.  Note that L
00200 *           is stored as a product of transformations and permutations.
00201 *
00202             DO 30 I = MIN( M-1, J ), J - JU, -1
00203                IL = MIN( KL, M-I )
00204                IF( IL.GT.0 ) THEN
00205                   IW = I - J + JU + 1
00206                   T = WORK( IW )
00207                   CALL DAXPY( IL, T, AFAC( KD+1, I ), 1, WORK( IW+1 ),
00208      $                        1 )
00209                   IP = IPIV( I )
00210                   IF( I.NE.IP ) THEN
00211                      IP = IP - J + JU + 1
00212                      WORK( IW ) = WORK( IP )
00213                      WORK( IP ) = T
00214                   END IF
00215                END IF
00216    30       CONTINUE
00217 *
00218 *           Subtract the corresponding column of A.
00219 *
00220             JUA = MIN( JU, KU )
00221             IF( JUA+JL+1.GT.0 )
00222      $         CALL DAXPY( JUA+JL+1, -ONE, A( KU+1-JUA, J ), 1,
00223      $                     WORK( JU+1-JUA ), 1 )
00224 *
00225 *           Compute the 1-norm of the column.
00226 *
00227             RESID = MAX( RESID, DASUM( JU+JL+1, WORK, 1 ) )
00228          END IF
00229    40 CONTINUE
00230 *
00231 *     Compute norm( L*U - A ) / ( N * norm(A) * EPS )
00232 *
00233       IF( ANORM.LE.ZERO ) THEN
00234          IF( RESID.NE.ZERO )
00235      $      RESID = ONE / EPS
00236       ELSE
00237          RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
00238       END IF
00239 *
00240       RETURN
00241 *
00242 *     End of DGBT01
00243 *
00244       END
 All Files Functions