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*> \date November 2011
106*
107*> \ingroup complex16_lin
108*
109*  =====================================================================
110      SUBROUTINE ZTPT01( 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      DOUBLE PRECISION   RCOND, RESID
121*     ..
122*     .. Array Arguments ..
123      DOUBLE PRECISION   RWORK( * )
124      COMPLEX*16         AINVP( * ), AP( * )
125*     ..
126*
127*  =====================================================================
128*
129*     .. Parameters ..
130      DOUBLE PRECISION   ZERO, ONE
131      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
132*     ..
133*     .. Local Scalars ..
134      LOGICAL            UNITD
135      INTEGER            J, JC
136      DOUBLE PRECISION   AINVNM, ANORM, EPS
137*     ..
138*     .. External Functions ..
139      LOGICAL            LSAME
140      DOUBLE PRECISION   DLAMCH, ZLANTP
141      EXTERNAL           LSAME, DLAMCH, ZLANTP
142*     ..
143*     .. External Subroutines ..
144      EXTERNAL           ZTPMV
145*     ..
146*     .. Intrinsic Functions ..
147      INTRINSIC          DBLE
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 = DLAMCH( 'Epsilon' )
162      ANORM = ZLANTP( '1', UPLO, DIAG, N, AP, RWORK )
163      AINVNM = ZLANTP( '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 ZTPMV( '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 ZTPMV( '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 = ZLANTP( '1', UPLO, 'Non-unit', N, AINVP, RWORK )
211*
212      RESID = ( ( RESID*RCOND ) / DBLE( N ) ) / EPS
213*
214      RETURN
215*
216*     End of ZTPT01
217*
218      END
219