 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, RWORK, 00002 $ INFO ) 00003 * 00004 * -- LAPACK routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 00008 * 00009 * Modified to call CLACN2 in place of CLACON, 10 Feb 03, SJH. 00010 * 00011 * .. Scalar Arguments .. 00012 CHARACTER UPLO 00013 INTEGER INFO, LDA, N 00014 REAL ANORM, RCOND 00015 * .. 00016 * .. Array Arguments .. 00017 REAL RWORK( * ) 00018 COMPLEX A( LDA, * ), WORK( * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * CPOCON estimates the reciprocal of the condition number (in the 00025 * 1-norm) of a complex Hermitian positive definite matrix using the 00026 * Cholesky factorization A = U**H*U or A = L*L**H computed by CPOTRF. 00027 * 00028 * An estimate is obtained for norm(inv(A)), and the reciprocal of the 00029 * condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))). 00030 * 00031 * Arguments 00032 * ========= 00033 * 00034 * UPLO (input) CHARACTER*1 00035 * = 'U': Upper triangle of A is stored; 00036 * = 'L': Lower triangle of A is stored. 00037 * 00038 * N (input) INTEGER 00039 * The order of the matrix A. N >= 0. 00040 * 00041 * A (input) COMPLEX array, dimension (LDA,N) 00042 * The triangular factor U or L from the Cholesky factorization 00043 * A = U**H*U or A = L*L**H, as computed by CPOTRF. 00044 * 00045 * LDA (input) INTEGER 00046 * The leading dimension of the array A. LDA >= max(1,N). 00047 * 00048 * ANORM (input) REAL 00049 * The 1-norm (or infinity-norm) of the Hermitian matrix A. 00050 * 00051 * RCOND (output) REAL 00052 * The reciprocal of the condition number of the matrix A, 00053 * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an 00054 * estimate of the 1-norm of inv(A) computed in this routine. 00055 * 00056 * WORK (workspace) COMPLEX array, dimension (2*N) 00057 * 00058 * RWORK (workspace) REAL array, dimension (N) 00059 * 00060 * INFO (output) INTEGER 00061 * = 0: successful exit 00062 * < 0: if INFO = -i, the i-th argument had an illegal value 00063 * 00064 * ===================================================================== 00065 * 00066 * .. Parameters .. 00067 REAL ONE, ZERO 00068 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00069 * .. 00070 * .. Local Scalars .. 00071 LOGICAL UPPER 00072 CHARACTER NORMIN 00073 INTEGER IX, KASE 00074 REAL AINVNM, SCALE, SCALEL, SCALEU, SMLNUM 00075 COMPLEX ZDUM 00076 * .. 00077 * .. Local Arrays .. 00078 INTEGER ISAVE( 3 ) 00079 * .. 00080 * .. External Functions .. 00081 LOGICAL LSAME 00082 INTEGER ICAMAX 00083 REAL SLAMCH 00084 EXTERNAL LSAME, ICAMAX, SLAMCH 00085 * .. 00086 * .. External Subroutines .. 00087 EXTERNAL CLACN2, CLATRS, CSRSCL, XERBLA 00088 * .. 00089 * .. Intrinsic Functions .. 00090 INTRINSIC ABS, AIMAG, MAX, REAL 00091 * .. 00092 * .. Statement Functions .. 00093 REAL CABS1 00094 * .. 00095 * .. Statement Function definitions .. 00096 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) ) 00097 * .. 00098 * .. Executable Statements .. 00099 * 00100 * Test the input parameters. 00101 * 00102 INFO = 0 00103 UPPER = LSAME( UPLO, 'U' ) 00104 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00105 INFO = -1 00106 ELSE IF( N.LT.0 ) THEN 00107 INFO = -2 00108 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00109 INFO = -4 00110 ELSE IF( ANORM.LT.ZERO ) THEN 00111 INFO = -5 00112 END IF 00113 IF( INFO.NE.0 ) THEN 00114 CALL XERBLA( 'CPOCON', -INFO ) 00115 RETURN 00116 END IF 00117 * 00118 * Quick return if possible 00119 * 00120 RCOND = ZERO 00121 IF( N.EQ.0 ) THEN 00122 RCOND = ONE 00123 RETURN 00124 ELSE IF( ANORM.EQ.ZERO ) THEN 00125 RETURN 00126 END IF 00127 * 00128 SMLNUM = SLAMCH( 'Safe minimum' ) 00129 * 00130 * Estimate the 1-norm of inv(A). 00131 * 00132 KASE = 0 00133 NORMIN = 'N' 00134 10 CONTINUE 00135 CALL CLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) 00136 IF( KASE.NE.0 ) THEN 00137 IF( UPPER ) THEN 00138 * 00139 * Multiply by inv(U**H). 00140 * 00141 CALL CLATRS( 'Upper', 'Conjugate transpose', 'Non-unit', 00142 $ NORMIN, N, A, LDA, WORK, SCALEL, RWORK, INFO ) 00143 NORMIN = 'Y' 00144 * 00145 * Multiply by inv(U). 00146 * 00147 CALL CLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, 00148 $ A, LDA, WORK, SCALEU, RWORK, INFO ) 00149 ELSE 00150 * 00151 * Multiply by inv(L). 00152 * 00153 CALL CLATRS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N, 00154 $ A, LDA, WORK, SCALEL, RWORK, INFO ) 00155 NORMIN = 'Y' 00156 * 00157 * Multiply by inv(L**H). 00158 * 00159 CALL CLATRS( 'Lower', 'Conjugate transpose', 'Non-unit', 00160 $ NORMIN, N, A, LDA, WORK, SCALEU, RWORK, INFO ) 00161 END IF 00162 * 00163 * Multiply by 1/SCALE if doing so will not cause overflow. 00164 * 00165 SCALE = SCALEL*SCALEU 00166 IF( SCALE.NE.ONE ) THEN 00167 IX = ICAMAX( N, WORK, 1 ) 00168 IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) 00169 $ GO TO 20 00170 CALL CSRSCL( N, SCALE, WORK, 1 ) 00171 END IF 00172 GO TO 10 00173 END IF 00174 * 00175 * Compute the estimate of the reciprocal condition number. 00176 * 00177 IF( AINVNM.NE.ZERO ) 00178 $ RCOND = ( ONE / AINVNM ) / ANORM 00179 * 00180 20 CONTINUE 00181 RETURN 00182 * 00183 * End of CPOCON 00184 * 00185 END
1.7.3 
