1*> \brief \b CPSTF2 computes the Cholesky factorization with complete pivoting of complex Hermitian positive semidefinite matrix.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CPSTF2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cpstf2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cpstf2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cpstf2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CPSTF2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
22*
23*       .. Scalar Arguments ..
24*       REAL               TOL
25*       INTEGER            INFO, LDA, N, RANK
26*       CHARACTER          UPLO
27*       ..
28*       .. Array Arguments ..
29*       COMPLEX            A( LDA, * )
30*       REAL               WORK( 2*N )
31*       INTEGER            PIV( N )
32*       ..
33*
34*
35*> \par Purpose:
36*  =============
37*>
38*> \verbatim
39*>
40*> CPSTF2 computes the Cholesky factorization with complete
41*> pivoting of a complex Hermitian positive semidefinite matrix A.
42*>
43*> The factorization has the form
44*>    P**T * A * P = U**H * U ,  if UPLO = 'U',
45*>    P**T * A * P = L  * L**H,  if UPLO = 'L',
46*> where U is an upper triangular matrix and L is lower triangular, and
47*> P is stored as vector PIV.
48*>
49*> This algorithm does not attempt to check that A is positive
50*> semidefinite. This version of the algorithm calls level 2 BLAS.
51*> \endverbatim
52*
53*  Arguments:
54*  ==========
55*
56*> \param[in] UPLO
57*> \verbatim
58*>          UPLO is CHARACTER*1
59*>          Specifies whether the upper or lower triangular part of the
60*>          symmetric matrix A is stored.
61*>          = 'U':  Upper triangular
62*>          = 'L':  Lower triangular
63*> \endverbatim
64*>
65*> \param[in] N
66*> \verbatim
67*>          N is INTEGER
68*>          The order of the matrix A.  N >= 0.
69*> \endverbatim
70*>
71*> \param[in,out] A
72*> \verbatim
73*>          A is COMPLEX array, dimension (LDA,N)
74*>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
75*>          n by n upper triangular part of A contains the upper
76*>          triangular part of the matrix A, and the strictly lower
77*>          triangular part of A is not referenced.  If UPLO = 'L', the
78*>          leading n by n lower triangular part of A contains the lower
79*>          triangular part of the matrix A, and the strictly upper
80*>          triangular part of A is not referenced.
81*>
82*>          On exit, if INFO = 0, the factor U or L from the Cholesky
83*>          factorization as above.
84*> \endverbatim
85*>
86*> \param[out] PIV
87*> \verbatim
88*>          PIV is INTEGER array, dimension (N)
89*>          PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
90*> \endverbatim
91*>
92*> \param[out] RANK
93*> \verbatim
94*>          RANK is INTEGER
95*>          The rank of A given by the number of steps the algorithm
96*>          completed.
97*> \endverbatim
98*>
99*> \param[in] TOL
100*> \verbatim
101*>          TOL is REAL
102*>          User defined tolerance. If TOL < 0, then N*U*MAX( A( K,K ) )
103*>          will be used. The algorithm terminates at the (K-1)st step
104*>          if the pivot <= TOL.
105*> \endverbatim
106*>
107*> \param[in] LDA
108*> \verbatim
109*>          LDA is INTEGER
110*>          The leading dimension of the array A.  LDA >= max(1,N).
111*> \endverbatim
112*>
113*> \param[out] WORK
114*> \verbatim
115*>          WORK is REAL array, dimension (2*N)
116*>          Work space.
117*> \endverbatim
118*>
119*> \param[out] INFO
120*> \verbatim
121*>          INFO is INTEGER
122*>          < 0: If INFO = -K, the K-th argument had an illegal value,
123*>          = 0: algorithm completed successfully, and
124*>          > 0: the matrix A is either rank deficient with computed rank
125*>               as returned in RANK, or is not positive semidefinite. See
126*>               Section 7 of LAPACK Working Note #161 for further
127*>               information.
128*> \endverbatim
129*
130*  Authors:
131*  ========
132*
133*> \author Univ. of Tennessee
134*> \author Univ. of California Berkeley
135*> \author Univ. of Colorado Denver
136*> \author NAG Ltd.
137*
138*> \ingroup complexOTHERcomputational
139*
140*  =====================================================================
141      SUBROUTINE CPSTF2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
142*
143*  -- LAPACK computational routine --
144*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
145*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
146*
147*     .. Scalar Arguments ..
148      REAL               TOL
149      INTEGER            INFO, LDA, N, RANK
150      CHARACTER          UPLO
151*     ..
152*     .. Array Arguments ..
153      COMPLEX            A( LDA, * )
154      REAL               WORK( 2*N )
155      INTEGER            PIV( N )
156*     ..
157*
158*  =====================================================================
159*
160*     .. Parameters ..
161      REAL               ONE, ZERO
162      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
163      COMPLEX            CONE
164      PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
165*     ..
166*     .. Local Scalars ..
167      COMPLEX            CTEMP
168      REAL               AJJ, SSTOP, STEMP
169      INTEGER            I, ITEMP, J, PVT
170      LOGICAL            UPPER
171*     ..
172*     .. External Functions ..
173      REAL               SLAMCH
174      LOGICAL            LSAME, SISNAN
175      EXTERNAL           SLAMCH, LSAME, SISNAN
176*     ..
177*     .. External Subroutines ..
178      EXTERNAL           CGEMV, CLACGV, CSSCAL, CSWAP, XERBLA
179*     ..
180*     .. Intrinsic Functions ..
181      INTRINSIC          CONJG, MAX, REAL, SQRT
182*     ..
183*     .. Executable Statements ..
184*
185*     Test the input parameters
186*
187      INFO = 0
188      UPPER = LSAME( UPLO, 'U' )
189      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
190         INFO = -1
191      ELSE IF( N.LT.0 ) THEN
192         INFO = -2
193      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
194         INFO = -4
195      END IF
196      IF( INFO.NE.0 ) THEN
197         CALL XERBLA( 'CPSTF2', -INFO )
198         RETURN
199      END IF
200*
201*     Quick return if possible
202*
203      IF( N.EQ.0 )
204     $   RETURN
205*
206*     Initialize PIV
207*
208      DO 100 I = 1, N
209         PIV( I ) = I
210  100 CONTINUE
211*
212*     Compute stopping value
213*
214      DO 110 I = 1, N
215         WORK( I ) = REAL( A( I, I ) )
216  110 CONTINUE
217      PVT = MAXLOC( WORK( 1:N ), 1 )
218      AJJ = REAL ( A( PVT, PVT ) )
219      IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN
220         RANK = 0
221         INFO = 1
222         GO TO 200
223      END IF
224*
225*     Compute stopping value if not supplied
226*
227      IF( TOL.LT.ZERO ) THEN
228         SSTOP = N * SLAMCH( 'Epsilon' ) * AJJ
229      ELSE
230         SSTOP = TOL
231      END IF
232*
233*     Set first half of WORK to zero, holds dot products
234*
235      DO 120 I = 1, N
236         WORK( I ) = 0
237  120 CONTINUE
238*
239      IF( UPPER ) THEN
240*
241*        Compute the Cholesky factorization P**T * A * P = U**H * U
242*
243         DO 150 J = 1, N
244*
245*        Find pivot, test for exit, else swap rows and columns
246*        Update dot products, compute possible pivots which are
247*        stored in the second half of WORK
248*
249            DO 130 I = J, N
250*
251               IF( J.GT.1 ) THEN
252                  WORK( I ) = WORK( I ) +
253     $                        REAL( CONJG( A( J-1, I ) )*
254     $                              A( J-1, I ) )
255               END IF
256               WORK( N+I ) = REAL( A( I, I ) ) - WORK( I )
257*
258  130       CONTINUE
259*
260            IF( J.GT.1 ) THEN
261               ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
262               PVT = ITEMP + J - 1
263               AJJ = WORK( N+PVT )
264               IF( AJJ.LE.SSTOP.OR.SISNAN( AJJ ) ) THEN
265                  A( J, J ) = AJJ
266                  GO TO 190
267               END IF
268            END IF
269*
270            IF( J.NE.PVT ) THEN
271*
272*              Pivot OK, so can now swap pivot rows and columns
273*
274               A( PVT, PVT ) = A( J, J )
275               CALL CSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
276               IF( PVT.LT.N )
277     $            CALL CSWAP( N-PVT, A( J, PVT+1 ), LDA,
278     $                        A( PVT, PVT+1 ), LDA )
279               DO 140 I = J + 1, PVT - 1
280                  CTEMP = CONJG( A( J, I ) )
281                  A( J, I ) = CONJG( A( I, PVT ) )
282                  A( I, PVT ) = CTEMP
283  140          CONTINUE
284               A( J, PVT ) = CONJG( A( J, PVT ) )
285*
286*              Swap dot products and PIV
287*
288               STEMP = WORK( J )
289               WORK( J ) = WORK( PVT )
290               WORK( PVT ) = STEMP
291               ITEMP = PIV( PVT )
292               PIV( PVT ) = PIV( J )
293               PIV( J ) = ITEMP
294            END IF
295*
296            AJJ = SQRT( AJJ )
297            A( J, J ) = AJJ
298*
299*           Compute elements J+1:N of row J
300*
301            IF( J.LT.N ) THEN
302               CALL CLACGV( J-1, A( 1, J ), 1 )
303               CALL CGEMV( 'Trans', J-1, N-J, -CONE, A( 1, J+1 ), LDA,
304     $                     A( 1, J ), 1, CONE, A( J, J+1 ), LDA )
305               CALL CLACGV( J-1, A( 1, J ), 1 )
306               CALL CSSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
307            END IF
308*
309  150    CONTINUE
310*
311      ELSE
312*
313*        Compute the Cholesky factorization P**T * A * P = L * L**H
314*
315         DO 180 J = 1, N
316*
317*        Find pivot, test for exit, else swap rows and columns
318*        Update dot products, compute possible pivots which are
319*        stored in the second half of WORK
320*
321            DO 160 I = J, N
322*
323               IF( J.GT.1 ) THEN
324                  WORK( I ) = WORK( I ) +
325     $                        REAL( CONJG( A( I, J-1 ) )*
326     $                              A( I, J-1 ) )
327               END IF
328               WORK( N+I ) = REAL( A( I, I ) ) - WORK( I )
329*
330  160       CONTINUE
331*
332            IF( J.GT.1 ) THEN
333               ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
334               PVT = ITEMP + J - 1
335               AJJ = WORK( N+PVT )
336               IF( AJJ.LE.SSTOP.OR.SISNAN( AJJ ) ) THEN
337                  A( J, J ) = AJJ
338                  GO TO 190
339               END IF
340            END IF
341*
342            IF( J.NE.PVT ) THEN
343*
344*              Pivot OK, so can now swap pivot rows and columns
345*
346               A( PVT, PVT ) = A( J, J )
347               CALL CSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
348               IF( PVT.LT.N )
349     $            CALL CSWAP( N-PVT, A( PVT+1, J ), 1, A( PVT+1, PVT ),
350     $                        1 )
351               DO 170 I = J + 1, PVT - 1
352                  CTEMP = CONJG( A( I, J ) )
353                  A( I, J ) = CONJG( A( PVT, I ) )
354                  A( PVT, I ) = CTEMP
355  170          CONTINUE
356               A( PVT, J ) = CONJG( A( PVT, J ) )
357*
358*              Swap dot products and PIV
359*
360               STEMP = WORK( J )
361               WORK( J ) = WORK( PVT )
362               WORK( PVT ) = STEMP
363               ITEMP = PIV( PVT )
364               PIV( PVT ) = PIV( J )
365               PIV( J ) = ITEMP
366            END IF
367*
368            AJJ = SQRT( AJJ )
369            A( J, J ) = AJJ
370*
371*           Compute elements J+1:N of column J
372*
373            IF( J.LT.N ) THEN
374               CALL CLACGV( J-1, A( J, 1 ), LDA )
375               CALL CGEMV( 'No Trans', N-J, J-1, -CONE, A( J+1, 1 ),
376     $                     LDA, A( J, 1 ), LDA, CONE, A( J+1, J ), 1 )
377               CALL CLACGV( J-1, A( J, 1 ), LDA )
378               CALL CSSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
379            END IF
380*
381  180    CONTINUE
382*
383      END IF
384*
385*     Ran to completion, A has full rank
386*
387      RANK = N
388*
389      GO TO 200
390  190 CONTINUE
391*
392*     Rank is number of steps completed.  Set INFO = 1 to signal
393*     that the factorization cannot be used to solve a system.
394*
395      RANK = J - 1
396      INFO = 1
397*
398  200 CONTINUE
399      RETURN
400*
401*     End of CPSTF2
402*
403      END
404