LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dpbcon.f
Go to the documentation of this file.
00001 *> \brief \b DPBCON
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DPBCON + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dpbcon.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dpbcon.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dpbcon.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DPBCON( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK,
00022 *                          IWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          UPLO
00026 *       INTEGER            INFO, KD, LDAB, N
00027 *       DOUBLE PRECISION   ANORM, RCOND
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       INTEGER            IWORK( * )
00031 *       DOUBLE PRECISION   AB( LDAB, * ), WORK( * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> DPBCON estimates the reciprocal of the condition number (in the
00041 *> 1-norm) of a real symmetric positive definite band matrix using the
00042 *> Cholesky factorization A = U**T*U or A = L*L**T computed by DPBTRF.
00043 *>
00044 *> An estimate is obtained for norm(inv(A)), and the reciprocal of the
00045 *> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
00046 *> \endverbatim
00047 *
00048 *  Arguments:
00049 *  ==========
00050 *
00051 *> \param[in] UPLO
00052 *> \verbatim
00053 *>          UPLO is CHARACTER*1
00054 *>          = 'U':  Upper triangular factor stored in AB;
00055 *>          = 'L':  Lower triangular factor stored in AB.
00056 *> \endverbatim
00057 *>
00058 *> \param[in] N
00059 *> \verbatim
00060 *>          N is INTEGER
00061 *>          The order of the matrix A.  N >= 0.
00062 *> \endverbatim
00063 *>
00064 *> \param[in] KD
00065 *> \verbatim
00066 *>          KD is INTEGER
00067 *>          The number of superdiagonals of the matrix A if UPLO = 'U',
00068 *>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
00069 *> \endverbatim
00070 *>
00071 *> \param[in] AB
00072 *> \verbatim
00073 *>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
00074 *>          The triangular factor U or L from the Cholesky factorization
00075 *>          A = U**T*U or A = L*L**T of the band matrix A, stored in the
00076 *>          first KD+1 rows of the array.  The j-th column of U or L is
00077 *>          stored in the j-th column of the array AB as follows:
00078 *>          if UPLO ='U', AB(kd+1+i-j,j) = U(i,j) for max(1,j-kd)<=i<=j;
00079 *>          if UPLO ='L', AB(1+i-j,j)    = L(i,j) for j<=i<=min(n,j+kd).
00080 *> \endverbatim
00081 *>
00082 *> \param[in] LDAB
00083 *> \verbatim
00084 *>          LDAB is INTEGER
00085 *>          The leading dimension of the array AB.  LDAB >= KD+1.
00086 *> \endverbatim
00087 *>
00088 *> \param[in] ANORM
00089 *> \verbatim
00090 *>          ANORM is DOUBLE PRECISION
00091 *>          The 1-norm (or infinity-norm) of the symmetric band matrix A.
00092 *> \endverbatim
00093 *>
00094 *> \param[out] RCOND
00095 *> \verbatim
00096 *>          RCOND is DOUBLE PRECISION
00097 *>          The reciprocal of the condition number of the matrix A,
00098 *>          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
00099 *>          estimate of the 1-norm of inv(A) computed in this routine.
00100 *> \endverbatim
00101 *>
00102 *> \param[out] WORK
00103 *> \verbatim
00104 *>          WORK is DOUBLE PRECISION array, dimension (3*N)
00105 *> \endverbatim
00106 *>
00107 *> \param[out] IWORK
00108 *> \verbatim
00109 *>          IWORK is INTEGER array, dimension (N)
00110 *> \endverbatim
00111 *>
00112 *> \param[out] INFO
00113 *> \verbatim
00114 *>          INFO is INTEGER
00115 *>          = 0:  successful exit
00116 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00117 *> \endverbatim
00118 *
00119 *  Authors:
00120 *  ========
00121 *
00122 *> \author Univ. of Tennessee 
00123 *> \author Univ. of California Berkeley 
00124 *> \author Univ. of Colorado Denver 
00125 *> \author NAG Ltd. 
00126 *
00127 *> \date November 2011
00128 *
00129 *> \ingroup doubleOTHERcomputational
00130 *
00131 *  =====================================================================
00132       SUBROUTINE DPBCON( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK,
00133      $                   IWORK, INFO )
00134 *
00135 *  -- LAPACK computational routine (version 3.4.0) --
00136 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00137 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00138 *     November 2011
00139 *
00140 *     .. Scalar Arguments ..
00141       CHARACTER          UPLO
00142       INTEGER            INFO, KD, LDAB, N
00143       DOUBLE PRECISION   ANORM, RCOND
00144 *     ..
00145 *     .. Array Arguments ..
00146       INTEGER            IWORK( * )
00147       DOUBLE PRECISION   AB( LDAB, * ), WORK( * )
00148 *     ..
00149 *
00150 *  =====================================================================
00151 *
00152 *     .. Parameters ..
00153       DOUBLE PRECISION   ONE, ZERO
00154       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00155 *     ..
00156 *     .. Local Scalars ..
00157       LOGICAL            UPPER
00158       CHARACTER          NORMIN
00159       INTEGER            IX, KASE
00160       DOUBLE PRECISION   AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
00161 *     ..
00162 *     .. Local Arrays ..
00163       INTEGER            ISAVE( 3 )
00164 *     ..
00165 *     .. External Functions ..
00166       LOGICAL            LSAME
00167       INTEGER            IDAMAX
00168       DOUBLE PRECISION   DLAMCH
00169       EXTERNAL           LSAME, IDAMAX, DLAMCH
00170 *     ..
00171 *     .. External Subroutines ..
00172       EXTERNAL           DLACN2, DLATBS, DRSCL, XERBLA
00173 *     ..
00174 *     .. Intrinsic Functions ..
00175       INTRINSIC          ABS
00176 *     ..
00177 *     .. Executable Statements ..
00178 *
00179 *     Test the input parameters.
00180 *
00181       INFO = 0
00182       UPPER = LSAME( UPLO, 'U' )
00183       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00184          INFO = -1
00185       ELSE IF( N.LT.0 ) THEN
00186          INFO = -2
00187       ELSE IF( KD.LT.0 ) THEN
00188          INFO = -3
00189       ELSE IF( LDAB.LT.KD+1 ) THEN
00190          INFO = -5
00191       ELSE IF( ANORM.LT.ZERO ) THEN
00192          INFO = -6
00193       END IF
00194       IF( INFO.NE.0 ) THEN
00195          CALL XERBLA( 'DPBCON', -INFO )
00196          RETURN
00197       END IF
00198 *
00199 *     Quick return if possible
00200 *
00201       RCOND = ZERO
00202       IF( N.EQ.0 ) THEN
00203          RCOND = ONE
00204          RETURN
00205       ELSE IF( ANORM.EQ.ZERO ) THEN
00206          RETURN
00207       END IF
00208 *
00209       SMLNUM = DLAMCH( 'Safe minimum' )
00210 *
00211 *     Estimate the 1-norm of the inverse.
00212 *
00213       KASE = 0
00214       NORMIN = 'N'
00215    10 CONTINUE
00216       CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
00217       IF( KASE.NE.0 ) THEN
00218          IF( UPPER ) THEN
00219 *
00220 *           Multiply by inv(U**T).
00221 *
00222             CALL DLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
00223      $                   KD, AB, LDAB, WORK, SCALEL, WORK( 2*N+1 ),
00224      $                   INFO )
00225             NORMIN = 'Y'
00226 *
00227 *           Multiply by inv(U).
00228 *
00229             CALL DLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
00230      $                   KD, AB, LDAB, WORK, SCALEU, WORK( 2*N+1 ),
00231      $                   INFO )
00232          ELSE
00233 *
00234 *           Multiply by inv(L).
00235 *
00236             CALL DLATBS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
00237      $                   KD, AB, LDAB, WORK, SCALEL, WORK( 2*N+1 ),
00238      $                   INFO )
00239             NORMIN = 'Y'
00240 *
00241 *           Multiply by inv(L**T).
00242 *
00243             CALL DLATBS( 'Lower', 'Transpose', 'Non-unit', NORMIN, N,
00244      $                   KD, AB, LDAB, WORK, SCALEU, WORK( 2*N+1 ),
00245      $                   INFO )
00246          END IF
00247 *
00248 *        Multiply by 1/SCALE if doing so will not cause overflow.
00249 *
00250          SCALE = SCALEL*SCALEU
00251          IF( SCALE.NE.ONE ) THEN
00252             IX = IDAMAX( N, WORK, 1 )
00253             IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
00254      $         GO TO 20
00255             CALL DRSCL( N, SCALE, WORK, 1 )
00256          END IF
00257          GO TO 10
00258       END IF
00259 *
00260 *     Compute the estimate of the reciprocal condition number.
00261 *
00262       IF( AINVNM.NE.ZERO )
00263      $   RCOND = ( ONE / AINVNM ) / ANORM
00264 *
00265    20 CONTINUE
00266 *
00267       RETURN
00268 *
00269 *     End of DPBCON
00270 *
00271       END
 All Files Functions