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