1 SUBROUTINE DPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK,
2 $ INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
10 *
11 * .. Scalar Arguments ..
12 CHARACTER UPLO
13 INTEGER INFO, LDA, N
14 DOUBLE PRECISION ANORM, RCOND
15 * ..
16 * .. Array Arguments ..
17 INTEGER IWORK( * )
18 DOUBLE PRECISION A( LDA, * ), WORK( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DPOCON estimates the reciprocal of the condition number (in the
25 * 1-norm) of a real symmetric positive definite matrix using the
26 * Cholesky factorization A = U**T*U or A = L*L**T computed by DPOTRF.
27 *
28 * An estimate is obtained for norm(inv(A)), and the reciprocal of the
29 * condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
30 *
31 * Arguments
32 * =========
33 *
34 * UPLO (input) CHARACTER*1
35 * = 'U': Upper triangle of A is stored;
36 * = 'L': Lower triangle of A is stored.
37 *
38 * N (input) INTEGER
39 * The order of the matrix A. N >= 0.
40 *
41 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
42 * The triangular factor U or L from the Cholesky factorization
43 * A = U**T*U or A = L*L**T, as computed by DPOTRF.
44 *
45 * LDA (input) INTEGER
46 * The leading dimension of the array A. LDA >= max(1,N).
47 *
48 * ANORM (input) DOUBLE PRECISION
49 * The 1-norm (or infinity-norm) of the symmetric matrix A.
50 *
51 * RCOND (output) DOUBLE PRECISION
52 * The reciprocal of the condition number of the matrix A,
53 * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
54 * estimate of the 1-norm of inv(A) computed in this routine.
55 *
56 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
57 *
58 * IWORK (workspace) INTEGER array, dimension (N)
59 *
60 * INFO (output) INTEGER
61 * = 0: successful exit
62 * < 0: if INFO = -i, the i-th argument had an illegal value
63 *
64 * =====================================================================
65 *
66 * .. Parameters ..
67 DOUBLE PRECISION ONE, ZERO
68 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
69 * ..
70 * .. Local Scalars ..
71 LOGICAL UPPER
72 CHARACTER NORMIN
73 INTEGER IX, KASE
74 DOUBLE PRECISION AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
75 * ..
76 * .. Local Arrays ..
77 INTEGER ISAVE( 3 )
78 * ..
79 * .. External Functions ..
80 LOGICAL LSAME
81 INTEGER IDAMAX
82 DOUBLE PRECISION DLAMCH
83 EXTERNAL LSAME, IDAMAX, DLAMCH
84 * ..
85 * .. External Subroutines ..
86 EXTERNAL DLACN2, DLATRS, DRSCL, XERBLA
87 * ..
88 * .. Intrinsic Functions ..
89 INTRINSIC ABS, MAX
90 * ..
91 * .. Executable Statements ..
92 *
93 * Test the input parameters.
94 *
95 INFO = 0
96 UPPER = LSAME( UPLO, 'U' )
97 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
98 INFO = -1
99 ELSE IF( N.LT.0 ) THEN
100 INFO = -2
101 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
102 INFO = -4
103 ELSE IF( ANORM.LT.ZERO ) THEN
104 INFO = -5
105 END IF
106 IF( INFO.NE.0 ) THEN
107 CALL XERBLA( 'DPOCON', -INFO )
108 RETURN
109 END IF
110 *
111 * Quick return if possible
112 *
113 RCOND = ZERO
114 IF( N.EQ.0 ) THEN
115 RCOND = ONE
116 RETURN
117 ELSE IF( ANORM.EQ.ZERO ) THEN
118 RETURN
119 END IF
120 *
121 SMLNUM = DLAMCH( 'Safe minimum' )
122 *
123 * Estimate the 1-norm of inv(A).
124 *
125 KASE = 0
126 NORMIN = 'N'
127 10 CONTINUE
128 CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
129 IF( KASE.NE.0 ) THEN
130 IF( UPPER ) THEN
131 *
132 * Multiply by inv(U**T).
133 *
134 CALL DLATRS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, A,
135 $ LDA, WORK, SCALEL, WORK( 2*N+1 ), INFO )
136 NORMIN = 'Y'
137 *
138 * Multiply by inv(U).
139 *
140 CALL DLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
141 $ A, LDA, WORK, SCALEU, WORK( 2*N+1 ), INFO )
142 ELSE
143 *
144 * Multiply by inv(L).
145 *
146 CALL DLATRS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
147 $ A, LDA, WORK, SCALEL, WORK( 2*N+1 ), INFO )
148 NORMIN = 'Y'
149 *
150 * Multiply by inv(L**T).
151 *
152 CALL DLATRS( 'Lower', 'Transpose', 'Non-unit', NORMIN, N, A,
153 $ LDA, WORK, SCALEU, WORK( 2*N+1 ), INFO )
154 END IF
155 *
156 * Multiply by 1/SCALE if doing so will not cause overflow.
157 *
158 SCALE = SCALEL*SCALEU
159 IF( SCALE.NE.ONE ) THEN
160 IX = IDAMAX( N, WORK, 1 )
161 IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
162 $ GO TO 20
163 CALL DRSCL( N, SCALE, WORK, 1 )
164 END IF
165 GO TO 10
166 END IF
167 *
168 * Compute the estimate of the reciprocal condition number.
169 *
170 IF( AINVNM.NE.ZERO )
171 $ RCOND = ( ONE / AINVNM ) / ANORM
172 *
173 20 CONTINUE
174 RETURN
175 *
176 * End of DPOCON
177 *
178 END
2 $ INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
10 *
11 * .. Scalar Arguments ..
12 CHARACTER UPLO
13 INTEGER INFO, LDA, N
14 DOUBLE PRECISION ANORM, RCOND
15 * ..
16 * .. Array Arguments ..
17 INTEGER IWORK( * )
18 DOUBLE PRECISION A( LDA, * ), WORK( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DPOCON estimates the reciprocal of the condition number (in the
25 * 1-norm) of a real symmetric positive definite matrix using the
26 * Cholesky factorization A = U**T*U or A = L*L**T computed by DPOTRF.
27 *
28 * An estimate is obtained for norm(inv(A)), and the reciprocal of the
29 * condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
30 *
31 * Arguments
32 * =========
33 *
34 * UPLO (input) CHARACTER*1
35 * = 'U': Upper triangle of A is stored;
36 * = 'L': Lower triangle of A is stored.
37 *
38 * N (input) INTEGER
39 * The order of the matrix A. N >= 0.
40 *
41 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
42 * The triangular factor U or L from the Cholesky factorization
43 * A = U**T*U or A = L*L**T, as computed by DPOTRF.
44 *
45 * LDA (input) INTEGER
46 * The leading dimension of the array A. LDA >= max(1,N).
47 *
48 * ANORM (input) DOUBLE PRECISION
49 * The 1-norm (or infinity-norm) of the symmetric matrix A.
50 *
51 * RCOND (output) DOUBLE PRECISION
52 * The reciprocal of the condition number of the matrix A,
53 * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
54 * estimate of the 1-norm of inv(A) computed in this routine.
55 *
56 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
57 *
58 * IWORK (workspace) INTEGER array, dimension (N)
59 *
60 * INFO (output) INTEGER
61 * = 0: successful exit
62 * < 0: if INFO = -i, the i-th argument had an illegal value
63 *
64 * =====================================================================
65 *
66 * .. Parameters ..
67 DOUBLE PRECISION ONE, ZERO
68 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
69 * ..
70 * .. Local Scalars ..
71 LOGICAL UPPER
72 CHARACTER NORMIN
73 INTEGER IX, KASE
74 DOUBLE PRECISION AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
75 * ..
76 * .. Local Arrays ..
77 INTEGER ISAVE( 3 )
78 * ..
79 * .. External Functions ..
80 LOGICAL LSAME
81 INTEGER IDAMAX
82 DOUBLE PRECISION DLAMCH
83 EXTERNAL LSAME, IDAMAX, DLAMCH
84 * ..
85 * .. External Subroutines ..
86 EXTERNAL DLACN2, DLATRS, DRSCL, XERBLA
87 * ..
88 * .. Intrinsic Functions ..
89 INTRINSIC ABS, MAX
90 * ..
91 * .. Executable Statements ..
92 *
93 * Test the input parameters.
94 *
95 INFO = 0
96 UPPER = LSAME( UPLO, 'U' )
97 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
98 INFO = -1
99 ELSE IF( N.LT.0 ) THEN
100 INFO = -2
101 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
102 INFO = -4
103 ELSE IF( ANORM.LT.ZERO ) THEN
104 INFO = -5
105 END IF
106 IF( INFO.NE.0 ) THEN
107 CALL XERBLA( 'DPOCON', -INFO )
108 RETURN
109 END IF
110 *
111 * Quick return if possible
112 *
113 RCOND = ZERO
114 IF( N.EQ.0 ) THEN
115 RCOND = ONE
116 RETURN
117 ELSE IF( ANORM.EQ.ZERO ) THEN
118 RETURN
119 END IF
120 *
121 SMLNUM = DLAMCH( 'Safe minimum' )
122 *
123 * Estimate the 1-norm of inv(A).
124 *
125 KASE = 0
126 NORMIN = 'N'
127 10 CONTINUE
128 CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
129 IF( KASE.NE.0 ) THEN
130 IF( UPPER ) THEN
131 *
132 * Multiply by inv(U**T).
133 *
134 CALL DLATRS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, A,
135 $ LDA, WORK, SCALEL, WORK( 2*N+1 ), INFO )
136 NORMIN = 'Y'
137 *
138 * Multiply by inv(U).
139 *
140 CALL DLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
141 $ A, LDA, WORK, SCALEU, WORK( 2*N+1 ), INFO )
142 ELSE
143 *
144 * Multiply by inv(L).
145 *
146 CALL DLATRS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
147 $ A, LDA, WORK, SCALEL, WORK( 2*N+1 ), INFO )
148 NORMIN = 'Y'
149 *
150 * Multiply by inv(L**T).
151 *
152 CALL DLATRS( 'Lower', 'Transpose', 'Non-unit', NORMIN, N, A,
153 $ LDA, WORK, SCALEU, WORK( 2*N+1 ), INFO )
154 END IF
155 *
156 * Multiply by 1/SCALE if doing so will not cause overflow.
157 *
158 SCALE = SCALEL*SCALEU
159 IF( SCALE.NE.ONE ) THEN
160 IX = IDAMAX( N, WORK, 1 )
161 IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
162 $ GO TO 20
163 CALL DRSCL( N, SCALE, WORK, 1 )
164 END IF
165 GO TO 10
166 END IF
167 *
168 * Compute the estimate of the reciprocal condition number.
169 *
170 IF( AINVNM.NE.ZERO )
171 $ RCOND = ( ONE / AINVNM ) / ANORM
172 *
173 20 CONTINUE
174 RETURN
175 *
176 * End of DPOCON
177 *
178 END