3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download DTPCON + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtpcon.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtpcon.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtpcon.f">
21 * SUBROUTINE DTPCON( NORM, UPLO, DIAG, N, AP, RCOND, WORK, IWORK,
24 * .. Scalar Arguments ..
25 * CHARACTER DIAG, NORM, UPLO
27 * DOUBLE PRECISION RCOND
29 * .. Array Arguments ..
31 * DOUBLE PRECISION AP( * ), WORK( * )
40 *> DTPCON estimates the reciprocal of the condition number of a packed
41 *> triangular matrix A, in either the 1-norm or the infinity-norm.
43 *> The norm of A is computed and an estimate is obtained for
44 *> norm(inv(A)), then the reciprocal of the condition number is
46 *> RCOND = 1 / ( norm(A) * norm(inv(A)) ).
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.
63 *> UPLO is CHARACTER*1
64 *> = 'U': A is upper triangular;
65 *> = 'L': A is lower triangular.
70 *> DIAG is CHARACTER*1
71 *> = 'N': A is non-unit triangular;
72 *> = 'U': A is unit triangular.
78 *> The order of the matrix A. N >= 0.
83 *> AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
84 *> The upper or lower triangular matrix A, packed columnwise in
85 *> a linear array. The j-th column of A is stored in the array
87 *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
88 *> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
89 *> If DIAG = 'U', the diagonal elements of A are not referenced
90 *> and are assumed to be 1.
95 *> RCOND is DOUBLE PRECISION
96 *> The reciprocal of the condition number of the matrix A,
97 *> computed as RCOND = 1/(norm(A) * norm(inv(A))).
102 *> WORK is DOUBLE PRECISION array, dimension (3*N)
107 *> IWORK is INTEGER array, dimension (N)
113 *> = 0: successful exit
114 *> < 0: if INFO = -i, the i-th argument had an illegal value
120 *> \author Univ. of Tennessee
121 *> \author Univ. of California Berkeley
122 *> \author Univ. of Colorado Denver
125 *> \date November 2011
127 *> \ingroup doubleOTHERcomputational
129 * =====================================================================
130 SUBROUTINE DTPCON( NORM, UPLO, DIAG, N, AP, RCOND, WORK, IWORK,
133 * -- LAPACK computational routine (version 3.4.0) --
134 * -- LAPACK is a software package provided by Univ. of Tennessee, --
135 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138 * .. Scalar Arguments ..
139 CHARACTER DIAG, NORM, UPLO
141 DOUBLE PRECISION RCOND
143 * .. Array Arguments ..
145 DOUBLE PRECISION AP( * ), WORK( * )
148 * =====================================================================
151 DOUBLE PRECISION ONE, ZERO
152 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
154 * .. Local Scalars ..
155 LOGICAL NOUNIT, ONENRM, UPPER
157 INTEGER IX, KASE, KASE1
158 DOUBLE PRECISION AINVNM, ANORM, SCALE, SMLNUM, XNORM
163 * .. External Functions ..
166 DOUBLE PRECISION DLAMCH, DLANTP
167 EXTERNAL LSAME, IDAMAX, DLAMCH, DLANTP
169 * .. External Subroutines ..
170 EXTERNAL DLACN2, DLATPS, DRSCL, XERBLA
172 * .. Intrinsic Functions ..
173 INTRINSIC ABS, DBLE, MAX
175 * .. Executable Statements ..
177 * Test the input parameters.
180 UPPER = LSAME( UPLO, 'U' )
181 ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
182 NOUNIT = LSAME( DIAG, 'N' )
184 IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
186 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
188 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
190 ELSE IF( N.LT.0 ) THEN
194 CALL XERBLA( 'DTPCON', -INFO )
198 * Quick return if possible
206 SMLNUM = DLAMCH( 'Safe minimum' )*DBLE( MAX( 1, N ) )
208 * Compute the norm of the triangular matrix A.
210 ANORM = DLANTP( NORM, UPLO, DIAG, N, AP, WORK )
212 * Continue only if ANORM > 0.
214 IF( ANORM.GT.ZERO ) THEN
216 * Estimate the norm of the inverse of A.
227 CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
229 IF( KASE.EQ.KASE1 ) THEN
231 * Multiply by inv(A).
233 CALL DLATPS( UPLO, 'No transpose', DIAG, NORMIN, N, AP,
234 $ WORK, SCALE, WORK( 2*N+1 ), INFO )
237 * Multiply by inv(A**T).
239 CALL DLATPS( UPLO, 'Transpose', DIAG, NORMIN, N, AP,
240 $ WORK, SCALE, WORK( 2*N+1 ), INFO )
244 * Multiply by 1/SCALE if doing so will not cause overflow.
246 IF( SCALE.NE.ONE ) THEN
247 IX = IDAMAX( N, WORK, 1 )
248 XNORM = ABS( WORK( IX ) )
249 IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
251 CALL DRSCL( N, SCALE, WORK, 1 )
256 * Compute the estimate of the reciprocal condition number.
259 $ RCOND = ( ONE / ANORM ) / AINVNM