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*> \ingroup complex_lin
106*
107*  =====================================================================
108      SUBROUTINE CTPT01( 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      REAL               RCOND, RESID
118*     ..
119*     .. Array Arguments ..
120      REAL               RWORK( * )
121      COMPLEX            AINVP( * ), AP( * )
122*     ..
123*
124*  =====================================================================
125*
126*     .. Parameters ..
127      REAL               ZERO, ONE
128      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
129*     ..
130*     .. Local Scalars ..
131      LOGICAL            UNITD
132      INTEGER            J, JC
133      REAL               AINVNM, ANORM, EPS
134*     ..
135*     .. External Functions ..
136      LOGICAL            LSAME
137      REAL               CLANTP, SLAMCH
138      EXTERNAL           LSAME, CLANTP, SLAMCH
139*     ..
140*     .. External Subroutines ..
141      EXTERNAL           CTPMV
142*     ..
143*     .. Intrinsic Functions ..
144      INTRINSIC          REAL
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 = SLAMCH( 'Epsilon' )
159      ANORM = CLANTP( '1', UPLO, DIAG, N, AP, RWORK )
160      AINVNM = CLANTP( '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 CTPMV( '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 CTPMV( '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 = CLANTP( '1', UPLO, 'Non-unit', N, AINVP, RWORK )
208*
209      RESID = ( ( RESID*RCOND ) / REAL( N ) ) / EPS
210*
211      RETURN
212*
213*     End of CTPT01
214*
215      END
216