![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b STPTRI 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download STPTRI + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stptri.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stptri.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stptri.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE STPTRI( UPLO, DIAG, N, AP, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER DIAG, UPLO 00025 * INTEGER INFO, N 00026 * .. 00027 * .. Array Arguments .. 00028 * REAL AP( * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> STPTRI computes the inverse of a real upper or lower triangular 00038 *> matrix A stored in packed format. 00039 *> \endverbatim 00040 * 00041 * Arguments: 00042 * ========== 00043 * 00044 *> \param[in] UPLO 00045 *> \verbatim 00046 *> UPLO is CHARACTER*1 00047 *> = 'U': A is upper triangular; 00048 *> = 'L': A is lower triangular. 00049 *> \endverbatim 00050 *> 00051 *> \param[in] DIAG 00052 *> \verbatim 00053 *> DIAG is CHARACTER*1 00054 *> = 'N': A is non-unit triangular; 00055 *> = 'U': A is unit triangular. 00056 *> \endverbatim 00057 *> 00058 *> \param[in] N 00059 *> \verbatim 00060 *> N is INTEGER 00061 *> The order of the matrix A. N >= 0. 00062 *> \endverbatim 00063 *> 00064 *> \param[in,out] AP 00065 *> \verbatim 00066 *> AP is REAL array, dimension (N*(N+1)/2) 00067 *> On entry, the upper or lower triangular matrix A, stored 00068 *> columnwise in a linear array. The j-th column of A is stored 00069 *> in the array AP as follows: 00070 *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 00071 *> if UPLO = 'L', AP(i + (j-1)*((2*n-j)/2) = A(i,j) for j<=i<=n. 00072 *> See below for further details. 00073 *> On exit, the (triangular) inverse of the original matrix, in 00074 *> the same packed storage format. 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, A(i,i) is exactly zero. The triangular 00083 *> matrix is singular and its inverse can not be computed. 00084 *> \endverbatim 00085 * 00086 * Authors: 00087 * ======== 00088 * 00089 *> \author Univ. of Tennessee 00090 *> \author Univ. of California Berkeley 00091 *> \author Univ. of Colorado Denver 00092 *> \author NAG Ltd. 00093 * 00094 *> \date November 2011 00095 * 00096 *> \ingroup realOTHERcomputational 00097 * 00098 *> \par Further Details: 00099 * ===================== 00100 *> 00101 *> \verbatim 00102 *> 00103 *> A triangular matrix A can be transferred to packed storage using one 00104 *> of the following program segments: 00105 *> 00106 *> UPLO = 'U': UPLO = 'L': 00107 *> 00108 *> JC = 1 JC = 1 00109 *> DO 2 J = 1, N DO 2 J = 1, N 00110 *> DO 1 I = 1, J DO 1 I = J, N 00111 *> AP(JC+I-1) = A(I,J) AP(JC+I-J) = A(I,J) 00112 *> 1 CONTINUE 1 CONTINUE 00113 *> JC = JC + J JC = JC + N - J + 1 00114 *> 2 CONTINUE 2 CONTINUE 00115 *> \endverbatim 00116 *> 00117 * ===================================================================== 00118 SUBROUTINE STPTRI( UPLO, DIAG, N, AP, INFO ) 00119 * 00120 * -- LAPACK computational routine (version 3.4.0) -- 00121 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00122 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00123 * November 2011 00124 * 00125 * .. Scalar Arguments .. 00126 CHARACTER DIAG, UPLO 00127 INTEGER INFO, N 00128 * .. 00129 * .. Array Arguments .. 00130 REAL AP( * ) 00131 * .. 00132 * 00133 * ===================================================================== 00134 * 00135 * .. Parameters .. 00136 REAL ONE, ZERO 00137 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00138 * .. 00139 * .. Local Scalars .. 00140 LOGICAL NOUNIT, UPPER 00141 INTEGER J, JC, JCLAST, JJ 00142 REAL AJJ 00143 * .. 00144 * .. External Functions .. 00145 LOGICAL LSAME 00146 EXTERNAL LSAME 00147 * .. 00148 * .. External Subroutines .. 00149 EXTERNAL SSCAL, STPMV, XERBLA 00150 * .. 00151 * .. Executable Statements .. 00152 * 00153 * Test the input parameters. 00154 * 00155 INFO = 0 00156 UPPER = LSAME( UPLO, 'U' ) 00157 NOUNIT = LSAME( DIAG, 'N' ) 00158 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00159 INFO = -1 00160 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN 00161 INFO = -2 00162 ELSE IF( N.LT.0 ) THEN 00163 INFO = -3 00164 END IF 00165 IF( INFO.NE.0 ) THEN 00166 CALL XERBLA( 'STPTRI', -INFO ) 00167 RETURN 00168 END IF 00169 * 00170 * Check for singularity if non-unit. 00171 * 00172 IF( NOUNIT ) THEN 00173 IF( UPPER ) THEN 00174 JJ = 0 00175 DO 10 INFO = 1, N 00176 JJ = JJ + INFO 00177 IF( AP( JJ ).EQ.ZERO ) 00178 $ RETURN 00179 10 CONTINUE 00180 ELSE 00181 JJ = 1 00182 DO 20 INFO = 1, N 00183 IF( AP( JJ ).EQ.ZERO ) 00184 $ RETURN 00185 JJ = JJ + N - INFO + 1 00186 20 CONTINUE 00187 END IF 00188 INFO = 0 00189 END IF 00190 * 00191 IF( UPPER ) THEN 00192 * 00193 * Compute inverse of upper triangular matrix. 00194 * 00195 JC = 1 00196 DO 30 J = 1, N 00197 IF( NOUNIT ) THEN 00198 AP( JC+J-1 ) = ONE / AP( JC+J-1 ) 00199 AJJ = -AP( JC+J-1 ) 00200 ELSE 00201 AJJ = -ONE 00202 END IF 00203 * 00204 * Compute elements 1:j-1 of j-th column. 00205 * 00206 CALL STPMV( 'Upper', 'No transpose', DIAG, J-1, AP, 00207 $ AP( JC ), 1 ) 00208 CALL SSCAL( J-1, AJJ, AP( JC ), 1 ) 00209 JC = JC + J 00210 30 CONTINUE 00211 * 00212 ELSE 00213 * 00214 * Compute inverse of lower triangular matrix. 00215 * 00216 JC = N*( N+1 ) / 2 00217 DO 40 J = N, 1, -1 00218 IF( NOUNIT ) THEN 00219 AP( JC ) = ONE / AP( JC ) 00220 AJJ = -AP( JC ) 00221 ELSE 00222 AJJ = -ONE 00223 END IF 00224 IF( J.LT.N ) THEN 00225 * 00226 * Compute elements j+1:n of j-th column. 00227 * 00228 CALL STPMV( 'Lower', 'No transpose', DIAG, N-J, 00229 $ AP( JCLAST ), AP( JC+1 ), 1 ) 00230 CALL SSCAL( N-J, AJJ, AP( JC+1 ), 1 ) 00231 END IF 00232 JCLAST = JC 00233 JC = JC - N + J - 2 00234 40 CONTINUE 00235 END IF 00236 * 00237 RETURN 00238 * 00239 * End of STPTRI 00240 * 00241 END