00001 SUBROUTINE ZHECON( UPLO, N, A, LDA, IPIV, ANORM, RCOND, WORK, 00002 $ INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH. 00010 * 00011 * .. Scalar Arguments .. 00012 CHARACTER UPLO 00013 INTEGER INFO, LDA, N 00014 DOUBLE PRECISION ANORM, RCOND 00015 * .. 00016 * .. Array Arguments .. 00017 INTEGER IPIV( * ) 00018 COMPLEX*16 A( LDA, * ), WORK( * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * ZHECON estimates the reciprocal of the condition number of a complex 00025 * Hermitian matrix A using the factorization A = U*D*U**H or 00026 * A = L*D*L**H computed by ZHETRF. 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 * Specifies whether the details of the factorization are stored 00036 * as an upper or lower triangular matrix. 00037 * = 'U': Upper triangular, form is A = U*D*U**H; 00038 * = 'L': Lower triangular, form is A = L*D*L**H. 00039 * 00040 * N (input) INTEGER 00041 * The order of the matrix A. N >= 0. 00042 * 00043 * A (input) COMPLEX*16 array, dimension (LDA,N) 00044 * The block diagonal matrix D and the multipliers used to 00045 * obtain the factor U or L as computed by ZHETRF. 00046 * 00047 * LDA (input) INTEGER 00048 * The leading dimension of the array A. LDA >= max(1,N). 00049 * 00050 * IPIV (input) INTEGER array, dimension (N) 00051 * Details of the interchanges and the block structure of D 00052 * as determined by ZHETRF. 00053 * 00054 * ANORM (input) DOUBLE PRECISION 00055 * The 1-norm of the original matrix A. 00056 * 00057 * RCOND (output) DOUBLE PRECISION 00058 * The reciprocal of the condition number of the matrix A, 00059 * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an 00060 * estimate of the 1-norm of inv(A) computed in this routine. 00061 * 00062 * WORK (workspace) COMPLEX*16 array, dimension (2*N) 00063 * 00064 * INFO (output) INTEGER 00065 * = 0: successful exit 00066 * < 0: if INFO = -i, the i-th argument had an illegal value 00067 * 00068 * ===================================================================== 00069 * 00070 * .. Parameters .. 00071 DOUBLE PRECISION ONE, ZERO 00072 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00073 * .. 00074 * .. Local Scalars .. 00075 LOGICAL UPPER 00076 INTEGER I, KASE 00077 DOUBLE PRECISION AINVNM 00078 * .. 00079 * .. Local Arrays .. 00080 INTEGER ISAVE( 3 ) 00081 * .. 00082 * .. External Functions .. 00083 LOGICAL LSAME 00084 EXTERNAL LSAME 00085 * .. 00086 * .. External Subroutines .. 00087 EXTERNAL XERBLA, ZHETRS, ZLACN2 00088 * .. 00089 * .. Intrinsic Functions .. 00090 INTRINSIC MAX 00091 * .. 00092 * .. Executable Statements .. 00093 * 00094 * Test the input parameters. 00095 * 00096 INFO = 0 00097 UPPER = LSAME( UPLO, 'U' ) 00098 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00099 INFO = -1 00100 ELSE IF( N.LT.0 ) THEN 00101 INFO = -2 00102 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00103 INFO = -4 00104 ELSE IF( ANORM.LT.ZERO ) THEN 00105 INFO = -6 00106 END IF 00107 IF( INFO.NE.0 ) THEN 00108 CALL XERBLA( 'ZHECON', -INFO ) 00109 RETURN 00110 END IF 00111 * 00112 * Quick return if possible 00113 * 00114 RCOND = ZERO 00115 IF( N.EQ.0 ) THEN 00116 RCOND = ONE 00117 RETURN 00118 ELSE IF( ANORM.LE.ZERO ) THEN 00119 RETURN 00120 END IF 00121 * 00122 * Check that the diagonal matrix D is nonsingular. 00123 * 00124 IF( UPPER ) THEN 00125 * 00126 * Upper triangular storage: examine D from bottom to top 00127 * 00128 DO 10 I = N, 1, -1 00129 IF( IPIV( I ).GT.0 .AND. A( I, I ).EQ.ZERO ) 00130 $ RETURN 00131 10 CONTINUE 00132 ELSE 00133 * 00134 * Lower triangular storage: examine D from top to bottom. 00135 * 00136 DO 20 I = 1, N 00137 IF( IPIV( I ).GT.0 .AND. A( I, I ).EQ.ZERO ) 00138 $ RETURN 00139 20 CONTINUE 00140 END IF 00141 * 00142 * Estimate the 1-norm of the inverse. 00143 * 00144 KASE = 0 00145 30 CONTINUE 00146 CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) 00147 IF( KASE.NE.0 ) THEN 00148 * 00149 * Multiply by inv(L*D*L') or inv(U*D*U'). 00150 * 00151 CALL ZHETRS( UPLO, N, 1, A, LDA, IPIV, WORK, N, INFO ) 00152 GO TO 30 00153 END IF 00154 * 00155 * Compute the estimate of the reciprocal condition number. 00156 * 00157 IF( AINVNM.NE.ZERO ) 00158 $ RCOND = ( ONE / AINVNM ) / ANORM 00159 * 00160 RETURN 00161 * 00162 * End of ZHECON 00163 * 00164 END
1.7.2