1*> \brief \b CTPCON
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CTPCON + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctpcon.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctpcon.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctpcon.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CTPCON( NORM, UPLO, DIAG, N, AP, RCOND, WORK, RWORK,
22*                          INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          DIAG, NORM, UPLO
26*       INTEGER            INFO, N
27*       REAL               RCOND
28*       ..
29*       .. Array Arguments ..
30*       REAL               RWORK( * )
31*       COMPLEX            AP( * ), WORK( * )
32*       ..
33*
34*
35*> \par Purpose:
36*  =============
37*>
38*> \verbatim
39*>
40*> CTPCON estimates the reciprocal of the condition number of a packed
41*> triangular matrix A, in either the 1-norm or the infinity-norm.
42*>
43*> The norm of A is computed and an estimate is obtained for
44*> norm(inv(A)), then the reciprocal of the condition number is
45*> computed as
46*>    RCOND = 1 / ( norm(A) * norm(inv(A)) ).
47*> \endverbatim
48*
49*  Arguments:
50*  ==========
51*
52*> \param[in] NORM
53*> \verbatim
54*>          NORM is CHARACTER*1
55*>          Specifies whether the 1-norm condition number or the
56*>          infinity-norm condition number is required:
57*>          = '1' or 'O':  1-norm;
58*>          = 'I':         Infinity-norm.
59*> \endverbatim
60*>
61*> \param[in] UPLO
62*> \verbatim
63*>          UPLO is CHARACTER*1
64*>          = 'U':  A is upper triangular;
65*>          = 'L':  A is lower triangular.
66*> \endverbatim
67*>
68*> \param[in] DIAG
69*> \verbatim
70*>          DIAG is CHARACTER*1
71*>          = 'N':  A is non-unit triangular;
72*>          = 'U':  A is unit triangular.
73*> \endverbatim
74*>
75*> \param[in] N
76*> \verbatim
77*>          N is INTEGER
78*>          The order of the matrix A.  N >= 0.
79*> \endverbatim
80*>
81*> \param[in] AP
82*> \verbatim
83*>          AP is COMPLEX array, dimension (N*(N+1)/2)
84*>          The upper or lower triangular matrix A, packed columnwise in
85*>          a linear array.  The j-th column of A is stored in the array
86*>          AP as follows:
87*>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
88*>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
89*>          If DIAG = 'U', the diagonal elements of A are not referenced
90*>          and are assumed to be 1.
91*> \endverbatim
92*>
93*> \param[out] RCOND
94*> \verbatim
95*>          RCOND is REAL
96*>          The reciprocal of the condition number of the matrix A,
97*>          computed as RCOND = 1/(norm(A) * norm(inv(A))).
98*> \endverbatim
99*>
100*> \param[out] WORK
101*> \verbatim
102*>          WORK is COMPLEX array, dimension (2*N)
103*> \endverbatim
104*>
105*> \param[out] RWORK
106*> \verbatim
107*>          RWORK is REAL array, dimension (N)
108*> \endverbatim
109*>
110*> \param[out] INFO
111*> \verbatim
112*>          INFO is INTEGER
113*>          = 0:  successful exit
114*>          < 0:  if INFO = -i, the i-th argument had an illegal value
115*> \endverbatim
116*
117*  Authors:
118*  ========
119*
120*> \author Univ. of Tennessee
121*> \author Univ. of California Berkeley
122*> \author Univ. of Colorado Denver
123*> \author NAG Ltd.
124*
125*> \ingroup complexOTHERcomputational
126*
127*  =====================================================================
128      SUBROUTINE CTPCON( NORM, UPLO, DIAG, N, AP, RCOND, WORK, RWORK,
129     $                   INFO )
130*
131*  -- LAPACK computational routine --
132*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
133*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
134*
135*     .. Scalar Arguments ..
136      CHARACTER          DIAG, NORM, UPLO
137      INTEGER            INFO, N
138      REAL               RCOND
139*     ..
140*     .. Array Arguments ..
141      REAL               RWORK( * )
142      COMPLEX            AP( * ), WORK( * )
143*     ..
144*
145*  =====================================================================
146*
147*     .. Parameters ..
148      REAL               ONE, ZERO
149      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
150*     ..
151*     .. Local Scalars ..
152      LOGICAL            NOUNIT, ONENRM, UPPER
153      CHARACTER          NORMIN
154      INTEGER            IX, KASE, KASE1
155      REAL               AINVNM, ANORM, SCALE, SMLNUM, XNORM
156      COMPLEX            ZDUM
157*     ..
158*     .. Local Arrays ..
159      INTEGER            ISAVE( 3 )
160*     ..
161*     .. External Functions ..
162      LOGICAL            LSAME
163      INTEGER            ICAMAX
164      REAL               CLANTP, SLAMCH
165      EXTERNAL           LSAME, ICAMAX, CLANTP, SLAMCH
166*     ..
167*     .. External Subroutines ..
168      EXTERNAL           CLACN2, CLATPS, CSRSCL, XERBLA
169*     ..
170*     .. Intrinsic Functions ..
171      INTRINSIC          ABS, AIMAG, MAX, REAL
172*     ..
173*     .. Statement Functions ..
174      REAL               CABS1
175*     ..
176*     .. Statement Function definitions ..
177      CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
178*     ..
179*     .. Executable Statements ..
180*
181*     Test the input parameters.
182*
183      INFO = 0
184      UPPER = LSAME( UPLO, 'U' )
185      ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
186      NOUNIT = LSAME( DIAG, 'N' )
187*
188      IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
189         INFO = -1
190      ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
191         INFO = -2
192      ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
193         INFO = -3
194      ELSE IF( N.LT.0 ) THEN
195         INFO = -4
196      END IF
197      IF( INFO.NE.0 ) THEN
198         CALL XERBLA( 'CTPCON', -INFO )
199         RETURN
200      END IF
201*
202*     Quick return if possible
203*
204      IF( N.EQ.0 ) THEN
205         RCOND = ONE
206         RETURN
207      END IF
208*
209      RCOND = ZERO
210      SMLNUM = SLAMCH( 'Safe minimum' )*REAL( MAX( 1, N ) )
211*
212*     Compute the norm of the triangular matrix A.
213*
214      ANORM = CLANTP( NORM, UPLO, DIAG, N, AP, RWORK )
215*
216*     Continue only if ANORM > 0.
217*
218      IF( ANORM.GT.ZERO ) THEN
219*
220*        Estimate the norm of the inverse of A.
221*
222         AINVNM = ZERO
223         NORMIN = 'N'
224         IF( ONENRM ) THEN
225            KASE1 = 1
226         ELSE
227            KASE1 = 2
228         END IF
229         KASE = 0
230   10    CONTINUE
231         CALL CLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
232         IF( KASE.NE.0 ) THEN
233            IF( KASE.EQ.KASE1 ) THEN
234*
235*              Multiply by inv(A).
236*
237               CALL CLATPS( UPLO, 'No transpose', DIAG, NORMIN, N, AP,
238     $                      WORK, SCALE, RWORK, INFO )
239            ELSE
240*
241*              Multiply by inv(A**H).
242*
243               CALL CLATPS( UPLO, 'Conjugate transpose', DIAG, NORMIN,
244     $                      N, AP, WORK, SCALE, RWORK, INFO )
245            END IF
246            NORMIN = 'Y'
247*
248*           Multiply by 1/SCALE if doing so will not cause overflow.
249*
250            IF( SCALE.NE.ONE ) THEN
251               IX = ICAMAX( N, WORK, 1 )
252               XNORM = CABS1( WORK( IX ) )
253               IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
254     $            GO TO 20
255               CALL CSRSCL( N, SCALE, WORK, 1 )
256            END IF
257            GO TO 10
258         END IF
259*
260*        Compute the estimate of the reciprocal condition number.
261*
262         IF( AINVNM.NE.ZERO )
263     $      RCOND = ( ONE / ANORM ) / AINVNM
264      END IF
265*
266   20 CONTINUE
267      RETURN
268*
269*     End of CTPCON
270*
271      END
272