LAPACK 3.3.0

cpptrf.f

Go to the documentation of this file.
00001       SUBROUTINE CPPTRF( UPLO, N, AP, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          UPLO
00010       INTEGER            INFO, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       COMPLEX            AP( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  CPPTRF computes the Cholesky factorization of a complex Hermitian
00020 *  positive definite matrix A stored in packed format.
00021 *
00022 *  The factorization has the form
00023 *     A = U**H * U,  if UPLO = 'U', or
00024 *     A = L  * L**H,  if UPLO = 'L',
00025 *  where U is an upper triangular matrix and L is lower triangular.
00026 *
00027 *  Arguments
00028 *  =========
00029 *
00030 *  UPLO    (input) CHARACTER*1
00031 *          = 'U':  Upper triangle of A is stored;
00032 *          = 'L':  Lower triangle of A is stored.
00033 *
00034 *  N       (input) INTEGER
00035 *          The order of the matrix A.  N >= 0.
00036 *
00037 *  AP      (input/output) COMPLEX array, dimension (N*(N+1)/2)
00038 *          On entry, the upper or lower triangle of the Hermitian matrix
00039 *          A, packed columnwise in a linear array.  The j-th column of A
00040 *          is stored in the array AP as follows:
00041 *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
00042 *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
00043 *          See below for further details.
00044 *
00045 *          On exit, if INFO = 0, the triangular factor U or L from the
00046 *          Cholesky factorization A = U**H*U or A = L*L**H, in the same
00047 *          storage format as A.
00048 *
00049 *  INFO    (output) INTEGER
00050 *          = 0:  successful exit
00051 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00052 *          > 0:  if INFO = i, the leading minor of order i is not
00053 *                positive definite, and the factorization could not be
00054 *                completed.
00055 *
00056 *  Further Details
00057 *  ===============
00058 *
00059 *  The packed storage scheme is illustrated by the following example
00060 *  when N = 4, UPLO = 'U':
00061 *
00062 *  Two-dimensional storage of the Hermitian matrix A:
00063 *
00064 *     a11 a12 a13 a14
00065 *         a22 a23 a24
00066 *             a33 a34     (aij = conjg(aji))
00067 *                 a44
00068 *
00069 *  Packed storage of the upper triangle of A:
00070 *
00071 *  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
00072 *
00073 *  =====================================================================
00074 *
00075 *     .. Parameters ..
00076       REAL               ZERO, ONE
00077       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00078 *     ..
00079 *     .. Local Scalars ..
00080       LOGICAL            UPPER
00081       INTEGER            J, JC, JJ
00082       REAL               AJJ
00083 *     ..
00084 *     .. External Functions ..
00085       LOGICAL            LSAME
00086       COMPLEX            CDOTC
00087       EXTERNAL           LSAME, CDOTC
00088 *     ..
00089 *     .. External Subroutines ..
00090       EXTERNAL           CHPR, CSSCAL, CTPSV, XERBLA
00091 *     ..
00092 *     .. Intrinsic Functions ..
00093       INTRINSIC          REAL, SQRT
00094 *     ..
00095 *     .. Executable Statements ..
00096 *
00097 *     Test the input parameters.
00098 *
00099       INFO = 0
00100       UPPER = LSAME( UPLO, 'U' )
00101       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00102          INFO = -1
00103       ELSE IF( N.LT.0 ) THEN
00104          INFO = -2
00105       END IF
00106       IF( INFO.NE.0 ) THEN
00107          CALL XERBLA( 'CPPTRF', -INFO )
00108          RETURN
00109       END IF
00110 *
00111 *     Quick return if possible
00112 *
00113       IF( N.EQ.0 )
00114      $   RETURN
00115 *
00116       IF( UPPER ) THEN
00117 *
00118 *        Compute the Cholesky factorization A = U'*U.
00119 *
00120          JJ = 0
00121          DO 10 J = 1, N
00122             JC = JJ + 1
00123             JJ = JJ + J
00124 *
00125 *           Compute elements 1:J-1 of column J.
00126 *
00127             IF( J.GT.1 )
00128      $         CALL CTPSV( 'Upper', 'Conjugate transpose', 'Non-unit',
00129      $                     J-1, AP, AP( JC ), 1 )
00130 *
00131 *           Compute U(J,J) and test for non-positive-definiteness.
00132 *
00133             AJJ = REAL( AP( JJ ) ) - CDOTC( J-1, AP( JC ), 1, AP( JC ),
00134      $            1 )
00135             IF( AJJ.LE.ZERO ) THEN
00136                AP( JJ ) = AJJ
00137                GO TO 30
00138             END IF
00139             AP( JJ ) = SQRT( AJJ )
00140    10    CONTINUE
00141       ELSE
00142 *
00143 *        Compute the Cholesky factorization A = L*L'.
00144 *
00145          JJ = 1
00146          DO 20 J = 1, N
00147 *
00148 *           Compute L(J,J) and test for non-positive-definiteness.
00149 *
00150             AJJ = REAL( AP( JJ ) )
00151             IF( AJJ.LE.ZERO ) THEN
00152                AP( JJ ) = AJJ
00153                GO TO 30
00154             END IF
00155             AJJ = SQRT( AJJ )
00156             AP( JJ ) = AJJ
00157 *
00158 *           Compute elements J+1:N of column J and update the trailing
00159 *           submatrix.
00160 *
00161             IF( J.LT.N ) THEN
00162                CALL CSSCAL( N-J, ONE / AJJ, AP( JJ+1 ), 1 )
00163                CALL CHPR( 'Lower', N-J, -ONE, AP( JJ+1 ), 1,
00164      $                    AP( JJ+N-J+1 ) )
00165                JJ = JJ + N - J + 1
00166             END IF
00167    20    CONTINUE
00168       END IF
00169       GO TO 40
00170 *
00171    30 CONTINUE
00172       INFO = J
00173 *
00174    40 CONTINUE
00175       RETURN
00176 *
00177 *     End of CPPTRF
00178 *
00179       END
 All Files Functions