![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZPPTRF 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZPPTRF + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zpptrf.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zpptrf.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zpptrf.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZPPTRF( UPLO, N, AP, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER UPLO 00025 * INTEGER INFO, N 00026 * .. 00027 * .. Array Arguments .. 00028 * COMPLEX*16 AP( * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> ZPPTRF computes the Cholesky factorization of a complex Hermitian 00038 *> positive definite matrix A stored in packed format. 00039 *> 00040 *> The factorization has the form 00041 *> A = U**H * U, if UPLO = 'U', or 00042 *> A = L * L**H, if UPLO = 'L', 00043 *> where U is an upper triangular matrix and L is lower triangular. 00044 *> \endverbatim 00045 * 00046 * Arguments: 00047 * ========== 00048 * 00049 *> \param[in] UPLO 00050 *> \verbatim 00051 *> UPLO is CHARACTER*1 00052 *> = 'U': Upper triangle of A is stored; 00053 *> = 'L': Lower triangle of A is stored. 00054 *> \endverbatim 00055 *> 00056 *> \param[in] N 00057 *> \verbatim 00058 *> N is INTEGER 00059 *> The order of the matrix A. N >= 0. 00060 *> \endverbatim 00061 *> 00062 *> \param[in,out] AP 00063 *> \verbatim 00064 *> AP is COMPLEX*16 array, dimension (N*(N+1)/2) 00065 *> On entry, the upper or lower triangle of the Hermitian matrix 00066 *> A, packed columnwise in a linear array. The j-th column of A 00067 *> is stored in the array AP as follows: 00068 *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 00069 *> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. 00070 *> See below for further details. 00071 *> 00072 *> On exit, if INFO = 0, the triangular factor U or L from the 00073 *> Cholesky factorization A = U**H*U or A = L*L**H, in the same 00074 *> storage format as A. 00075 *> \endverbatim 00076 *> 00077 *> \param[out] INFO 00078 *> \verbatim 00079 *> INFO is INTEGER 00080 *> = 0: successful exit 00081 *> < 0: if INFO = -i, the i-th argument had an illegal value 00082 *> > 0: if INFO = i, the leading minor of order i is not 00083 *> positive definite, and the factorization could not be 00084 *> completed. 00085 *> \endverbatim 00086 * 00087 * Authors: 00088 * ======== 00089 * 00090 *> \author Univ. of Tennessee 00091 *> \author Univ. of California Berkeley 00092 *> \author Univ. of Colorado Denver 00093 *> \author NAG Ltd. 00094 * 00095 *> \date November 2011 00096 * 00097 *> \ingroup complex16OTHERcomputational 00098 * 00099 *> \par Further Details: 00100 * ===================== 00101 *> 00102 *> \verbatim 00103 *> 00104 *> The packed storage scheme is illustrated by the following example 00105 *> when N = 4, UPLO = 'U': 00106 *> 00107 *> Two-dimensional storage of the Hermitian matrix A: 00108 *> 00109 *> a11 a12 a13 a14 00110 *> a22 a23 a24 00111 *> a33 a34 (aij = conjg(aji)) 00112 *> a44 00113 *> 00114 *> Packed storage of the upper triangle of A: 00115 *> 00116 *> AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ] 00117 *> \endverbatim 00118 *> 00119 * ===================================================================== 00120 SUBROUTINE ZPPTRF( UPLO, N, AP, INFO ) 00121 * 00122 * -- LAPACK computational routine (version 3.4.0) -- 00123 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00124 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00125 * November 2011 00126 * 00127 * .. Scalar Arguments .. 00128 CHARACTER UPLO 00129 INTEGER INFO, N 00130 * .. 00131 * .. Array Arguments .. 00132 COMPLEX*16 AP( * ) 00133 * .. 00134 * 00135 * ===================================================================== 00136 * 00137 * .. Parameters .. 00138 DOUBLE PRECISION ZERO, ONE 00139 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00140 * .. 00141 * .. Local Scalars .. 00142 LOGICAL UPPER 00143 INTEGER J, JC, JJ 00144 DOUBLE PRECISION AJJ 00145 * .. 00146 * .. External Functions .. 00147 LOGICAL LSAME 00148 COMPLEX*16 ZDOTC 00149 EXTERNAL LSAME, ZDOTC 00150 * .. 00151 * .. External Subroutines .. 00152 EXTERNAL XERBLA, ZDSCAL, ZHPR, ZTPSV 00153 * .. 00154 * .. Intrinsic Functions .. 00155 INTRINSIC DBLE, SQRT 00156 * .. 00157 * .. Executable Statements .. 00158 * 00159 * Test the input parameters. 00160 * 00161 INFO = 0 00162 UPPER = LSAME( UPLO, 'U' ) 00163 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00164 INFO = -1 00165 ELSE IF( N.LT.0 ) THEN 00166 INFO = -2 00167 END IF 00168 IF( INFO.NE.0 ) THEN 00169 CALL XERBLA( 'ZPPTRF', -INFO ) 00170 RETURN 00171 END IF 00172 * 00173 * Quick return if possible 00174 * 00175 IF( N.EQ.0 ) 00176 $ RETURN 00177 * 00178 IF( UPPER ) THEN 00179 * 00180 * Compute the Cholesky factorization A = U**H * U. 00181 * 00182 JJ = 0 00183 DO 10 J = 1, N 00184 JC = JJ + 1 00185 JJ = JJ + J 00186 * 00187 * Compute elements 1:J-1 of column J. 00188 * 00189 IF( J.GT.1 ) 00190 $ CALL ZTPSV( 'Upper', 'Conjugate transpose', 'Non-unit', 00191 $ J-1, AP, AP( JC ), 1 ) 00192 * 00193 * Compute U(J,J) and test for non-positive-definiteness. 00194 * 00195 AJJ = DBLE( AP( JJ ) ) - ZDOTC( J-1, AP( JC ), 1, AP( JC ), 00196 $ 1 ) 00197 IF( AJJ.LE.ZERO ) THEN 00198 AP( JJ ) = AJJ 00199 GO TO 30 00200 END IF 00201 AP( JJ ) = SQRT( AJJ ) 00202 10 CONTINUE 00203 ELSE 00204 * 00205 * Compute the Cholesky factorization A = L * L**H. 00206 * 00207 JJ = 1 00208 DO 20 J = 1, N 00209 * 00210 * Compute L(J,J) and test for non-positive-definiteness. 00211 * 00212 AJJ = DBLE( AP( JJ ) ) 00213 IF( AJJ.LE.ZERO ) THEN 00214 AP( JJ ) = AJJ 00215 GO TO 30 00216 END IF 00217 AJJ = SQRT( AJJ ) 00218 AP( JJ ) = AJJ 00219 * 00220 * Compute elements J+1:N of column J and update the trailing 00221 * submatrix. 00222 * 00223 IF( J.LT.N ) THEN 00224 CALL ZDSCAL( N-J, ONE / AJJ, AP( JJ+1 ), 1 ) 00225 CALL ZHPR( 'Lower', N-J, -ONE, AP( JJ+1 ), 1, 00226 $ AP( JJ+N-J+1 ) ) 00227 JJ = JJ + N - J + 1 00228 END IF 00229 20 CONTINUE 00230 END IF 00231 GO TO 40 00232 * 00233 30 CONTINUE 00234 INFO = J 00235 * 00236 40 CONTINUE 00237 RETURN 00238 * 00239 * End of ZPPTRF 00240 * 00241 END