![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CPPTRI 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CPPTRI + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cpptri.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cpptri.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cpptri.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CPPTRI( UPLO, N, AP, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER UPLO 00025 * INTEGER INFO, N 00026 * .. 00027 * .. Array Arguments .. 00028 * COMPLEX AP( * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> CPPTRI computes the inverse of a complex Hermitian positive definite 00038 *> matrix A using the Cholesky factorization A = U**H*U or A = L*L**H 00039 *> computed by CPPTRF. 00040 *> \endverbatim 00041 * 00042 * Arguments: 00043 * ========== 00044 * 00045 *> \param[in] UPLO 00046 *> \verbatim 00047 *> UPLO is CHARACTER*1 00048 *> = 'U': Upper triangular factor is stored in AP; 00049 *> = 'L': Lower triangular factor is stored in AP. 00050 *> \endverbatim 00051 *> 00052 *> \param[in] N 00053 *> \verbatim 00054 *> N is INTEGER 00055 *> The order of the matrix A. N >= 0. 00056 *> \endverbatim 00057 *> 00058 *> \param[in,out] AP 00059 *> \verbatim 00060 *> AP is COMPLEX array, dimension (N*(N+1)/2) 00061 *> On entry, the triangular factor U or L from the Cholesky 00062 *> factorization A = U**H*U or A = L*L**H, packed columnwise as 00063 *> a linear array. The j-th column of U or L is stored in the 00064 *> array AP as follows: 00065 *> if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j; 00066 *> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n. 00067 *> 00068 *> On exit, the upper or lower triangle of the (Hermitian) 00069 *> inverse of A, overwriting the input factor U or L. 00070 *> \endverbatim 00071 *> 00072 *> \param[out] INFO 00073 *> \verbatim 00074 *> INFO is INTEGER 00075 *> = 0: successful exit 00076 *> < 0: if INFO = -i, the i-th argument had an illegal value 00077 *> > 0: if INFO = i, the (i,i) element of the factor U or L is 00078 *> zero, and the inverse could not be computed. 00079 *> \endverbatim 00080 * 00081 * Authors: 00082 * ======== 00083 * 00084 *> \author Univ. of Tennessee 00085 *> \author Univ. of California Berkeley 00086 *> \author Univ. of Colorado Denver 00087 *> \author NAG Ltd. 00088 * 00089 *> \date November 2011 00090 * 00091 *> \ingroup complexOTHERcomputational 00092 * 00093 * ===================================================================== 00094 SUBROUTINE CPPTRI( UPLO, N, AP, INFO ) 00095 * 00096 * -- LAPACK computational routine (version 3.4.0) -- 00097 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00098 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00099 * November 2011 00100 * 00101 * .. Scalar Arguments .. 00102 CHARACTER UPLO 00103 INTEGER INFO, N 00104 * .. 00105 * .. Array Arguments .. 00106 COMPLEX AP( * ) 00107 * .. 00108 * 00109 * ===================================================================== 00110 * 00111 * .. Parameters .. 00112 REAL ONE 00113 PARAMETER ( ONE = 1.0E+0 ) 00114 * .. 00115 * .. Local Scalars .. 00116 LOGICAL UPPER 00117 INTEGER J, JC, JJ, JJN 00118 REAL AJJ 00119 * .. 00120 * .. External Functions .. 00121 LOGICAL LSAME 00122 COMPLEX CDOTC 00123 EXTERNAL LSAME, CDOTC 00124 * .. 00125 * .. External Subroutines .. 00126 EXTERNAL CHPR, CSSCAL, CTPMV, CTPTRI, XERBLA 00127 * .. 00128 * .. Intrinsic Functions .. 00129 INTRINSIC REAL 00130 * .. 00131 * .. Executable Statements .. 00132 * 00133 * Test the input parameters. 00134 * 00135 INFO = 0 00136 UPPER = LSAME( UPLO, 'U' ) 00137 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00138 INFO = -1 00139 ELSE IF( N.LT.0 ) THEN 00140 INFO = -2 00141 END IF 00142 IF( INFO.NE.0 ) THEN 00143 CALL XERBLA( 'CPPTRI', -INFO ) 00144 RETURN 00145 END IF 00146 * 00147 * Quick return if possible 00148 * 00149 IF( N.EQ.0 ) 00150 $ RETURN 00151 * 00152 * Invert the triangular Cholesky factor U or L. 00153 * 00154 CALL CTPTRI( UPLO, 'Non-unit', N, AP, INFO ) 00155 IF( INFO.GT.0 ) 00156 $ RETURN 00157 IF( UPPER ) THEN 00158 * 00159 * Compute the product inv(U) * inv(U)**H. 00160 * 00161 JJ = 0 00162 DO 10 J = 1, N 00163 JC = JJ + 1 00164 JJ = JJ + J 00165 IF( J.GT.1 ) 00166 $ CALL CHPR( 'Upper', J-1, ONE, AP( JC ), 1, AP ) 00167 AJJ = AP( JJ ) 00168 CALL CSSCAL( J, AJJ, AP( JC ), 1 ) 00169 10 CONTINUE 00170 * 00171 ELSE 00172 * 00173 * Compute the product inv(L)**H * inv(L). 00174 * 00175 JJ = 1 00176 DO 20 J = 1, N 00177 JJN = JJ + N - J + 1 00178 AP( JJ ) = REAL( CDOTC( N-J+1, AP( JJ ), 1, AP( JJ ), 1 ) ) 00179 IF( J.LT.N ) 00180 $ CALL CTPMV( 'Lower', 'Conjugate transpose', 'Non-unit', 00181 $ N-J, AP( JJN ), AP( JJ+1 ), 1 ) 00182 JJ = JJN 00183 20 CONTINUE 00184 END IF 00185 * 00186 RETURN 00187 * 00188 * End of CPPTRI 00189 * 00190 END