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