LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dpbt01.f
Go to the documentation of this file.
00001 *> \brief \b DPBT01
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 DPBT01( UPLO, N, KD, A, LDA, AFAC, LDAFAC, RWORK,
00012 *                          RESID )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       CHARACTER          UPLO
00016 *       INTEGER            KD, LDA, LDAFAC, N
00017 *       DOUBLE PRECISION   RESID
00018 *       ..
00019 *       .. Array Arguments ..
00020 *       DOUBLE PRECISION   A( LDA, * ), AFAC( LDAFAC, * ), RWORK( * )
00021 *       ..
00022 *  
00023 *
00024 *> \par Purpose:
00025 *  =============
00026 *>
00027 *> \verbatim
00028 *>
00029 *> DPBT01 reconstructs a symmetric positive definite band matrix A from
00030 *> its L*L' or U'*U factorization and computes the residual
00031 *>    norm( L*L' - A ) / ( N * norm(A) * EPS ) or
00032 *>    norm( U'*U - A ) / ( N * norm(A) * EPS ),
00033 *> where EPS is the machine epsilon, L' is the conjugate transpose of
00034 *> L, and U' is the conjugate transpose of U.
00035 *> \endverbatim
00036 *
00037 *  Arguments:
00038 *  ==========
00039 *
00040 *> \param[in] UPLO
00041 *> \verbatim
00042 *>          UPLO is CHARACTER*1
00043 *>          Specifies whether the upper or lower triangular part of the
00044 *>          symmetric matrix A is stored:
00045 *>          = 'U':  Upper triangular
00046 *>          = 'L':  Lower triangular
00047 *> \endverbatim
00048 *>
00049 *> \param[in] N
00050 *> \verbatim
00051 *>          N is INTEGER
00052 *>          The number of rows and columns of the matrix A.  N >= 0.
00053 *> \endverbatim
00054 *>
00055 *> \param[in] KD
00056 *> \verbatim
00057 *>          KD is INTEGER
00058 *>          The number of super-diagonals of the matrix A if UPLO = 'U',
00059 *>          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0.
00060 *> \endverbatim
00061 *>
00062 *> \param[in] A
00063 *> \verbatim
00064 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
00065 *>          The original symmetric band matrix A.  If UPLO = 'U', the
00066 *>          upper triangular part of A is stored as a band matrix; if
00067 *>          UPLO = 'L', the lower triangular part of A is stored.  The
00068 *>          columns of the appropriate triangle are stored in the columns
00069 *>          of A and the diagonals of the triangle are stored in the rows
00070 *>          of A.  See DPBTRF for further details.
00071 *> \endverbatim
00072 *>
00073 *> \param[in] LDA
00074 *> \verbatim
00075 *>          LDA is INTEGER.
00076 *>          The leading dimension of the array A.  LDA >= max(1,KD+1).
00077 *> \endverbatim
00078 *>
00079 *> \param[in] AFAC
00080 *> \verbatim
00081 *>          AFAC is DOUBLE PRECISION array, dimension (LDAFAC,N)
00082 *>          The factored form of the matrix A.  AFAC contains the factor
00083 *>          L or U from the L*L' or U'*U factorization in band storage
00084 *>          format, as computed by DPBTRF.
00085 *> \endverbatim
00086 *>
00087 *> \param[in] LDAFAC
00088 *> \verbatim
00089 *>          LDAFAC is INTEGER
00090 *>          The leading dimension of the array AFAC.
00091 *>          LDAFAC >= max(1,KD+1).
00092 *> \endverbatim
00093 *>
00094 *> \param[out] RWORK
00095 *> \verbatim
00096 *>          RWORK is DOUBLE PRECISION array, dimension (N)
00097 *> \endverbatim
00098 *>
00099 *> \param[out] RESID
00100 *> \verbatim
00101 *>          RESID is DOUBLE PRECISION
00102 *>          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS )
00103 *>          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS )
00104 *> \endverbatim
00105 *
00106 *  Authors:
00107 *  ========
00108 *
00109 *> \author Univ. of Tennessee 
00110 *> \author Univ. of California Berkeley 
00111 *> \author Univ. of Colorado Denver 
00112 *> \author NAG Ltd. 
00113 *
00114 *> \date November 2011
00115 *
00116 *> \ingroup double_lin
00117 *
00118 *  =====================================================================
00119       SUBROUTINE DPBT01( UPLO, N, KD, A, LDA, AFAC, LDAFAC, RWORK,
00120      $                   RESID )
00121 *
00122 *  -- LAPACK test routine (version 3.4.0) --
00123 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00124 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00125 *     November 2011
00126 *
00127 *     .. Scalar Arguments ..
00128       CHARACTER          UPLO
00129       INTEGER            KD, LDA, LDAFAC, N
00130       DOUBLE PRECISION   RESID
00131 *     ..
00132 *     .. Array Arguments ..
00133       DOUBLE PRECISION   A( LDA, * ), AFAC( LDAFAC, * ), RWORK( * )
00134 *     ..
00135 *
00136 *  =====================================================================
00137 *
00138 *
00139 *     .. Parameters ..
00140       DOUBLE PRECISION   ZERO, ONE
00141       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00142 *     ..
00143 *     .. Local Scalars ..
00144       INTEGER            I, J, K, KC, KLEN, ML, MU
00145       DOUBLE PRECISION   ANORM, EPS, T
00146 *     ..
00147 *     .. External Functions ..
00148       LOGICAL            LSAME
00149       DOUBLE PRECISION   DDOT, DLAMCH, DLANSB
00150       EXTERNAL           LSAME, DDOT, DLAMCH, DLANSB
00151 *     ..
00152 *     .. External Subroutines ..
00153       EXTERNAL           DSCAL, DSYR, DTRMV
00154 *     ..
00155 *     .. Intrinsic Functions ..
00156       INTRINSIC          DBLE, MAX, MIN
00157 *     ..
00158 *     .. Executable Statements ..
00159 *
00160 *     Quick exit if N = 0.
00161 *
00162       IF( N.LE.0 ) THEN
00163          RESID = ZERO
00164          RETURN
00165       END IF
00166 *
00167 *     Exit with RESID = 1/EPS if ANORM = 0.
00168 *
00169       EPS = DLAMCH( 'Epsilon' )
00170       ANORM = DLANSB( '1', UPLO, N, KD, A, LDA, RWORK )
00171       IF( ANORM.LE.ZERO ) THEN
00172          RESID = ONE / EPS
00173          RETURN
00174       END IF
00175 *
00176 *     Compute the product U'*U, overwriting U.
00177 *
00178       IF( LSAME( UPLO, 'U' ) ) THEN
00179          DO 10 K = N, 1, -1
00180             KC = MAX( 1, KD+2-K )
00181             KLEN = KD + 1 - KC
00182 *
00183 *           Compute the (K,K) element of the result.
00184 *
00185             T = DDOT( KLEN+1, AFAC( KC, K ), 1, AFAC( KC, K ), 1 )
00186             AFAC( KD+1, K ) = T
00187 *
00188 *           Compute the rest of column K.
00189 *
00190             IF( KLEN.GT.0 )
00191      $         CALL DTRMV( 'Upper', 'Transpose', 'Non-unit', KLEN,
00192      $                     AFAC( KD+1, K-KLEN ), LDAFAC-1,
00193      $                     AFAC( KC, K ), 1 )
00194 *
00195    10    CONTINUE
00196 *
00197 *     UPLO = 'L':  Compute the product L*L', overwriting L.
00198 *
00199       ELSE
00200          DO 20 K = N, 1, -1
00201             KLEN = MIN( KD, N-K )
00202 *
00203 *           Add a multiple of column K of the factor L to each of
00204 *           columns K+1 through N.
00205 *
00206             IF( KLEN.GT.0 )
00207      $         CALL DSYR( 'Lower', KLEN, ONE, AFAC( 2, K ), 1,
00208      $                    AFAC( 1, K+1 ), LDAFAC-1 )
00209 *
00210 *           Scale column K by the diagonal element.
00211 *
00212             T = AFAC( 1, K )
00213             CALL DSCAL( KLEN+1, T, AFAC( 1, K ), 1 )
00214 *
00215    20    CONTINUE
00216       END IF
00217 *
00218 *     Compute the difference  L*L' - A  or  U'*U - A.
00219 *
00220       IF( LSAME( UPLO, 'U' ) ) THEN
00221          DO 40 J = 1, N
00222             MU = MAX( 1, KD+2-J )
00223             DO 30 I = MU, KD + 1
00224                AFAC( I, J ) = AFAC( I, J ) - A( I, J )
00225    30       CONTINUE
00226    40    CONTINUE
00227       ELSE
00228          DO 60 J = 1, N
00229             ML = MIN( KD+1, N-J+1 )
00230             DO 50 I = 1, ML
00231                AFAC( I, J ) = AFAC( I, J ) - A( I, J )
00232    50       CONTINUE
00233    60    CONTINUE
00234       END IF
00235 *
00236 *     Compute norm( L*L' - A ) / ( N * norm(A) * EPS )
00237 *
00238       RESID = DLANSB( 'I', UPLO, N, KD, AFAC, LDAFAC, RWORK )
00239 *
00240       RESID = ( ( RESID / DBLE( N ) ) / ANORM ) / EPS
00241 *
00242       RETURN
00243 *
00244 *     End of DPBT01
00245 *
00246       END
 All Files Functions