LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
ctpt01.f
Go to the documentation of this file.
1 *> \brief \b CTPT01
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 CTPT01( UPLO, DIAG, N, AP, AINVP, RCOND, RWORK, RESID )
12 *
13 * .. Scalar Arguments ..
14 * CHARACTER DIAG, UPLO
15 * INTEGER N
16 * REAL RCOND, RESID
17 * ..
18 * .. Array Arguments ..
19 * REAL RWORK( * )
20 * COMPLEX AINVP( * ), AP( * )
21 * ..
22 *
23 *
24 *> \par Purpose:
25 * =============
26 *>
27 *> \verbatim
28 *>
29 *> CTPT01 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 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 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 REAL
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 REAL array, dimension (N)
89 *> \endverbatim
90 *>
91 *> \param[out] RESID
92 *> \verbatim
93 *> RESID is REAL
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 *> \date November 2011
106 *
107 *> \ingroup complex_lin
108 *
109 * =====================================================================
110  SUBROUTINE ctpt01( UPLO, DIAG, N, AP, AINVP, RCOND, RWORK, RESID )
111 *
112 * -- LAPACK test routine (version 3.4.0) --
113 * -- LAPACK is a software package provided by Univ. of Tennessee, --
114 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
115 * November 2011
116 *
117 * .. Scalar Arguments ..
118  CHARACTER diag, uplo
119  INTEGER n
120  REAL rcond, resid
121 * ..
122 * .. Array Arguments ..
123  REAL rwork( * )
124  COMPLEX ainvp( * ), ap( * )
125 * ..
126 *
127 * =====================================================================
128 *
129 * .. Parameters ..
130  REAL zero, one
131  parameter( zero = 0.0e+0, one = 1.0e+0 )
132 * ..
133 * .. Local Scalars ..
134  LOGICAL unitd
135  INTEGER j, jc
136  REAL ainvnm, anorm, eps
137 * ..
138 * .. External Functions ..
139  LOGICAL lsame
140  REAL clantp, slamch
141  EXTERNAL lsame, clantp, slamch
142 * ..
143 * .. External Subroutines ..
144  EXTERNAL ctpmv
145 * ..
146 * .. Intrinsic Functions ..
147  INTRINSIC real
148 * ..
149 * .. Executable Statements ..
150 *
151 * Quick exit if N = 0.
152 *
153  IF( n.LE.0 ) THEN
154  rcond = one
155  resid = zero
156  return
157  END IF
158 *
159 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
160 *
161  eps = slamch( 'Epsilon' )
162  anorm = clantp( '1', uplo, diag, n, ap, rwork )
163  ainvnm = clantp( '1', uplo, diag, n, ainvp, rwork )
164  IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
165  rcond = zero
166  resid = one / eps
167  return
168  END IF
169  rcond = ( one / anorm ) / ainvnm
170 *
171 * Compute A * AINV, overwriting AINV.
172 *
173  unitd = lsame( diag, 'U' )
174  IF( lsame( uplo, 'U' ) ) THEN
175  jc = 1
176  DO 10 j = 1, n
177  IF( unitd )
178  \$ ainvp( jc+j-1 ) = one
179 *
180 * Form the j-th column of A*AINV.
181 *
182  CALL ctpmv( 'Upper', 'No transpose', diag, j, ap,
183  \$ ainvp( jc ), 1 )
184 *
185 * Subtract 1 from the diagonal to form A*AINV - I.
186 *
187  ainvp( jc+j-1 ) = ainvp( jc+j-1 ) - one
188  jc = jc + j
189  10 continue
190  ELSE
191  jc = 1
192  DO 20 j = 1, n
193  IF( unitd )
194  \$ ainvp( jc ) = one
195 *
196 * Form the j-th column of A*AINV.
197 *
198  CALL ctpmv( 'Lower', 'No transpose', diag, n-j+1, ap( jc ),
199  \$ ainvp( jc ), 1 )
200 *
201 * Subtract 1 from the diagonal to form A*AINV - I.
202 *
203  ainvp( jc ) = ainvp( jc ) - one
204  jc = jc + n - j + 1
205  20 continue
206  END IF
207 *
208 * Compute norm(A*AINV - I) / (N * norm(A) * norm(AINV) * EPS)
209 *
210  resid = clantp( '1', uplo, 'Non-unit', n, ainvp, rwork )
211 *
212  resid = ( ( resid*rcond ) / REAL( N ) ) / eps
213 *
214  return
215 *
216 * End of CTPT01
217 *
218  END