![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SPOCON 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SPOCON + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/spocon.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/spocon.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/spocon.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK, 00022 * INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER UPLO 00026 * INTEGER INFO, LDA, N 00027 * REAL ANORM, RCOND 00028 * .. 00029 * .. Array Arguments .. 00030 * INTEGER IWORK( * ) 00031 * REAL A( LDA, * ), WORK( * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> SPOCON estimates the reciprocal of the condition number (in the 00041 *> 1-norm) of a real symmetric positive definite matrix using the 00042 *> Cholesky factorization A = U**T*U or A = L*L**T computed by SPOTRF. 00043 *> 00044 *> An estimate is obtained for norm(inv(A)), and the reciprocal of the 00045 *> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))). 00046 *> \endverbatim 00047 * 00048 * Arguments: 00049 * ========== 00050 * 00051 *> \param[in] UPLO 00052 *> \verbatim 00053 *> UPLO is CHARACTER*1 00054 *> = 'U': Upper triangle of A is stored; 00055 *> = 'L': Lower triangle of A is stored. 00056 *> \endverbatim 00057 *> 00058 *> \param[in] N 00059 *> \verbatim 00060 *> N is INTEGER 00061 *> The order of the matrix A. N >= 0. 00062 *> \endverbatim 00063 *> 00064 *> \param[in] A 00065 *> \verbatim 00066 *> A is REAL array, dimension (LDA,N) 00067 *> The triangular factor U or L from the Cholesky factorization 00068 *> A = U**T*U or A = L*L**T, as computed by SPOTRF. 00069 *> \endverbatim 00070 *> 00071 *> \param[in] LDA 00072 *> \verbatim 00073 *> LDA is INTEGER 00074 *> The leading dimension of the array A. LDA >= max(1,N). 00075 *> \endverbatim 00076 *> 00077 *> \param[in] ANORM 00078 *> \verbatim 00079 *> ANORM is REAL 00080 *> The 1-norm (or infinity-norm) of the symmetric matrix A. 00081 *> \endverbatim 00082 *> 00083 *> \param[out] RCOND 00084 *> \verbatim 00085 *> RCOND is REAL 00086 *> The reciprocal of the condition number of the matrix A, 00087 *> computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an 00088 *> estimate of the 1-norm of inv(A) computed in this routine. 00089 *> \endverbatim 00090 *> 00091 *> \param[out] WORK 00092 *> \verbatim 00093 *> WORK is REAL array, dimension (3*N) 00094 *> \endverbatim 00095 *> 00096 *> \param[out] IWORK 00097 *> \verbatim 00098 *> IWORK is INTEGER array, dimension (N) 00099 *> \endverbatim 00100 *> 00101 *> \param[out] INFO 00102 *> \verbatim 00103 *> INFO is INTEGER 00104 *> = 0: successful exit 00105 *> < 0: if INFO = -i, the i-th argument had an illegal value 00106 *> \endverbatim 00107 * 00108 * Authors: 00109 * ======== 00110 * 00111 *> \author Univ. of Tennessee 00112 *> \author Univ. of California Berkeley 00113 *> \author Univ. of Colorado Denver 00114 *> \author NAG Ltd. 00115 * 00116 *> \date November 2011 00117 * 00118 *> \ingroup realPOcomputational 00119 * 00120 * ===================================================================== 00121 SUBROUTINE SPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK, 00122 $ INFO ) 00123 * 00124 * -- LAPACK computational routine (version 3.4.0) -- 00125 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00126 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00127 * November 2011 00128 * 00129 * .. Scalar Arguments .. 00130 CHARACTER UPLO 00131 INTEGER INFO, LDA, N 00132 REAL ANORM, RCOND 00133 * .. 00134 * .. Array Arguments .. 00135 INTEGER IWORK( * ) 00136 REAL A( LDA, * ), WORK( * ) 00137 * .. 00138 * 00139 * ===================================================================== 00140 * 00141 * .. Parameters .. 00142 REAL ONE, ZERO 00143 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00144 * .. 00145 * .. Local Scalars .. 00146 LOGICAL UPPER 00147 CHARACTER NORMIN 00148 INTEGER IX, KASE 00149 REAL AINVNM, SCALE, SCALEL, SCALEU, SMLNUM 00150 * .. 00151 * .. Local Arrays .. 00152 INTEGER ISAVE( 3 ) 00153 * .. 00154 * .. External Functions .. 00155 LOGICAL LSAME 00156 INTEGER ISAMAX 00157 REAL SLAMCH 00158 EXTERNAL LSAME, ISAMAX, SLAMCH 00159 * .. 00160 * .. External Subroutines .. 00161 EXTERNAL SLACN2, SLATRS, SRSCL, XERBLA 00162 * .. 00163 * .. Intrinsic Functions .. 00164 INTRINSIC ABS, MAX 00165 * .. 00166 * .. Executable Statements .. 00167 * 00168 * Test the input parameters. 00169 * 00170 INFO = 0 00171 UPPER = LSAME( UPLO, 'U' ) 00172 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00173 INFO = -1 00174 ELSE IF( N.LT.0 ) THEN 00175 INFO = -2 00176 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00177 INFO = -4 00178 ELSE IF( ANORM.LT.ZERO ) THEN 00179 INFO = -5 00180 END IF 00181 IF( INFO.NE.0 ) THEN 00182 CALL XERBLA( 'SPOCON', -INFO ) 00183 RETURN 00184 END IF 00185 * 00186 * Quick return if possible 00187 * 00188 RCOND = ZERO 00189 IF( N.EQ.0 ) THEN 00190 RCOND = ONE 00191 RETURN 00192 ELSE IF( ANORM.EQ.ZERO ) THEN 00193 RETURN 00194 END IF 00195 * 00196 SMLNUM = SLAMCH( 'Safe minimum' ) 00197 * 00198 * Estimate the 1-norm of inv(A). 00199 * 00200 KASE = 0 00201 NORMIN = 'N' 00202 10 CONTINUE 00203 CALL SLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) 00204 IF( KASE.NE.0 ) THEN 00205 IF( UPPER ) THEN 00206 * 00207 * Multiply by inv(U**T). 00208 * 00209 CALL SLATRS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, A, 00210 $ LDA, WORK, SCALEL, WORK( 2*N+1 ), INFO ) 00211 NORMIN = 'Y' 00212 * 00213 * Multiply by inv(U). 00214 * 00215 CALL SLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, 00216 $ A, LDA, WORK, SCALEU, WORK( 2*N+1 ), INFO ) 00217 ELSE 00218 * 00219 * Multiply by inv(L). 00220 * 00221 CALL SLATRS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N, 00222 $ A, LDA, WORK, SCALEL, WORK( 2*N+1 ), INFO ) 00223 NORMIN = 'Y' 00224 * 00225 * Multiply by inv(L**T). 00226 * 00227 CALL SLATRS( 'Lower', 'Transpose', 'Non-unit', NORMIN, N, A, 00228 $ LDA, WORK, SCALEU, WORK( 2*N+1 ), INFO ) 00229 END IF 00230 * 00231 * Multiply by 1/SCALE if doing so will not cause overflow. 00232 * 00233 SCALE = SCALEL*SCALEU 00234 IF( SCALE.NE.ONE ) THEN 00235 IX = ISAMAX( N, WORK, 1 ) 00236 IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) 00237 $ GO TO 20 00238 CALL SRSCL( N, SCALE, WORK, 1 ) 00239 END IF 00240 GO TO 10 00241 END IF 00242 * 00243 * Compute the estimate of the reciprocal condition number. 00244 * 00245 IF( AINVNM.NE.ZERO ) 00246 $ RCOND = ( ONE / AINVNM ) / ANORM 00247 * 00248 20 CONTINUE 00249 RETURN 00250 * 00251 * End of SPOCON 00252 * 00253 END