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