1*> \brief \b STPCON
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download STPCON + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stpcon.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stpcon.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stpcon.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE STPCON( NORM, UPLO, DIAG, N, AP, RCOND, WORK, IWORK,
22*                          INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          DIAG, NORM, UPLO
26*       INTEGER            INFO, N
27*       REAL               RCOND
28*       ..
29*       .. Array Arguments ..
30*       INTEGER            IWORK( * )
31*       REAL               AP( * ), WORK( * )
32*       ..
33*
34*
35*> \par Purpose:
36*  =============
37*>
38*> \verbatim
39*>
40*> STPCON 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 REAL 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 REAL array, dimension (3*N)
103*> \endverbatim
104*>
105*> \param[out] IWORK
106*> \verbatim
107*>          IWORK is INTEGER 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*> \date November 2011
126*
127*> \ingroup realOTHERcomputational
128*
129*  =====================================================================
130      SUBROUTINE STPCON( NORM, UPLO, DIAG, N, AP, RCOND, WORK, IWORK,
131     $                   INFO )
132*
133*  -- LAPACK computational routine (version 3.4.0) --
134*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
135*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
136*     November 2011
137*
138*     .. Scalar Arguments ..
139      CHARACTER          DIAG, NORM, UPLO
140      INTEGER            INFO, N
141      REAL               RCOND
142*     ..
143*     .. Array Arguments ..
144      INTEGER            IWORK( * )
145      REAL               AP( * ), WORK( * )
146*     ..
147*
148*  =====================================================================
149*
150*     .. Parameters ..
151      REAL               ONE, ZERO
152      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
153*     ..
154*     .. Local Scalars ..
155      LOGICAL            NOUNIT, ONENRM, UPPER
156      CHARACTER          NORMIN
157      INTEGER            IX, KASE, KASE1
158      REAL               AINVNM, ANORM, SCALE, SMLNUM, XNORM
159*     ..
160*     .. Local Arrays ..
161      INTEGER            ISAVE( 3 )
162*     ..
163*     .. External Functions ..
164      LOGICAL            LSAME
165      INTEGER            ISAMAX
166      REAL               SLAMCH, SLANTP
167      EXTERNAL           LSAME, ISAMAX, SLAMCH, SLANTP
168*     ..
169*     .. External Subroutines ..
170      EXTERNAL           SLACN2, SLATPS, SRSCL, XERBLA
171*     ..
172*     .. Intrinsic Functions ..
173      INTRINSIC          ABS, MAX, REAL
174*     ..
175*     .. Executable Statements ..
176*
177*     Test the input parameters.
178*
179      INFO = 0
180      UPPER = LSAME( UPLO, 'U' )
181      ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
182      NOUNIT = LSAME( DIAG, 'N' )
183*
184      IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
185         INFO = -1
186      ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
187         INFO = -2
188      ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
189         INFO = -3
190      ELSE IF( N.LT.0 ) THEN
191         INFO = -4
192      END IF
193      IF( INFO.NE.0 ) THEN
194         CALL XERBLA( 'STPCON', -INFO )
195         RETURN
196      END IF
197*
198*     Quick return if possible
199*
200      IF( N.EQ.0 ) THEN
201         RCOND = ONE
202         RETURN
203      END IF
204*
205      RCOND = ZERO
206      SMLNUM = SLAMCH( 'Safe minimum' )*REAL( MAX( 1, N ) )
207*
208*     Compute the norm of the triangular matrix A.
209*
210      ANORM = SLANTP( NORM, UPLO, DIAG, N, AP, WORK )
211*
212*     Continue only if ANORM > 0.
213*
214      IF( ANORM.GT.ZERO ) THEN
215*
216*        Estimate the norm of the inverse of A.
217*
218         AINVNM = ZERO
219         NORMIN = 'N'
220         IF( ONENRM ) THEN
221            KASE1 = 1
222         ELSE
223            KASE1 = 2
224         END IF
225         KASE = 0
226   10    CONTINUE
227         CALL SLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
228         IF( KASE.NE.0 ) THEN
229            IF( KASE.EQ.KASE1 ) THEN
230*
231*              Multiply by inv(A).
232*
233               CALL SLATPS( UPLO, 'No transpose', DIAG, NORMIN, N, AP,
234     $                      WORK, SCALE, WORK( 2*N+1 ), INFO )
235            ELSE
236*
237*              Multiply by inv(A**T).
238*
239               CALL SLATPS( UPLO, 'Transpose', DIAG, NORMIN, N, AP,
240     $                      WORK, SCALE, WORK( 2*N+1 ), INFO )
241            END IF
242            NORMIN = 'Y'
243*
244*           Multiply by 1/SCALE if doing so will not cause overflow.
245*
246            IF( SCALE.NE.ONE ) THEN
247               IX = ISAMAX( N, WORK, 1 )
248               XNORM = ABS( WORK( IX ) )
249               IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
250     $            GO TO 20
251               CALL SRSCL( N, SCALE, WORK, 1 )
252            END IF
253            GO TO 10
254         END IF
255*
256*        Compute the estimate of the reciprocal condition number.
257*
258         IF( AINVNM.NE.ZERO )
259     $      RCOND = ( ONE / ANORM ) / AINVNM
260      END IF
261*
262   20 CONTINUE
263      RETURN
264*
265*     End of STPCON
266*
267      END
268