1 SUBROUTINE DPTCON( N, D, E, ANORM, RCOND, WORK, 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 * .. Scalar Arguments ..
9 INTEGER INFO, N
10 DOUBLE PRECISION ANORM, RCOND
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION D( * ), E( * ), WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DPTCON computes the reciprocal of the condition number (in the
20 * 1-norm) of a real symmetric positive definite tridiagonal matrix
21 * using the factorization A = L*D*L**T or A = U**T*D*U computed by
22 * DPTTRF.
23 *
24 * Norm(inv(A)) is computed by a direct method, and the reciprocal of
25 * the condition number is computed as
26 * RCOND = 1 / (ANORM * norm(inv(A))).
27 *
28 * Arguments
29 * =========
30 *
31 * N (input) INTEGER
32 * The order of the matrix A. N >= 0.
33 *
34 * D (input) DOUBLE PRECISION array, dimension (N)
35 * The n diagonal elements of the diagonal matrix D from the
36 * factorization of A, as computed by DPTTRF.
37 *
38 * E (input) DOUBLE PRECISION array, dimension (N-1)
39 * The (n-1) off-diagonal elements of the unit bidiagonal factor
40 * U or L from the factorization of A, as computed by DPTTRF.
41 *
42 * ANORM (input) DOUBLE PRECISION
43 * The 1-norm of the original matrix A.
44 *
45 * RCOND (output) DOUBLE PRECISION
46 * The reciprocal of the condition number of the matrix A,
47 * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is the
48 * 1-norm of inv(A) computed in this routine.
49 *
50 * WORK (workspace) DOUBLE PRECISION array, dimension (N)
51 *
52 * INFO (output) INTEGER
53 * = 0: successful exit
54 * < 0: if INFO = -i, the i-th argument had an illegal value
55 *
56 * Further Details
57 * ===============
58 *
59 * The method used is described in Nicholas J. Higham, "Efficient
60 * Algorithms for Computing the Condition Number of a Tridiagonal
61 * Matrix", SIAM J. Sci. Stat. Comput., Vol. 7, No. 1, January 1986.
62 *
63 * =====================================================================
64 *
65 * .. Parameters ..
66 DOUBLE PRECISION ONE, ZERO
67 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
68 * ..
69 * .. Local Scalars ..
70 INTEGER I, IX
71 DOUBLE PRECISION AINVNM
72 * ..
73 * .. External Functions ..
74 INTEGER IDAMAX
75 EXTERNAL IDAMAX
76 * ..
77 * .. External Subroutines ..
78 EXTERNAL XERBLA
79 * ..
80 * .. Intrinsic Functions ..
81 INTRINSIC ABS
82 * ..
83 * .. Executable Statements ..
84 *
85 * Test the input arguments.
86 *
87 INFO = 0
88 IF( N.LT.0 ) THEN
89 INFO = -1
90 ELSE IF( ANORM.LT.ZERO ) THEN
91 INFO = -4
92 END IF
93 IF( INFO.NE.0 ) THEN
94 CALL XERBLA( 'DPTCON', -INFO )
95 RETURN
96 END IF
97 *
98 * Quick return if possible
99 *
100 RCOND = ZERO
101 IF( N.EQ.0 ) THEN
102 RCOND = ONE
103 RETURN
104 ELSE IF( ANORM.EQ.ZERO ) THEN
105 RETURN
106 END IF
107 *
108 * Check that D(1:N) is positive.
109 *
110 DO 10 I = 1, N
111 IF( D( I ).LE.ZERO )
112 $ RETURN
113 10 CONTINUE
114 *
115 * Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
116 *
117 * m(i,j) = abs(A(i,j)), i = j,
118 * m(i,j) = -abs(A(i,j)), i .ne. j,
119 *
120 * and e = [ 1, 1, ..., 1 ]**T. Note M(A) = M(L)*D*M(L)**T.
121 *
122 * Solve M(L) * x = e.
123 *
124 WORK( 1 ) = ONE
125 DO 20 I = 2, N
126 WORK( I ) = ONE + WORK( I-1 )*ABS( E( I-1 ) )
127 20 CONTINUE
128 *
129 * Solve D * M(L)**T * x = b.
130 *
131 WORK( N ) = WORK( N ) / D( N )
132 DO 30 I = N - 1, 1, -1
133 WORK( I ) = WORK( I ) / D( I ) + WORK( I+1 )*ABS( E( I ) )
134 30 CONTINUE
135 *
136 * Compute AINVNM = max(x(i)), 1<=i<=n.
137 *
138 IX = IDAMAX( N, WORK, 1 )
139 AINVNM = ABS( WORK( IX ) )
140 *
141 * Compute the reciprocal condition number.
142 *
143 IF( AINVNM.NE.ZERO )
144 $ RCOND = ( ONE / AINVNM ) / ANORM
145 *
146 RETURN
147 *
148 * End of DPTCON
149 *
150 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 * .. Scalar Arguments ..
9 INTEGER INFO, N
10 DOUBLE PRECISION ANORM, RCOND
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION D( * ), E( * ), WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DPTCON computes the reciprocal of the condition number (in the
20 * 1-norm) of a real symmetric positive definite tridiagonal matrix
21 * using the factorization A = L*D*L**T or A = U**T*D*U computed by
22 * DPTTRF.
23 *
24 * Norm(inv(A)) is computed by a direct method, and the reciprocal of
25 * the condition number is computed as
26 * RCOND = 1 / (ANORM * norm(inv(A))).
27 *
28 * Arguments
29 * =========
30 *
31 * N (input) INTEGER
32 * The order of the matrix A. N >= 0.
33 *
34 * D (input) DOUBLE PRECISION array, dimension (N)
35 * The n diagonal elements of the diagonal matrix D from the
36 * factorization of A, as computed by DPTTRF.
37 *
38 * E (input) DOUBLE PRECISION array, dimension (N-1)
39 * The (n-1) off-diagonal elements of the unit bidiagonal factor
40 * U or L from the factorization of A, as computed by DPTTRF.
41 *
42 * ANORM (input) DOUBLE PRECISION
43 * The 1-norm of the original matrix A.
44 *
45 * RCOND (output) DOUBLE PRECISION
46 * The reciprocal of the condition number of the matrix A,
47 * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is the
48 * 1-norm of inv(A) computed in this routine.
49 *
50 * WORK (workspace) DOUBLE PRECISION array, dimension (N)
51 *
52 * INFO (output) INTEGER
53 * = 0: successful exit
54 * < 0: if INFO = -i, the i-th argument had an illegal value
55 *
56 * Further Details
57 * ===============
58 *
59 * The method used is described in Nicholas J. Higham, "Efficient
60 * Algorithms for Computing the Condition Number of a Tridiagonal
61 * Matrix", SIAM J. Sci. Stat. Comput., Vol. 7, No. 1, January 1986.
62 *
63 * =====================================================================
64 *
65 * .. Parameters ..
66 DOUBLE PRECISION ONE, ZERO
67 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
68 * ..
69 * .. Local Scalars ..
70 INTEGER I, IX
71 DOUBLE PRECISION AINVNM
72 * ..
73 * .. External Functions ..
74 INTEGER IDAMAX
75 EXTERNAL IDAMAX
76 * ..
77 * .. External Subroutines ..
78 EXTERNAL XERBLA
79 * ..
80 * .. Intrinsic Functions ..
81 INTRINSIC ABS
82 * ..
83 * .. Executable Statements ..
84 *
85 * Test the input arguments.
86 *
87 INFO = 0
88 IF( N.LT.0 ) THEN
89 INFO = -1
90 ELSE IF( ANORM.LT.ZERO ) THEN
91 INFO = -4
92 END IF
93 IF( INFO.NE.0 ) THEN
94 CALL XERBLA( 'DPTCON', -INFO )
95 RETURN
96 END IF
97 *
98 * Quick return if possible
99 *
100 RCOND = ZERO
101 IF( N.EQ.0 ) THEN
102 RCOND = ONE
103 RETURN
104 ELSE IF( ANORM.EQ.ZERO ) THEN
105 RETURN
106 END IF
107 *
108 * Check that D(1:N) is positive.
109 *
110 DO 10 I = 1, N
111 IF( D( I ).LE.ZERO )
112 $ RETURN
113 10 CONTINUE
114 *
115 * Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
116 *
117 * m(i,j) = abs(A(i,j)), i = j,
118 * m(i,j) = -abs(A(i,j)), i .ne. j,
119 *
120 * and e = [ 1, 1, ..., 1 ]**T. Note M(A) = M(L)*D*M(L)**T.
121 *
122 * Solve M(L) * x = e.
123 *
124 WORK( 1 ) = ONE
125 DO 20 I = 2, N
126 WORK( I ) = ONE + WORK( I-1 )*ABS( E( I-1 ) )
127 20 CONTINUE
128 *
129 * Solve D * M(L)**T * x = b.
130 *
131 WORK( N ) = WORK( N ) / D( N )
132 DO 30 I = N - 1, 1, -1
133 WORK( I ) = WORK( I ) / D( I ) + WORK( I+1 )*ABS( E( I ) )
134 30 CONTINUE
135 *
136 * Compute AINVNM = max(x(i)), 1<=i<=n.
137 *
138 IX = IDAMAX( N, WORK, 1 )
139 AINVNM = ABS( WORK( IX ) )
140 *
141 * Compute the reciprocal condition number.
142 *
143 IF( AINVNM.NE.ZERO )
144 $ RCOND = ( ONE / AINVNM ) / ANORM
145 *
146 RETURN
147 *
148 * End of DPTCON
149 *
150 END