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