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