1*> \brief \b ZPPTRF
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZPPTRF + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zpptrf.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zpptrf.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zpptrf.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE ZPPTRF( UPLO, N, AP, INFO )
22*
23*       .. Scalar Arguments ..
24*       CHARACTER          UPLO
25*       INTEGER            INFO, N
26*       ..
27*       .. Array Arguments ..
28*       COMPLEX*16         AP( * )
29*       ..
30*
31*
32*> \par Purpose:
33*  =============
34*>
35*> \verbatim
36*>
37*> ZPPTRF computes the Cholesky factorization of a complex Hermitian
38*> positive definite matrix A stored in packed format.
39*>
40*> The factorization has the form
41*>    A = U**H * U,  if UPLO = 'U', or
42*>    A = L  * L**H,  if UPLO = 'L',
43*> where U is an upper triangular matrix and L is lower triangular.
44*> \endverbatim
45*
46*  Arguments:
47*  ==========
48*
49*> \param[in] UPLO
50*> \verbatim
51*>          UPLO is CHARACTER*1
52*>          = 'U':  Upper triangle of A is stored;
53*>          = 'L':  Lower triangle of A is stored.
54*> \endverbatim
55*>
56*> \param[in] N
57*> \verbatim
58*>          N is INTEGER
59*>          The order of the matrix A.  N >= 0.
60*> \endverbatim
61*>
62*> \param[in,out] AP
63*> \verbatim
64*>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
65*>          On entry, the upper or lower triangle of the Hermitian matrix
66*>          A, packed columnwise in a linear array.  The j-th column of A
67*>          is stored in the array AP as follows:
68*>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
69*>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
70*>          See below for further details.
71*>
72*>          On exit, if INFO = 0, the triangular factor U or L from the
73*>          Cholesky factorization A = U**H*U or A = L*L**H, in the same
74*>          storage format as A.
75*> \endverbatim
76*>
77*> \param[out] INFO
78*> \verbatim
79*>          INFO is INTEGER
80*>          = 0:  successful exit
81*>          < 0:  if INFO = -i, the i-th argument had an illegal value
82*>          > 0:  if INFO = i, the leading minor of order i is not
83*>                positive definite, and the factorization could not be
84*>                completed.
85*> \endverbatim
86*
87*  Authors:
88*  ========
89*
90*> \author Univ. of Tennessee
91*> \author Univ. of California Berkeley
92*> \author Univ. of Colorado Denver
93*> \author NAG Ltd.
94*
95*> \ingroup complex16OTHERcomputational
96*
97*> \par Further Details:
98*  =====================
99*>
100*> \verbatim
101*>
102*>  The packed storage scheme is illustrated by the following example
103*>  when N = 4, UPLO = 'U':
104*>
105*>  Two-dimensional storage of the Hermitian matrix A:
106*>
107*>     a11 a12 a13 a14
108*>         a22 a23 a24
109*>             a33 a34     (aij = conjg(aji))
110*>                 a44
111*>
112*>  Packed storage of the upper triangle of A:
113*>
114*>  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
115*> \endverbatim
116*>
117*  =====================================================================
118      SUBROUTINE ZPPTRF( UPLO, N, AP, INFO )
119*
120*  -- LAPACK computational routine --
121*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
122*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
123*
124*     .. Scalar Arguments ..
125      CHARACTER          UPLO
126      INTEGER            INFO, N
127*     ..
128*     .. Array Arguments ..
129      COMPLEX*16         AP( * )
130*     ..
131*
132*  =====================================================================
133*
134*     .. Parameters ..
135      DOUBLE PRECISION   ZERO, ONE
136      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
137*     ..
138*     .. Local Scalars ..
139      LOGICAL            UPPER
140      INTEGER            J, JC, JJ
141      DOUBLE PRECISION   AJJ
142*     ..
143*     .. External Functions ..
144      LOGICAL            LSAME
145      COMPLEX*16         ZDOTC
146      EXTERNAL           LSAME, ZDOTC
147*     ..
148*     .. External Subroutines ..
149      EXTERNAL           XERBLA, ZDSCAL, ZHPR, ZTPSV
150*     ..
151*     .. Intrinsic Functions ..
152      INTRINSIC          DBLE, SQRT
153*     ..
154*     .. Executable Statements ..
155*
156*     Test the input parameters.
157*
158      INFO = 0
159      UPPER = LSAME( UPLO, 'U' )
160      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
161         INFO = -1
162      ELSE IF( N.LT.0 ) THEN
163         INFO = -2
164      END IF
165      IF( INFO.NE.0 ) THEN
166         CALL XERBLA( 'ZPPTRF', -INFO )
167         RETURN
168      END IF
169*
170*     Quick return if possible
171*
172      IF( N.EQ.0 )
173     $   RETURN
174*
175      IF( UPPER ) THEN
176*
177*        Compute the Cholesky factorization A = U**H * U.
178*
179         JJ = 0
180         DO 10 J = 1, N
181            JC = JJ + 1
182            JJ = JJ + J
183*
184*           Compute elements 1:J-1 of column J.
185*
186            IF( J.GT.1 )
187     $         CALL ZTPSV( 'Upper', 'Conjugate transpose', 'Non-unit',
188     $                     J-1, AP, AP( JC ), 1 )
189*
190*           Compute U(J,J) and test for non-positive-definiteness.
191*
192            AJJ = DBLE( AP( JJ ) ) - DBLE( ZDOTC( J-1,
193     $            AP( JC ), 1, AP( JC ), 1 ) )
194            IF( AJJ.LE.ZERO ) THEN
195               AP( JJ ) = AJJ
196               GO TO 30
197            END IF
198            AP( JJ ) = SQRT( AJJ )
199   10    CONTINUE
200      ELSE
201*
202*        Compute the Cholesky factorization A = L * L**H.
203*
204         JJ = 1
205         DO 20 J = 1, N
206*
207*           Compute L(J,J) and test for non-positive-definiteness.
208*
209            AJJ = DBLE( AP( JJ ) )
210            IF( AJJ.LE.ZERO ) THEN
211               AP( JJ ) = AJJ
212               GO TO 30
213            END IF
214            AJJ = SQRT( AJJ )
215            AP( JJ ) = AJJ
216*
217*           Compute elements J+1:N of column J and update the trailing
218*           submatrix.
219*
220            IF( J.LT.N ) THEN
221               CALL ZDSCAL( N-J, ONE / AJJ, AP( JJ+1 ), 1 )
222               CALL ZHPR( 'Lower', N-J, -ONE, AP( JJ+1 ), 1,
223     $                    AP( JJ+N-J+1 ) )
224               JJ = JJ + N - J + 1
225            END IF
226   20    CONTINUE
227      END IF
228      GO TO 40
229*
230   30 CONTINUE
231      INFO = J
232*
233   40 CONTINUE
234      RETURN
235*
236*     End of ZPPTRF
237*
238      END
239