Annotation of rpl/lapack/lapack/dgbcon.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DGBCON
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DGBCON + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgbcon.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgbcon.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgbcon.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
! 22: * WORK, IWORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * CHARACTER NORM
! 26: * INTEGER INFO, KL, KU, LDAB, N
! 27: * DOUBLE PRECISION ANORM, RCOND
! 28: * ..
! 29: * .. Array Arguments ..
! 30: * INTEGER IPIV( * ), IWORK( * )
! 31: * DOUBLE PRECISION AB( LDAB, * ), WORK( * )
! 32: * ..
! 33: *
! 34: *
! 35: *> \par Purpose:
! 36: * =============
! 37: *>
! 38: *> \verbatim
! 39: *>
! 40: *> DGBCON estimates the reciprocal of the condition number of a real
! 41: *> general band matrix A, in either the 1-norm or the infinity-norm,
! 42: *> using the LU factorization computed by DGBTRF.
! 43: *>
! 44: *> An estimate is obtained for norm(inv(A)), and the reciprocal of the
! 45: *> condition number is computed as
! 46: *> RCOND = 1 / ( norm(A) * norm(inv(A)) ).
! 47: *> \endverbatim
! 48: *
! 49: * Arguments:
! 50: * ==========
! 51: *
! 52: *> \param[in] NORM
! 53: *> \verbatim
! 54: *> NORM is CHARACTER*1
! 55: *> Specifies whether the 1-norm condition number or the
! 56: *> infinity-norm condition number is required:
! 57: *> = '1' or 'O': 1-norm;
! 58: *> = 'I': Infinity-norm.
! 59: *> \endverbatim
! 60: *>
! 61: *> \param[in] N
! 62: *> \verbatim
! 63: *> N is INTEGER
! 64: *> The order of the matrix A. N >= 0.
! 65: *> \endverbatim
! 66: *>
! 67: *> \param[in] KL
! 68: *> \verbatim
! 69: *> KL is INTEGER
! 70: *> The number of subdiagonals within the band of A. KL >= 0.
! 71: *> \endverbatim
! 72: *>
! 73: *> \param[in] KU
! 74: *> \verbatim
! 75: *> KU is INTEGER
! 76: *> The number of superdiagonals within the band of A. KU >= 0.
! 77: *> \endverbatim
! 78: *>
! 79: *> \param[in] AB
! 80: *> \verbatim
! 81: *> AB is DOUBLE PRECISION array, dimension (LDAB,N)
! 82: *> Details of the LU factorization of the band matrix A, as
! 83: *> computed by DGBTRF. U is stored as an upper triangular band
! 84: *> matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
! 85: *> the multipliers used during the factorization are stored in
! 86: *> rows KL+KU+2 to 2*KL+KU+1.
! 87: *> \endverbatim
! 88: *>
! 89: *> \param[in] LDAB
! 90: *> \verbatim
! 91: *> LDAB is INTEGER
! 92: *> The leading dimension of the array AB. LDAB >= 2*KL+KU+1.
! 93: *> \endverbatim
! 94: *>
! 95: *> \param[in] IPIV
! 96: *> \verbatim
! 97: *> IPIV is INTEGER array, dimension (N)
! 98: *> The pivot indices; for 1 <= i <= N, row i of the matrix was
! 99: *> interchanged with row IPIV(i).
! 100: *> \endverbatim
! 101: *>
! 102: *> \param[in] ANORM
! 103: *> \verbatim
! 104: *> ANORM is DOUBLE PRECISION
! 105: *> If NORM = '1' or 'O', the 1-norm of the original matrix A.
! 106: *> If NORM = 'I', the infinity-norm of the original matrix A.
! 107: *> \endverbatim
! 108: *>
! 109: *> \param[out] RCOND
! 110: *> \verbatim
! 111: *> RCOND is DOUBLE PRECISION
! 112: *> The reciprocal of the condition number of the matrix A,
! 113: *> computed as RCOND = 1/(norm(A) * norm(inv(A))).
! 114: *> \endverbatim
! 115: *>
! 116: *> \param[out] WORK
! 117: *> \verbatim
! 118: *> WORK is DOUBLE PRECISION array, dimension (3*N)
! 119: *> \endverbatim
! 120: *>
! 121: *> \param[out] IWORK
! 122: *> \verbatim
! 123: *> IWORK is INTEGER array, dimension (N)
! 124: *> \endverbatim
! 125: *>
! 126: *> \param[out] INFO
! 127: *> \verbatim
! 128: *> INFO is INTEGER
! 129: *> = 0: successful exit
! 130: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 131: *> \endverbatim
! 132: *
! 133: * Authors:
! 134: * ========
! 135: *
! 136: *> \author Univ. of Tennessee
! 137: *> \author Univ. of California Berkeley
! 138: *> \author Univ. of Colorado Denver
! 139: *> \author NAG Ltd.
! 140: *
! 141: *> \date November 2011
! 142: *
! 143: *> \ingroup doubleGBcomputational
! 144: *
! 145: * =====================================================================
1.1 bertrand 146: SUBROUTINE DGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND,
147: $ WORK, IWORK, INFO )
148: *
1.9 ! bertrand 149: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 150: * -- LAPACK is a software package provided by Univ. of Tennessee, --
151: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 152: * November 2011
1.1 bertrand 153: *
154: * .. Scalar Arguments ..
155: CHARACTER NORM
156: INTEGER INFO, KL, KU, LDAB, N
157: DOUBLE PRECISION ANORM, RCOND
158: * ..
159: * .. Array Arguments ..
160: INTEGER IPIV( * ), IWORK( * )
161: DOUBLE PRECISION AB( LDAB, * ), WORK( * )
162: * ..
163: *
164: * =====================================================================
165: *
166: * .. Parameters ..
167: DOUBLE PRECISION ONE, ZERO
168: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
169: * ..
170: * .. Local Scalars ..
171: LOGICAL LNOTI, ONENRM
172: CHARACTER NORMIN
173: INTEGER IX, J, JP, KASE, KASE1, KD, LM
174: DOUBLE PRECISION AINVNM, SCALE, SMLNUM, T
175: * ..
176: * .. Local Arrays ..
177: INTEGER ISAVE( 3 )
178: * ..
179: * .. External Functions ..
180: LOGICAL LSAME
181: INTEGER IDAMAX
182: DOUBLE PRECISION DDOT, DLAMCH
183: EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH
184: * ..
185: * .. External Subroutines ..
186: EXTERNAL DAXPY, DLACN2, DLATBS, DRSCL, XERBLA
187: * ..
188: * .. Intrinsic Functions ..
189: INTRINSIC ABS, MIN
190: * ..
191: * .. Executable Statements ..
192: *
193: * Test the input parameters.
194: *
195: INFO = 0
196: ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
197: IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
198: INFO = -1
199: ELSE IF( N.LT.0 ) THEN
200: INFO = -2
201: ELSE IF( KL.LT.0 ) THEN
202: INFO = -3
203: ELSE IF( KU.LT.0 ) THEN
204: INFO = -4
205: ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN
206: INFO = -6
207: ELSE IF( ANORM.LT.ZERO ) THEN
208: INFO = -8
209: END IF
210: IF( INFO.NE.0 ) THEN
211: CALL XERBLA( 'DGBCON', -INFO )
212: RETURN
213: END IF
214: *
215: * Quick return if possible
216: *
217: RCOND = ZERO
218: IF( N.EQ.0 ) THEN
219: RCOND = ONE
220: RETURN
221: ELSE IF( ANORM.EQ.ZERO ) THEN
222: RETURN
223: END IF
224: *
225: SMLNUM = DLAMCH( 'Safe minimum' )
226: *
227: * Estimate the norm of inv(A).
228: *
229: AINVNM = ZERO
230: NORMIN = 'N'
231: IF( ONENRM ) THEN
232: KASE1 = 1
233: ELSE
234: KASE1 = 2
235: END IF
236: KD = KL + KU + 1
237: LNOTI = KL.GT.0
238: KASE = 0
239: 10 CONTINUE
240: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
241: IF( KASE.NE.0 ) THEN
242: IF( KASE.EQ.KASE1 ) THEN
243: *
244: * Multiply by inv(L).
245: *
246: IF( LNOTI ) THEN
247: DO 20 J = 1, N - 1
248: LM = MIN( KL, N-J )
249: JP = IPIV( J )
250: T = WORK( JP )
251: IF( JP.NE.J ) THEN
252: WORK( JP ) = WORK( J )
253: WORK( J ) = T
254: END IF
255: CALL DAXPY( LM, -T, AB( KD+1, J ), 1, WORK( J+1 ), 1 )
256: 20 CONTINUE
257: END IF
258: *
259: * Multiply by inv(U).
260: *
261: CALL DLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
262: $ KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
263: $ INFO )
264: ELSE
265: *
1.8 bertrand 266: * Multiply by inv(U**T).
1.1 bertrand 267: *
268: CALL DLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N,
269: $ KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ),
270: $ INFO )
271: *
1.8 bertrand 272: * Multiply by inv(L**T).
1.1 bertrand 273: *
274: IF( LNOTI ) THEN
275: DO 30 J = N - 1, 1, -1
276: LM = MIN( KL, N-J )
277: WORK( J ) = WORK( J ) - DDOT( LM, AB( KD+1, J ), 1,
278: $ WORK( J+1 ), 1 )
279: JP = IPIV( J )
280: IF( JP.NE.J ) THEN
281: T = WORK( JP )
282: WORK( JP ) = WORK( J )
283: WORK( J ) = T
284: END IF
285: 30 CONTINUE
286: END IF
287: END IF
288: *
289: * Divide X by 1/SCALE if doing so will not cause overflow.
290: *
291: NORMIN = 'Y'
292: IF( SCALE.NE.ONE ) THEN
293: IX = IDAMAX( N, WORK, 1 )
294: IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
295: $ GO TO 40
296: CALL DRSCL( N, SCALE, WORK, 1 )
297: END IF
298: GO TO 10
299: END IF
300: *
301: * Compute the estimate of the reciprocal condition number.
302: *
303: IF( AINVNM.NE.ZERO )
304: $ RCOND = ( ONE / AINVNM ) / ANORM
305: *
306: 40 CONTINUE
307: RETURN
308: *
309: * End of DGBCON
310: *
311: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgbcon.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack