![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLANEG 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLANEG + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaneg.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaneg.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaneg.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * INTEGER FUNCTION SLANEG( N, D, LLD, SIGMA, PIVMIN, R ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER N, R 00025 * REAL PIVMIN, SIGMA 00026 * .. 00027 * .. Array Arguments .. 00028 * REAL D( * ), LLD( * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> SLANEG computes the Sturm count, the number of negative pivots 00038 *> encountered while factoring tridiagonal T - sigma I = L D L^T. 00039 *> This implementation works directly on the factors without forming 00040 *> the tridiagonal matrix T. The Sturm count is also the number of 00041 *> eigenvalues of T less than sigma. 00042 *> 00043 *> This routine is called from SLARRB. 00044 *> 00045 *> The current routine does not use the PIVMIN parameter but rather 00046 *> requires IEEE-754 propagation of Infinities and NaNs. This 00047 *> routine also has no input range restrictions but does require 00048 *> default exception handling such that x/0 produces Inf when x is 00049 *> non-zero, and Inf/Inf produces NaN. For more information, see: 00050 *> 00051 *> Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in 00052 *> Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on 00053 *> Scientific Computing, v28, n5, 2006. DOI 10.1137/050641624 00054 *> (Tech report version in LAWN 172 with the same title.) 00055 *> \endverbatim 00056 * 00057 * Arguments: 00058 * ========== 00059 * 00060 *> \param[in] N 00061 *> \verbatim 00062 *> N is INTEGER 00063 *> The order of the matrix. 00064 *> \endverbatim 00065 *> 00066 *> \param[in] D 00067 *> \verbatim 00068 *> D is REAL array, dimension (N) 00069 *> The N diagonal elements of the diagonal matrix D. 00070 *> \endverbatim 00071 *> 00072 *> \param[in] LLD 00073 *> \verbatim 00074 *> LLD is REAL array, dimension (N-1) 00075 *> The (N-1) elements L(i)*L(i)*D(i). 00076 *> \endverbatim 00077 *> 00078 *> \param[in] SIGMA 00079 *> \verbatim 00080 *> SIGMA is REAL 00081 *> Shift amount in T - sigma I = L D L^T. 00082 *> \endverbatim 00083 *> 00084 *> \param[in] PIVMIN 00085 *> \verbatim 00086 *> PIVMIN is REAL 00087 *> The minimum pivot in the Sturm sequence. May be used 00088 *> when zero pivots are encountered on non-IEEE-754 00089 *> architectures. 00090 *> \endverbatim 00091 *> 00092 *> \param[in] R 00093 *> \verbatim 00094 *> R is INTEGER 00095 *> The twist index for the twisted factorization that is used 00096 *> for the negcount. 00097 *> \endverbatim 00098 * 00099 * Authors: 00100 * ======== 00101 * 00102 *> \author Univ. of Tennessee 00103 *> \author Univ. of California Berkeley 00104 *> \author Univ. of Colorado Denver 00105 *> \author NAG Ltd. 00106 * 00107 *> \date November 2011 00108 * 00109 *> \ingroup auxOTHERauxiliary 00110 * 00111 *> \par Contributors: 00112 * ================== 00113 *> 00114 *> Osni Marques, LBNL/NERSC, USA \n 00115 *> Christof Voemel, University of California, Berkeley, USA \n 00116 *> Jason Riedy, University of California, Berkeley, USA \n 00117 *> 00118 * ===================================================================== 00119 INTEGER FUNCTION SLANEG( N, D, LLD, SIGMA, PIVMIN, R ) 00120 * 00121 * -- LAPACK auxiliary routine (version 3.4.0) -- 00122 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00123 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00124 * November 2011 00125 * 00126 * .. Scalar Arguments .. 00127 INTEGER N, R 00128 REAL PIVMIN, SIGMA 00129 * .. 00130 * .. Array Arguments .. 00131 REAL D( * ), LLD( * ) 00132 * .. 00133 * 00134 * ===================================================================== 00135 * 00136 * .. Parameters .. 00137 REAL ZERO, ONE 00138 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00139 * Some architectures propagate Infinities and NaNs very slowly, so 00140 * the code computes counts in BLKLEN chunks. Then a NaN can 00141 * propagate at most BLKLEN columns before being detected. This is 00142 * not a general tuning parameter; it needs only to be just large 00143 * enough that the overhead is tiny in common cases. 00144 INTEGER BLKLEN 00145 PARAMETER ( BLKLEN = 128 ) 00146 * .. 00147 * .. Local Scalars .. 00148 INTEGER BJ, J, NEG1, NEG2, NEGCNT 00149 REAL BSAV, DMINUS, DPLUS, GAMMA, P, T, TMP 00150 LOGICAL SAWNAN 00151 * .. 00152 * .. Intrinsic Functions .. 00153 INTRINSIC MIN, MAX 00154 * .. 00155 * .. External Functions .. 00156 LOGICAL SISNAN 00157 EXTERNAL SISNAN 00158 * .. 00159 * .. Executable Statements .. 00160 00161 NEGCNT = 0 00162 00163 * I) upper part: L D L^T - SIGMA I = L+ D+ L+^T 00164 T = -SIGMA 00165 DO 210 BJ = 1, R-1, BLKLEN 00166 NEG1 = 0 00167 BSAV = T 00168 DO 21 J = BJ, MIN(BJ+BLKLEN-1, R-1) 00169 DPLUS = D( J ) + T 00170 IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1 00171 TMP = T / DPLUS 00172 T = TMP * LLD( J ) - SIGMA 00173 21 CONTINUE 00174 SAWNAN = SISNAN( T ) 00175 * Run a slower version of the above loop if a NaN is detected. 00176 * A NaN should occur only with a zero pivot after an infinite 00177 * pivot. In that case, substituting 1 for T/DPLUS is the 00178 * correct limit. 00179 IF( SAWNAN ) THEN 00180 NEG1 = 0 00181 T = BSAV 00182 DO 22 J = BJ, MIN(BJ+BLKLEN-1, R-1) 00183 DPLUS = D( J ) + T 00184 IF( DPLUS.LT.ZERO ) NEG1 = NEG1 + 1 00185 TMP = T / DPLUS 00186 IF (SISNAN(TMP)) TMP = ONE 00187 T = TMP * LLD(J) - SIGMA 00188 22 CONTINUE 00189 END IF 00190 NEGCNT = NEGCNT + NEG1 00191 210 CONTINUE 00192 * 00193 * II) lower part: L D L^T - SIGMA I = U- D- U-^T 00194 P = D( N ) - SIGMA 00195 DO 230 BJ = N-1, R, -BLKLEN 00196 NEG2 = 0 00197 BSAV = P 00198 DO 23 J = BJ, MAX(BJ-BLKLEN+1, R), -1 00199 DMINUS = LLD( J ) + P 00200 IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1 00201 TMP = P / DMINUS 00202 P = TMP * D( J ) - SIGMA 00203 23 CONTINUE 00204 SAWNAN = SISNAN( P ) 00205 * As above, run a slower version that substitutes 1 for Inf/Inf. 00206 * 00207 IF( SAWNAN ) THEN 00208 NEG2 = 0 00209 P = BSAV 00210 DO 24 J = BJ, MAX(BJ-BLKLEN+1, R), -1 00211 DMINUS = LLD( J ) + P 00212 IF( DMINUS.LT.ZERO ) NEG2 = NEG2 + 1 00213 TMP = P / DMINUS 00214 IF (SISNAN(TMP)) TMP = ONE 00215 P = TMP * D(J) - SIGMA 00216 24 CONTINUE 00217 END IF 00218 NEGCNT = NEGCNT + NEG2 00219 230 CONTINUE 00220 * 00221 * III) Twist index 00222 * T was shifted by SIGMA initially. 00223 GAMMA = (T + SIGMA) + P 00224 IF( GAMMA.LT.ZERO ) NEGCNT = NEGCNT+1 00225 00226 SLANEG = NEGCNT 00227 END