LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
ztpt01.f
Go to the documentation of this file.
1 *> \brief \b ZTPT01
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE ZTPT01( UPLO, DIAG, N, AP, AINVP, RCOND, RWORK, RESID )
12 *
13 * .. Scalar Arguments ..
14 * CHARACTER DIAG, UPLO
15 * INTEGER N
16 * DOUBLE PRECISION RCOND, RESID
17 * ..
18 * .. Array Arguments ..
19 * DOUBLE PRECISION RWORK( * )
20 * COMPLEX*16 AINVP( * ), AP( * )
21 * ..
22 *
23 *
24 *> \par Purpose:
25 * =============
26 *>
27 *> \verbatim
28 *>
29 *> ZTPT01 computes the residual for a triangular matrix A times its
30 *> inverse when A is stored in packed format:
31 *> RESID = norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS ),
32 *> where EPS is the machine epsilon.
33 *> \endverbatim
34 *
35 * Arguments:
36 * ==========
37 *
38 *> \param[in] UPLO
39 *> \verbatim
40 *> UPLO is CHARACTER*1
41 *> Specifies whether the matrix A is upper or lower triangular.
42 *> = 'U': Upper triangular
43 *> = 'L': Lower triangular
44 *> \endverbatim
45 *>
46 *> \param[in] DIAG
47 *> \verbatim
48 *> DIAG is CHARACTER*1
49 *> Specifies whether or not the matrix A is unit triangular.
50 *> = 'N': Non-unit triangular
51 *> = 'U': Unit triangular
52 *> \endverbatim
53 *>
54 *> \param[in] N
55 *> \verbatim
56 *> N is INTEGER
57 *> The order of the matrix A. N >= 0.
58 *> \endverbatim
59 *>
60 *> \param[in] AP
61 *> \verbatim
62 *> AP is COMPLEX*16 array, dimension (N*(N+1)/2)
63 *> The original upper or lower triangular matrix A, packed
64 *> columnwise in a linear array. The j-th column of A is stored
65 *> in the array AP as follows:
66 *> if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j;
67 *> if UPLO = 'L',
68 *> AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n.
69 *> \endverbatim
70 *>
71 *> \param[in] AINVP
72 *> \verbatim
73 *> AINVP is COMPLEX*16 array, dimension (N*(N+1)/2)
74 *> On entry, the (triangular) inverse of the matrix A, packed
75 *> columnwise in a linear array as in AP.
76 *> On exit, the contents of AINVP are destroyed.
77 *> \endverbatim
78 *>
79 *> \param[out] RCOND
80 *> \verbatim
81 *> RCOND is DOUBLE PRECISION
82 *> The reciprocal condition number of A, computed as
83 *> 1/(norm(A) * norm(AINV)).
84 *> \endverbatim
85 *>
86 *> \param[out] RWORK
87 *> \verbatim
88 *> RWORK is DOUBLE PRECISION array, dimension (N)
89 *> \endverbatim
90 *>
91 *> \param[out] RESID
92 *> \verbatim
93 *> RESID is DOUBLE PRECISION
94 *> norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS )
95 *> \endverbatim
96 *
97 * Authors:
98 * ========
99 *
100 *> \author Univ. of Tennessee
101 *> \author Univ. of California Berkeley
102 *> \author Univ. of Colorado Denver
103 *> \author NAG Ltd.
104 *
105 *> \ingroup complex16_lin
106 *
107 * =====================================================================
108  SUBROUTINE ztpt01( UPLO, DIAG, N, AP, AINVP, RCOND, RWORK, RESID )
109 *
110 * -- LAPACK test routine --
111 * -- LAPACK is a software package provided by Univ. of Tennessee, --
112 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
113 *
114 * .. Scalar Arguments ..
115  CHARACTER DIAG, UPLO
116  INTEGER N
117  DOUBLE PRECISION RCOND, RESID
118 * ..
119 * .. Array Arguments ..
120  DOUBLE PRECISION RWORK( * )
121  COMPLEX*16 AINVP( * ), AP( * )
122 * ..
123 *
124 * =====================================================================
125 *
126 * .. Parameters ..
127  DOUBLE PRECISION ZERO, ONE
128  parameter( zero = 0.0d+0, one = 1.0d+0 )
129 * ..
130 * .. Local Scalars ..
131  LOGICAL UNITD
132  INTEGER J, JC
133  DOUBLE PRECISION AINVNM, ANORM, EPS
134 * ..
135 * .. External Functions ..
136  LOGICAL LSAME
137  DOUBLE PRECISION DLAMCH, ZLANTP
138  EXTERNAL lsame, dlamch, zlantp
139 * ..
140 * .. External Subroutines ..
141  EXTERNAL ztpmv
142 * ..
143 * .. Intrinsic Functions ..
144  INTRINSIC dble
145 * ..
146 * .. Executable Statements ..
147 *
148 * Quick exit if N = 0.
149 *
150  IF( n.LE.0 ) THEN
151  rcond = one
152  resid = zero
153  RETURN
154  END IF
155 *
156 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
157 *
158  eps = dlamch( 'Epsilon' )
159  anorm = zlantp( '1', uplo, diag, n, ap, rwork )
160  ainvnm = zlantp( '1', uplo, diag, n, ainvp, rwork )
161  IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
162  rcond = zero
163  resid = one / eps
164  RETURN
165  END IF
166  rcond = ( one / anorm ) / ainvnm
167 *
168 * Compute A * AINV, overwriting AINV.
169 *
170  unitd = lsame( diag, 'U' )
171  IF( lsame( uplo, 'U' ) ) THEN
172  jc = 1
173  DO 10 j = 1, n
174  IF( unitd )
175  $ ainvp( jc+j-1 ) = one
176 *
177 * Form the j-th column of A*AINV.
178 *
179  CALL ztpmv( 'Upper', 'No transpose', diag, j, ap,
180  $ ainvp( jc ), 1 )
181 *
182 * Subtract 1 from the diagonal to form A*AINV - I.
183 *
184  ainvp( jc+j-1 ) = ainvp( jc+j-1 ) - one
185  jc = jc + j
186  10 CONTINUE
187  ELSE
188  jc = 1
189  DO 20 j = 1, n
190  IF( unitd )
191  $ ainvp( jc ) = one
192 *
193 * Form the j-th column of A*AINV.
194 *
195  CALL ztpmv( 'Lower', 'No transpose', diag, n-j+1, ap( jc ),
196  $ ainvp( jc ), 1 )
197 *
198 * Subtract 1 from the diagonal to form A*AINV - I.
199 *
200  ainvp( jc ) = ainvp( jc ) - one
201  jc = jc + n - j + 1
202  20 CONTINUE
203  END IF
204 *
205 * Compute norm(A*AINV - I) / (N * norm(A) * norm(AINV) * EPS)
206 *
207  resid = zlantp( '1', uplo, 'Non-unit', n, ainvp, rwork )
208 *
209  resid = ( ( resid*rcond ) / dble( n ) ) / eps
210 *
211  RETURN
212 *
213 * End of ZTPT01
214 *
215  END
subroutine ztpmv(UPLO, TRANS, DIAG, N, AP, X, INCX)
ZTPMV
Definition: ztpmv.f:142
subroutine ztpt01(UPLO, DIAG, N, AP, AINVP, RCOND, RWORK, RESID)
ZTPT01
Definition: ztpt01.f:109