![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CPTCON 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CPTCON + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cptcon.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cptcon.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cptcon.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CPTCON( N, D, E, ANORM, RCOND, RWORK, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER INFO, N 00025 * REAL ANORM, RCOND 00026 * .. 00027 * .. Array Arguments .. 00028 * REAL D( * ), RWORK( * ) 00029 * COMPLEX E( * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> CPTCON computes the reciprocal of the condition number (in the 00039 *> 1-norm) of a complex Hermitian positive definite tridiagonal matrix 00040 *> using the factorization A = L*D*L**H or A = U**H*D*U computed by 00041 *> CPTTRF. 00042 *> 00043 *> Norm(inv(A)) is computed by a direct method, and the reciprocal of 00044 *> the condition number is computed as 00045 *> RCOND = 1 / (ANORM * norm(inv(A))). 00046 *> \endverbatim 00047 * 00048 * Arguments: 00049 * ========== 00050 * 00051 *> \param[in] N 00052 *> \verbatim 00053 *> N is INTEGER 00054 *> The order of the matrix A. N >= 0. 00055 *> \endverbatim 00056 *> 00057 *> \param[in] D 00058 *> \verbatim 00059 *> D is REAL array, dimension (N) 00060 *> The n diagonal elements of the diagonal matrix D from the 00061 *> factorization of A, as computed by CPTTRF. 00062 *> \endverbatim 00063 *> 00064 *> \param[in] E 00065 *> \verbatim 00066 *> E is COMPLEX array, dimension (N-1) 00067 *> The (n-1) off-diagonal elements of the unit bidiagonal factor 00068 *> U or L from the factorization of A, as computed by CPTTRF. 00069 *> \endverbatim 00070 *> 00071 *> \param[in] ANORM 00072 *> \verbatim 00073 *> ANORM is REAL 00074 *> The 1-norm of the original matrix A. 00075 *> \endverbatim 00076 *> 00077 *> \param[out] RCOND 00078 *> \verbatim 00079 *> RCOND is REAL 00080 *> The reciprocal of the condition number of the matrix A, 00081 *> computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is the 00082 *> 1-norm of inv(A) computed in this routine. 00083 *> \endverbatim 00084 *> 00085 *> \param[out] RWORK 00086 *> \verbatim 00087 *> RWORK is REAL array, dimension (N) 00088 *> \endverbatim 00089 *> 00090 *> \param[out] INFO 00091 *> \verbatim 00092 *> INFO is INTEGER 00093 *> = 0: successful exit 00094 *> < 0: if INFO = -i, the i-th argument had an illegal value 00095 *> \endverbatim 00096 * 00097 * Authors: 00098 * ======== 00099 * 00100 *> \author Univ. of Tennessee 00101 *> \author Univ. of California Berkeley 00102 *> \author Univ. of Colorado Denver 00103 *> \author NAG Ltd. 00104 * 00105 *> \date November 2011 00106 * 00107 *> \ingroup complexOTHERcomputational 00108 * 00109 *> \par Further Details: 00110 * ===================== 00111 *> 00112 *> \verbatim 00113 *> 00114 *> The method used is described in Nicholas J. Higham, "Efficient 00115 *> Algorithms for Computing the Condition Number of a Tridiagonal 00116 *> Matrix", SIAM J. Sci. Stat. Comput., Vol. 7, No. 1, January 1986. 00117 *> \endverbatim 00118 *> 00119 * ===================================================================== 00120 SUBROUTINE CPTCON( N, D, E, ANORM, RCOND, RWORK, INFO ) 00121 * 00122 * -- LAPACK computational routine (version 3.4.0) -- 00123 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00124 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00125 * November 2011 00126 * 00127 * .. Scalar Arguments .. 00128 INTEGER INFO, N 00129 REAL ANORM, RCOND 00130 * .. 00131 * .. Array Arguments .. 00132 REAL D( * ), RWORK( * ) 00133 COMPLEX E( * ) 00134 * .. 00135 * 00136 * ===================================================================== 00137 * 00138 * .. Parameters .. 00139 REAL ONE, ZERO 00140 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00141 * .. 00142 * .. Local Scalars .. 00143 INTEGER I, IX 00144 REAL AINVNM 00145 * .. 00146 * .. External Functions .. 00147 INTEGER ISAMAX 00148 EXTERNAL ISAMAX 00149 * .. 00150 * .. External Subroutines .. 00151 EXTERNAL XERBLA 00152 * .. 00153 * .. Intrinsic Functions .. 00154 INTRINSIC ABS 00155 * .. 00156 * .. Executable Statements .. 00157 * 00158 * Test the input arguments. 00159 * 00160 INFO = 0 00161 IF( N.LT.0 ) THEN 00162 INFO = -1 00163 ELSE IF( ANORM.LT.ZERO ) THEN 00164 INFO = -4 00165 END IF 00166 IF( INFO.NE.0 ) THEN 00167 CALL XERBLA( 'CPTCON', -INFO ) 00168 RETURN 00169 END IF 00170 * 00171 * Quick return if possible 00172 * 00173 RCOND = ZERO 00174 IF( N.EQ.0 ) THEN 00175 RCOND = ONE 00176 RETURN 00177 ELSE IF( ANORM.EQ.ZERO ) THEN 00178 RETURN 00179 END IF 00180 * 00181 * Check that D(1:N) is positive. 00182 * 00183 DO 10 I = 1, N 00184 IF( D( I ).LE.ZERO ) 00185 $ RETURN 00186 10 CONTINUE 00187 * 00188 * Solve M(A) * x = e, where M(A) = (m(i,j)) is given by 00189 * 00190 * m(i,j) = abs(A(i,j)), i = j, 00191 * m(i,j) = -abs(A(i,j)), i .ne. j, 00192 * 00193 * and e = [ 1, 1, ..., 1 ]**T. Note M(A) = M(L)*D*M(L)**H. 00194 * 00195 * Solve M(L) * x = e. 00196 * 00197 RWORK( 1 ) = ONE 00198 DO 20 I = 2, N 00199 RWORK( I ) = ONE + RWORK( I-1 )*ABS( E( I-1 ) ) 00200 20 CONTINUE 00201 * 00202 * Solve D * M(L)**H * x = b. 00203 * 00204 RWORK( N ) = RWORK( N ) / D( N ) 00205 DO 30 I = N - 1, 1, -1 00206 RWORK( I ) = RWORK( I ) / D( I ) + RWORK( I+1 )*ABS( E( I ) ) 00207 30 CONTINUE 00208 * 00209 * Compute AINVNM = max(x(i)), 1<=i<=n. 00210 * 00211 IX = ISAMAX( N, RWORK, 1 ) 00212 AINVNM = ABS( RWORK( IX ) ) 00213 * 00214 * Compute the reciprocal condition number. 00215 * 00216 IF( AINVNM.NE.ZERO ) 00217 $ RCOND = ( ONE / AINVNM ) / ANORM 00218 * 00219 RETURN 00220 * 00221 * End of CPTCON 00222 * 00223 END