LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zgbcon.f
Go to the documentation of this file.
00001 *> \brief \b ZGBCON
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZGBCON + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgbcon.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgbcon.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgbcon.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
00022 *                          WORK, RWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          NORM
00026 *       INTEGER            INFO, KL, KU, LDAB, N
00027 *       DOUBLE PRECISION   ANORM, RCOND
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       INTEGER            IPIV( * )
00031 *       DOUBLE PRECISION   RWORK( * )
00032 *       COMPLEX*16         AB( LDAB, * ), WORK( * )
00033 *       ..
00034 *  
00035 *
00036 *> \par Purpose:
00037 *  =============
00038 *>
00039 *> \verbatim
00040 *>
00041 *> ZGBCON estimates the reciprocal of the condition number of a complex
00042 *> general band matrix A, in either the 1-norm or the infinity-norm,
00043 *> using the LU factorization computed by ZGBTRF.
00044 *>
00045 *> An estimate is obtained for norm(inv(A)), and the reciprocal of the
00046 *> condition number is computed as
00047 *>    RCOND = 1 / ( norm(A) * norm(inv(A)) ).
00048 *> \endverbatim
00049 *
00050 *  Arguments:
00051 *  ==========
00052 *
00053 *> \param[in] NORM
00054 *> \verbatim
00055 *>          NORM is CHARACTER*1
00056 *>          Specifies whether the 1-norm condition number or the
00057 *>          infinity-norm condition number is required:
00058 *>          = '1' or 'O':  1-norm;
00059 *>          = 'I':         Infinity-norm.
00060 *> \endverbatim
00061 *>
00062 *> \param[in] N
00063 *> \verbatim
00064 *>          N is INTEGER
00065 *>          The order of the matrix A.  N >= 0.
00066 *> \endverbatim
00067 *>
00068 *> \param[in] KL
00069 *> \verbatim
00070 *>          KL is INTEGER
00071 *>          The number of subdiagonals within the band of A.  KL >= 0.
00072 *> \endverbatim
00073 *>
00074 *> \param[in] KU
00075 *> \verbatim
00076 *>          KU is INTEGER
00077 *>          The number of superdiagonals within the band of A.  KU >= 0.
00078 *> \endverbatim
00079 *>
00080 *> \param[in] AB
00081 *> \verbatim
00082 *>          AB is COMPLEX*16 array, dimension (LDAB,N)
00083 *>          Details of the LU factorization of the band matrix A, as
00084 *>          computed by ZGBTRF.  U is stored as an upper triangular band
00085 *>          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
00086 *>          the multipliers used during the factorization are stored in
00087 *>          rows KL+KU+2 to 2*KL+KU+1.
00088 *> \endverbatim
00089 *>
00090 *> \param[in] LDAB
00091 *> \verbatim
00092 *>          LDAB is INTEGER
00093 *>          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
00094 *> \endverbatim
00095 *>
00096 *> \param[in] IPIV
00097 *> \verbatim
00098 *>          IPIV is INTEGER array, dimension (N)
00099 *>          The pivot indices; for 1 <= i <= N, row i of the matrix was
00100 *>          interchanged with row IPIV(i).
00101 *> \endverbatim
00102 *>
00103 *> \param[in] ANORM
00104 *> \verbatim
00105 *>          ANORM is DOUBLE PRECISION
00106 *>          If NORM = '1' or 'O', the 1-norm of the original matrix A.
00107 *>          If NORM = 'I', the infinity-norm of the original matrix A.
00108 *> \endverbatim
00109 *>
00110 *> \param[out] RCOND
00111 *> \verbatim
00112 *>          RCOND is DOUBLE PRECISION
00113 *>          The reciprocal of the condition number of the matrix A,
00114 *>          computed as RCOND = 1/(norm(A) * norm(inv(A))).
00115 *> \endverbatim
00116 *>
00117 *> \param[out] WORK
00118 *> \verbatim
00119 *>          WORK is COMPLEX*16 array, dimension (2*N)
00120 *> \endverbatim
00121 *>
00122 *> \param[out] RWORK
00123 *> \verbatim
00124 *>          RWORK is DOUBLE PRECISION array, dimension (N)
00125 *> \endverbatim
00126 *>
00127 *> \param[out] INFO
00128 *> \verbatim
00129 *>          INFO is INTEGER
00130 *>          = 0:  successful exit
00131 *>          < 0: if INFO = -i, the i-th argument had an illegal value
00132 *> \endverbatim
00133 *
00134 *  Authors:
00135 *  ========
00136 *
00137 *> \author Univ. of Tennessee 
00138 *> \author Univ. of California Berkeley 
00139 *> \author Univ. of Colorado Denver 
00140 *> \author NAG Ltd. 
00141 *
00142 *> \date November 2011
00143 *
00144 *> \ingroup complex16GBcomputational
00145 *
00146 *  =====================================================================
00147       SUBROUTINE ZGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
00148      $                   WORK, RWORK, INFO )
00149 *
00150 *  -- LAPACK computational routine (version 3.4.0) --
00151 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00152 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00153 *     November 2011
00154 *
00155 *     .. Scalar Arguments ..
00156       CHARACTER          NORM
00157       INTEGER            INFO, KL, KU, LDAB, N
00158       DOUBLE PRECISION   ANORM, RCOND
00159 *     ..
00160 *     .. Array Arguments ..
00161       INTEGER            IPIV( * )
00162       DOUBLE PRECISION   RWORK( * )
00163       COMPLEX*16         AB( LDAB, * ), WORK( * )
00164 *     ..
00165 *
00166 *  =====================================================================
00167 *
00168 *     .. Parameters ..
00169       DOUBLE PRECISION   ONE, ZERO
00170       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00171 *     ..
00172 *     .. Local Scalars ..
00173       LOGICAL            LNOTI, ONENRM
00174       CHARACTER          NORMIN
00175       INTEGER            IX, J, JP, KASE, KASE1, KD, LM
00176       DOUBLE PRECISION   AINVNM, SCALE, SMLNUM
00177       COMPLEX*16         T, ZDUM
00178 *     ..
00179 *     .. Local Arrays ..
00180       INTEGER            ISAVE( 3 )
00181 *     ..
00182 *     .. External Functions ..
00183       LOGICAL            LSAME
00184       INTEGER            IZAMAX
00185       DOUBLE PRECISION   DLAMCH
00186       COMPLEX*16         ZDOTC
00187       EXTERNAL           LSAME, IZAMAX, DLAMCH, ZDOTC
00188 *     ..
00189 *     .. External Subroutines ..
00190       EXTERNAL           XERBLA, ZAXPY, ZDRSCL, ZLACN2, ZLATBS
00191 *     ..
00192 *     .. Intrinsic Functions ..
00193       INTRINSIC          ABS, DBLE, DIMAG, MIN
00194 *     ..
00195 *     .. Statement Functions ..
00196       DOUBLE PRECISION   CABS1
00197 *     ..
00198 *     .. Statement Function definitions ..
00199       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
00200 *     ..
00201 *     .. Executable Statements ..
00202 *
00203 *     Test the input parameters.
00204 *
00205       INFO = 0
00206       ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
00207       IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
00208          INFO = -1
00209       ELSE IF( N.LT.0 ) THEN
00210          INFO = -2
00211       ELSE IF( KL.LT.0 ) THEN
00212          INFO = -3
00213       ELSE IF( KU.LT.0 ) THEN
00214          INFO = -4
00215       ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN
00216          INFO = -6
00217       ELSE IF( ANORM.LT.ZERO ) THEN
00218          INFO = -8
00219       END IF
00220       IF( INFO.NE.0 ) THEN
00221          CALL XERBLA( 'ZGBCON', -INFO )
00222          RETURN
00223       END IF
00224 *
00225 *     Quick return if possible
00226 *
00227       RCOND = ZERO
00228       IF( N.EQ.0 ) THEN
00229          RCOND = ONE
00230          RETURN
00231       ELSE IF( ANORM.EQ.ZERO ) THEN
00232          RETURN
00233       END IF
00234 *
00235       SMLNUM = DLAMCH( 'Safe minimum' )
00236 *
00237 *     Estimate the norm of inv(A).
00238 *
00239       AINVNM = ZERO
00240       NORMIN = 'N'
00241       IF( ONENRM ) THEN
00242          KASE1 = 1
00243       ELSE
00244          KASE1 = 2
00245       END IF
00246       KD = KL + KU + 1
00247       LNOTI = KL.GT.0
00248       KASE = 0
00249    10 CONTINUE
00250       CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
00251       IF( KASE.NE.0 ) THEN
00252          IF( KASE.EQ.KASE1 ) THEN
00253 *
00254 *           Multiply by inv(L).
00255 *
00256             IF( LNOTI ) THEN
00257                DO 20 J = 1, N - 1
00258                   LM = MIN( KL, N-J )
00259                   JP = IPIV( J )
00260                   T = WORK( JP )
00261                   IF( JP.NE.J ) THEN
00262                      WORK( JP ) = WORK( J )
00263                      WORK( J ) = T
00264                   END IF
00265                   CALL ZAXPY( LM, -T, AB( KD+1, J ), 1, WORK( J+1 ), 1 )
00266    20          CONTINUE
00267             END IF
00268 *
00269 *           Multiply by inv(U).
00270 *
00271             CALL ZLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
00272      $                   KL+KU, AB, LDAB, WORK, SCALE, RWORK, INFO )
00273          ELSE
00274 *
00275 *           Multiply by inv(U**H).
00276 *
00277             CALL ZLATBS( 'Upper', 'Conjugate transpose', 'Non-unit',
00278      $                   NORMIN, N, KL+KU, AB, LDAB, WORK, SCALE, RWORK,
00279      $                   INFO )
00280 *
00281 *           Multiply by inv(L**H).
00282 *
00283             IF( LNOTI ) THEN
00284                DO 30 J = N - 1, 1, -1
00285                   LM = MIN( KL, N-J )
00286                   WORK( J ) = WORK( J ) - ZDOTC( LM, AB( KD+1, J ), 1,
00287      $                        WORK( J+1 ), 1 )
00288                   JP = IPIV( J )
00289                   IF( JP.NE.J ) THEN
00290                      T = WORK( JP )
00291                      WORK( JP ) = WORK( J )
00292                      WORK( J ) = T
00293                   END IF
00294    30          CONTINUE
00295             END IF
00296          END IF
00297 *
00298 *        Divide X by 1/SCALE if doing so will not cause overflow.
00299 *
00300          NORMIN = 'Y'
00301          IF( SCALE.NE.ONE ) THEN
00302             IX = IZAMAX( N, WORK, 1 )
00303             IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
00304      $         GO TO 40
00305             CALL ZDRSCL( N, SCALE, WORK, 1 )
00306          END IF
00307          GO TO 10
00308       END IF
00309 *
00310 *     Compute the estimate of the reciprocal condition number.
00311 *
00312       IF( AINVNM.NE.ZERO )
00313      $   RCOND = ( ONE / AINVNM ) / ANORM
00314 *
00315    40 CONTINUE
00316       RETURN
00317 *
00318 *     End of ZGBCON
00319 *
00320       END
 All Files Functions