LAPACK: SRC/cgbcon.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
00002      $                   WORK, RWORK, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     Modified to call CLACN2 in place of CLACON, 10 Feb 03, SJH.
00010 *
00011 *     .. Scalar Arguments ..
00012       CHARACTER          NORM
00013       INTEGER            INFO, KL, KU, LDAB, N
00014       REAL               ANORM, RCOND
00015 *     ..
00016 *     .. Array Arguments ..
00017       INTEGER            IPIV( * )
00018       REAL               RWORK( * )
00019       COMPLEX            AB( LDAB, * ), WORK( * )
00020 *     ..
00021 *
00022 *  Purpose
00023 *  =======
00024 *
00025 *  CGBCON estimates the reciprocal of the condition number of a complex
00026 *  general band matrix A, in either the 1-norm or the infinity-norm,
00027 *  using the LU factorization computed by CGBTRF.
00028 *
00029 *  An estimate is obtained for norm(inv(A)), and the reciprocal of the
00030 *  condition number is computed as
00031 *     RCOND = 1 / ( norm(A) * norm(inv(A)) ).
00032 *
00033 *  Arguments
00034 *  =========
00035 *
00036 *  NORM    (input) CHARACTER*1
00037 *          Specifies whether the 1-norm condition number or the
00038 *          infinity-norm condition number is required:
00039 *          = '1' or 'O':  1-norm;
00040 *          = 'I':         Infinity-norm.
00041 *
00042 *  N       (input) INTEGER
00043 *          The order of the matrix A.  N >= 0.
00044 *
00045 *  KL      (input) INTEGER
00046 *          The number of subdiagonals within the band of A.  KL >= 0.
00047 *
00048 *  KU      (input) INTEGER
00049 *          The number of superdiagonals within the band of A.  KU >= 0.
00050 *
00051 *  AB      (input) COMPLEX array, dimension (LDAB,N)
00052 *          Details of the LU factorization of the band matrix A, as
00053 *          computed by CGBTRF.  U is stored as an upper triangular band
00054 *          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
00055 *          the multipliers used during the factorization are stored in
00056 *          rows KL+KU+2 to 2*KL+KU+1.
00057 *
00058 *  LDAB    (input) INTEGER
00059 *          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
00060 *
00061 *  IPIV    (input) INTEGER array, dimension (N)
00062 *          The pivot indices; for 1 <= i <= N, row i of the matrix was
00063 *          interchanged with row IPIV(i).
00064 *
00065 *  ANORM   (input) REAL
00066 *          If NORM = '1' or 'O', the 1-norm of the original matrix A.
00067 *          If NORM = 'I', the infinity-norm of the original matrix A.
00068 *
00069 *  RCOND   (output) REAL
00070 *          The reciprocal of the condition number of the matrix A,
00071 *          computed as RCOND = 1/(norm(A) * norm(inv(A))).
00072 *
00073 *  WORK    (workspace) COMPLEX array, dimension (2*N)
00074 *
00075 *  RWORK   (workspace) REAL array, dimension (N)
00076 *
00077 *  INFO    (output) INTEGER
00078 *          = 0:  successful exit
00079 *          < 0: if INFO = -i, the i-th argument had an illegal value
00080 *
00081 *  =====================================================================
00082 *
00083 *     .. Parameters ..
00084       REAL               ONE, ZERO
00085       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00086 *     ..
00087 *     .. Local Scalars ..
00088       LOGICAL            LNOTI, ONENRM
00089       CHARACTER          NORMIN
00090       INTEGER            IX, J, JP, KASE, KASE1, KD, LM
00091       REAL               AINVNM, SCALE, SMLNUM
00092       COMPLEX            T, ZDUM
00093 *     ..
00094 *     .. Local Arrays ..
00095       INTEGER            ISAVE( 3 )
00096 *     ..
00097 *     .. External Functions ..
00098       LOGICAL            LSAME
00099       INTEGER            ICAMAX
00100       REAL               SLAMCH
00101       COMPLEX            CDOTC
00102       EXTERNAL           LSAME, ICAMAX, SLAMCH, CDOTC
00103 *     ..
00104 *     .. External Subroutines ..
00105       EXTERNAL           CAXPY, CLACN2, CLATBS, CSRSCL, XERBLA
00106 *     ..
00107 *     .. Intrinsic Functions ..
00108       INTRINSIC          ABS, AIMAG, MIN, REAL
00109 *     ..
00110 *     .. Statement Functions ..
00111       REAL               CABS1
00112 *     ..
00113 *     .. Statement Function definitions ..
00114       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
00115 *     ..
00116 *     .. Executable Statements ..
00117 *
00118 *     Test the input parameters.
00119 *
00120       INFO = 0
00121       ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
00122       IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
00123          INFO = -1
00124       ELSE IF( N.LT.0 ) THEN
00125          INFO = -2
00126       ELSE IF( KL.LT.0 ) THEN
00127          INFO = -3
00128       ELSE IF( KU.LT.0 ) THEN
00129          INFO = -4
00130       ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN
00131          INFO = -6
00132       ELSE IF( ANORM.LT.ZERO ) THEN
00133          INFO = -8
00134       END IF
00135       IF( INFO.NE.0 ) THEN
00136          CALL XERBLA( 'CGBCON', -INFO )
00137          RETURN
00138       END IF
00139 *
00140 *     Quick return if possible
00141 *
00142       RCOND = ZERO
00143       IF( N.EQ.0 ) THEN
00144          RCOND = ONE
00145          RETURN
00146       ELSE IF( ANORM.EQ.ZERO ) THEN
00147          RETURN
00148       END IF
00149 *
00150       SMLNUM = SLAMCH( 'Safe minimum' )
00151 *
00152 *     Estimate the norm of inv(A).
00153 *
00154       AINVNM = ZERO
00155       NORMIN = 'N'
00156       IF( ONENRM ) THEN
00157          KASE1 = 1
00158       ELSE
00159          KASE1 = 2
00160       END IF
00161       KD = KL + KU + 1
00162       LNOTI = KL.GT.0
00163       KASE = 0
00164    10 CONTINUE
00165       CALL CLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
00166       IF( KASE.NE.0 ) THEN
00167          IF( KASE.EQ.KASE1 ) THEN
00168 *
00169 *           Multiply by inv(L).
00170 *
00171             IF( LNOTI ) THEN
00172                DO 20 J = 1, N - 1
00173                   LM = MIN( KL, N-J )
00174                   JP = IPIV( J )
00175                   T = WORK( JP )
00176                   IF( JP.NE.J ) THEN
00177                      WORK( JP ) = WORK( J )
00178                      WORK( J ) = T
00179                   END IF
00180                   CALL CAXPY( LM, -T, AB( KD+1, J ), 1, WORK( J+1 ), 1 )
00181    20          CONTINUE
00182             END IF
00183 *
00184 *           Multiply by inv(U).
00185 *
00186             CALL CLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
00187      $                   KL+KU, AB, LDAB, WORK, SCALE, RWORK, INFO )
00188          ELSE
00189 *
00190 *           Multiply by inv(U').
00191 *
00192             CALL CLATBS( 'Upper', 'Conjugate transpose', 'Non-unit',
00193      $                   NORMIN, N, KL+KU, AB, LDAB, WORK, SCALE, RWORK,
00194      $                   INFO )
00195 *
00196 *           Multiply by inv(L').
00197 *
00198             IF( LNOTI ) THEN
00199                DO 30 J = N - 1, 1, -1
00200                   LM = MIN( KL, N-J )
00201                   WORK( J ) = WORK( J ) - CDOTC( LM, AB( KD+1, J ), 1,
00202      $                        WORK( J+1 ), 1 )
00203                   JP = IPIV( J )
00204                   IF( JP.NE.J ) THEN
00205                      T = WORK( JP )
00206                      WORK( JP ) = WORK( J )
00207                      WORK( J ) = T
00208                   END IF
00209    30          CONTINUE
00210             END IF
00211          END IF
00212 *
00213 *        Divide X by 1/SCALE if doing so will not cause overflow.
00214 *
00215          NORMIN = 'Y'
00216          IF( SCALE.NE.ONE ) THEN
00217             IX = ICAMAX( N, WORK, 1 )
00218             IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
00219      $         GO TO 40
00220             CALL CSRSCL( N, SCALE, WORK, 1 )
00221          END IF
00222          GO TO 10
00223       END IF
00224 *
00225 *     Compute the estimate of the reciprocal condition number.
00226 *
00227       IF( AINVNM.NE.ZERO )
00228      $   RCOND = ( ONE / AINVNM ) / ANORM
00229 *
00230    40 CONTINUE
00231       RETURN
00232 *
00233 *     End of CGBCON
00234 *
00235       END