3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download SPPTRF + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/spptrf.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/spptrf.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/spptrf.f">
21 * SUBROUTINE SPPTRF( UPLO, N, AP, INFO )
23 * .. Scalar Arguments ..
27 * .. Array Arguments ..
37 *> SPPTRF computes the Cholesky factorization of a real symmetric
38 *> positive definite matrix A stored in packed format.
40 *> The factorization has the form
41 *> A = U**T * U, if UPLO = 'U', or
42 *> A = L * L**T, if UPLO = 'L',
43 *> where U is an upper triangular matrix and L is lower triangular.
51 *> UPLO is CHARACTER*1
52 *> = 'U': Upper triangle of A is stored;
53 *> = 'L': Lower triangle of A is stored.
59 *> The order of the matrix A. N >= 0.
64 *> AP is REAL array, dimension (N*(N+1)/2)
65 *> On entry, the upper or lower triangle of the symmetric 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.
72 *> On exit, if INFO = 0, the triangular factor U or L from the
73 *> Cholesky factorization A = U**T*U or A = L*L**T, in the same
74 *> storage format as A.
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
90 *> \author Univ. of Tennessee
91 *> \author Univ. of California Berkeley
92 *> \author Univ. of Colorado Denver
95 *> \date November 2011
97 *> \ingroup realOTHERcomputational
99 *> \par Further Details:
100 * =====================
104 *> The packed storage scheme is illustrated by the following example
105 *> when N = 4, UPLO = 'U':
107 *> Two-dimensional storage of the symmetric matrix A:
111 *> a33 a34 (aij = aji)
114 *> Packed storage of the upper triangle of A:
116 *> AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
119 * =====================================================================
120 SUBROUTINE SPPTRF( UPLO, N, AP, INFO )
122 * -- LAPACK computational routine (version 3.4.0) --
123 * -- LAPACK is a software package provided by Univ. of Tennessee, --
124 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
127 * .. Scalar Arguments ..
131 * .. Array Arguments ..
135 * =====================================================================
139 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
141 * .. Local Scalars ..
146 * .. External Functions ..
151 * .. External Subroutines ..
152 EXTERNAL SSCAL, SSPR, STPSV, XERBLA
154 * .. Intrinsic Functions ..
157 * .. Executable Statements ..
159 * Test the input parameters.
162 UPPER = LSAME( UPLO, 'U' )
163 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
165 ELSE IF( N.LT.0 ) THEN
169 CALL XERBLA( 'SPPTRF', -INFO )
173 * Quick return if possible
180 * Compute the Cholesky factorization A = U**T*U.
187 * Compute elements 1:J-1 of column J.
190 $ CALL STPSV( 'Upper', 'Transpose', 'Non-unit', J-1, AP,
193 * Compute U(J,J) and test for non-positive-definiteness.
195 AJJ = AP( JJ ) - SDOT( J-1, AP( JC ), 1, AP( JC ), 1 )
196 IF( AJJ.LE.ZERO ) THEN
200 AP( JJ ) = SQRT( AJJ )
204 * Compute the Cholesky factorization A = L*L**T.
209 * Compute L(J,J) and test for non-positive-definiteness.
212 IF( AJJ.LE.ZERO ) THEN
219 * Compute elements J+1:N of column J and update the trailing
223 CALL SSCAL( N-J, ONE / AJJ, AP( JJ+1 ), 1 )
224 CALL SSPR( 'Lower', N-J, -ONE, AP( JJ+1 ), 1,