![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZPOCON 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZPOCON + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zpocon.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zpocon.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zpocon.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, RWORK, 00022 * INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER UPLO 00026 * INTEGER INFO, LDA, N 00027 * DOUBLE PRECISION ANORM, RCOND 00028 * .. 00029 * .. Array Arguments .. 00030 * DOUBLE PRECISION RWORK( * ) 00031 * COMPLEX*16 A( LDA, * ), WORK( * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> ZPOCON estimates the reciprocal of the condition number (in the 00041 *> 1-norm) of a complex Hermitian positive definite matrix using the 00042 *> Cholesky factorization A = U**H*U or A = L*L**H computed by ZPOTRF. 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 COMPLEX*16 array, dimension (LDA,N) 00067 *> The triangular factor U or L from the Cholesky factorization 00068 *> A = U**H*U or A = L*L**H, as computed by ZPOTRF. 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 DOUBLE PRECISION 00080 *> The 1-norm (or infinity-norm) of the Hermitian matrix A. 00081 *> \endverbatim 00082 *> 00083 *> \param[out] RCOND 00084 *> \verbatim 00085 *> RCOND is DOUBLE PRECISION 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 COMPLEX*16 array, dimension (2*N) 00094 *> \endverbatim 00095 *> 00096 *> \param[out] RWORK 00097 *> \verbatim 00098 *> RWORK is DOUBLE PRECISION 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 complex16POcomputational 00119 * 00120 * ===================================================================== 00121 SUBROUTINE ZPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, RWORK, 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 DOUBLE PRECISION ANORM, RCOND 00133 * .. 00134 * .. Array Arguments .. 00135 DOUBLE PRECISION RWORK( * ) 00136 COMPLEX*16 A( LDA, * ), WORK( * ) 00137 * .. 00138 * 00139 * ===================================================================== 00140 * 00141 * .. Parameters .. 00142 DOUBLE PRECISION ONE, ZERO 00143 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00144 * .. 00145 * .. Local Scalars .. 00146 LOGICAL UPPER 00147 CHARACTER NORMIN 00148 INTEGER IX, KASE 00149 DOUBLE PRECISION AINVNM, SCALE, SCALEL, SCALEU, SMLNUM 00150 COMPLEX*16 ZDUM 00151 * .. 00152 * .. Local Arrays .. 00153 INTEGER ISAVE( 3 ) 00154 * .. 00155 * .. External Functions .. 00156 LOGICAL LSAME 00157 INTEGER IZAMAX 00158 DOUBLE PRECISION DLAMCH 00159 EXTERNAL LSAME, IZAMAX, DLAMCH 00160 * .. 00161 * .. External Subroutines .. 00162 EXTERNAL XERBLA, ZDRSCL, ZLACN2, ZLATRS 00163 * .. 00164 * .. Intrinsic Functions .. 00165 INTRINSIC ABS, DBLE, DIMAG, MAX 00166 * .. 00167 * .. Statement Functions .. 00168 DOUBLE PRECISION CABS1 00169 * .. 00170 * .. Statement Function definitions .. 00171 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 00172 * .. 00173 * .. Executable Statements .. 00174 * 00175 * Test the input parameters. 00176 * 00177 INFO = 0 00178 UPPER = LSAME( UPLO, 'U' ) 00179 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00180 INFO = -1 00181 ELSE IF( N.LT.0 ) THEN 00182 INFO = -2 00183 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00184 INFO = -4 00185 ELSE IF( ANORM.LT.ZERO ) THEN 00186 INFO = -5 00187 END IF 00188 IF( INFO.NE.0 ) THEN 00189 CALL XERBLA( 'ZPOCON', -INFO ) 00190 RETURN 00191 END IF 00192 * 00193 * Quick return if possible 00194 * 00195 RCOND = ZERO 00196 IF( N.EQ.0 ) THEN 00197 RCOND = ONE 00198 RETURN 00199 ELSE IF( ANORM.EQ.ZERO ) THEN 00200 RETURN 00201 END IF 00202 * 00203 SMLNUM = DLAMCH( 'Safe minimum' ) 00204 * 00205 * Estimate the 1-norm of inv(A). 00206 * 00207 KASE = 0 00208 NORMIN = 'N' 00209 10 CONTINUE 00210 CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) 00211 IF( KASE.NE.0 ) THEN 00212 IF( UPPER ) THEN 00213 * 00214 * Multiply by inv(U**H). 00215 * 00216 CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit', 00217 $ NORMIN, N, A, LDA, WORK, SCALEL, RWORK, INFO ) 00218 NORMIN = 'Y' 00219 * 00220 * Multiply by inv(U). 00221 * 00222 CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, 00223 $ A, LDA, WORK, SCALEU, RWORK, INFO ) 00224 ELSE 00225 * 00226 * Multiply by inv(L). 00227 * 00228 CALL ZLATRS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N, 00229 $ A, LDA, WORK, SCALEL, RWORK, INFO ) 00230 NORMIN = 'Y' 00231 * 00232 * Multiply by inv(L**H). 00233 * 00234 CALL ZLATRS( 'Lower', 'Conjugate transpose', 'Non-unit', 00235 $ NORMIN, N, A, LDA, WORK, SCALEU, RWORK, INFO ) 00236 END IF 00237 * 00238 * Multiply by 1/SCALE if doing so will not cause overflow. 00239 * 00240 SCALE = SCALEL*SCALEU 00241 IF( SCALE.NE.ONE ) THEN 00242 IX = IZAMAX( N, WORK, 1 ) 00243 IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) 00244 $ GO TO 20 00245 CALL ZDRSCL( N, SCALE, WORK, 1 ) 00246 END IF 00247 GO TO 10 00248 END IF 00249 * 00250 * Compute the estimate of the reciprocal condition number. 00251 * 00252 IF( AINVNM.NE.ZERO ) 00253 $ RCOND = ( ONE / AINVNM ) / ANORM 00254 * 00255 20 CONTINUE 00256 RETURN 00257 * 00258 * End of ZPOCON 00259 * 00260 END