00001 /* Ergo, version 3.7, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2018 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, 00004 * and Anastasia Kruchinina. 00005 * 00006 * This program is free software: you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation, either version 3 of the License, or 00009 * (at your option) any later version. 00010 * 00011 * This program is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00018 * 00019 * Primary academic reference: 00020 * Ergo: An open-source program for linear-scaling electronic structure 00021 * calculations, 00022 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia 00023 * Kruchinina, 00024 * SoftwareX 7, 107 (2018), 00025 * <http://dx.doi.org/10.1016/j.softx.2018.03.005> 00026 * 00027 * For further information about Ergo, see <http://www.ergoscf.org>. 00028 */ 00029 00030 /* This file belongs to the template_lapack part of the Ergo source 00031 * code. The source files in the template_lapack directory are modified 00032 * versions of files originally distributed as CLAPACK, see the 00033 * Copyright/license notice in the file template_lapack/COPYING. 00034 */ 00035 00036 00037 #ifndef TEMPLATE_LAPACK_LARRB_HEADER 00038 #define TEMPLATE_LAPACK_LARRB_HEADER 00039 00040 00041 #include "template_lapack_laneg.h" 00042 00043 00044 template<class Treal> 00045 int template_lapack_larrb(integer *n, Treal *d__, Treal *lld, 00046 integer *ifirst, integer *ilast, Treal *rtol1, Treal *rtol2, 00047 integer *offset, Treal *w, Treal *wgap, Treal *werr, 00048 Treal *work, integer *iwork, Treal *pivmin, Treal * 00049 spdiam, integer *twist, integer *info) 00050 { 00051 /* System generated locals */ 00052 integer i__1; 00053 Treal d__1, d__2; 00054 00055 /* Local variables */ 00056 integer i__, k, r__, i1, ii, ip; 00057 Treal gap, mid, tmp, back, lgap, rgap, left; 00058 integer iter, nint, prev, next; 00059 Treal cvrgd, right, width; 00060 integer negcnt; 00061 Treal mnwdth; 00062 integer olnint, maxitr; 00063 00064 00065 /* -- LAPACK auxiliary routine (version 3.2) -- */ 00066 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00067 /* November 2006 */ 00068 00069 /* .. Scalar Arguments .. */ 00070 /* .. */ 00071 /* .. Array Arguments .. */ 00072 /* .. */ 00073 00074 /* Purpose */ 00075 /* ======= */ 00076 00077 /* Given the relatively robust representation(RRR) L D L^T, DLARRB */ 00078 /* does "limited" bisection to refine the eigenvalues of L D L^T, */ 00079 /* W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial */ 00080 /* guesses for these eigenvalues are input in W, the corresponding estimate */ 00081 /* of the error in these guesses and their gaps are input in WERR */ 00082 /* and WGAP, respectively. During bisection, intervals */ 00083 /* [left, right] are maintained by storing their mid-points and */ 00084 /* semi-widths in the arrays W and WERR respectively. */ 00085 00086 /* Arguments */ 00087 /* ========= */ 00088 00089 /* N (input) INTEGER */ 00090 /* The order of the matrix. */ 00091 00092 /* D (input) DOUBLE PRECISION array, dimension (N) */ 00093 /* The N diagonal elements of the diagonal matrix D. */ 00094 00095 /* LLD (input) DOUBLE PRECISION array, dimension (N-1) */ 00096 /* The (N-1) elements L(i)*L(i)*D(i). */ 00097 00098 /* IFIRST (input) INTEGER */ 00099 /* The index of the first eigenvalue to be computed. */ 00100 00101 /* ILAST (input) INTEGER */ 00102 /* The index of the last eigenvalue to be computed. */ 00103 00104 /* RTOL1 (input) DOUBLE PRECISION */ 00105 /* RTOL2 (input) DOUBLE PRECISION */ 00106 /* Tolerance for the convergence of the bisection intervals. */ 00107 /* An interval [LEFT,RIGHT] has converged if */ 00108 /* RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) */ 00109 /* where GAP is the (estimated) distance to the nearest */ 00110 /* eigenvalue. */ 00111 00112 /* OFFSET (input) INTEGER */ 00113 /* Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET */ 00114 /* through ILAST-OFFSET elements of these arrays are to be used. */ 00115 00116 /* W (input/output) DOUBLE PRECISION array, dimension (N) */ 00117 /* On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are */ 00118 /* estimates of the eigenvalues of L D L^T indexed IFIRST throug */ 00119 /* ILAST. */ 00120 /* On output, these estimates are refined. */ 00121 00122 /* WGAP (input/output) DOUBLE PRECISION array, dimension (N-1) */ 00123 /* On input, the (estimated) gaps between consecutive */ 00124 /* eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between */ 00125 /* eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST */ 00126 /* then WGAP(IFIRST-OFFSET) must be set to ZERO. */ 00127 /* On output, these gaps are refined. */ 00128 00129 /* WERR (input/output) DOUBLE PRECISION array, dimension (N) */ 00130 /* On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are */ 00131 /* the errors in the estimates of the corresponding elements in W. */ 00132 /* On output, these errors are refined. */ 00133 00134 /* WORK (workspace) DOUBLE PRECISION array, dimension (2*N) */ 00135 /* Workspace. */ 00136 00137 /* IWORK (workspace) INTEGER array, dimension (2*N) */ 00138 /* Workspace. */ 00139 00140 /* PIVMIN (input) DOUBLE PRECISION */ 00141 /* The minimum pivot in the Sturm sequence. */ 00142 00143 /* SPDIAM (input) DOUBLE PRECISION */ 00144 /* The spectral diameter of the matrix. */ 00145 00146 /* TWIST (input) INTEGER */ 00147 /* The twist index for the twisted factorization that is used */ 00148 /* for the negcount. */ 00149 /* TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T */ 00150 /* TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T */ 00151 /* TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r) */ 00152 00153 /* INFO (output) INTEGER */ 00154 /* Error flag. */ 00155 00156 /* Further Details */ 00157 /* =============== */ 00158 00159 /* Based on contributions by */ 00160 /* Beresford Parlett, University of California, Berkeley, USA */ 00161 /* Jim Demmel, University of California, Berkeley, USA */ 00162 /* Inderjit Dhillon, University of Texas, Austin, USA */ 00163 /* Osni Marques, LBNL/NERSC, USA */ 00164 /* Christof Voemel, University of California, Berkeley, USA */ 00165 00166 /* ===================================================================== */ 00167 00168 /* .. Parameters .. */ 00169 /* .. */ 00170 /* .. Local Scalars .. */ 00171 /* .. */ 00172 /* .. External Functions .. */ 00173 00174 /* .. */ 00175 /* .. Intrinsic Functions .. */ 00176 /* .. */ 00177 /* .. Executable Statements .. */ 00178 00179 /* Parameter adjustments */ 00180 --iwork; 00181 --work; 00182 --werr; 00183 --wgap; 00184 --w; 00185 --lld; 00186 --d__; 00187 00188 /* Function Body */ 00189 *info = 0; 00190 00191 maxitr = (integer) ((template_blas_log(*spdiam + *pivmin) - template_blas_log(*pivmin)) / template_blas_log(2.)) + 00192 2; 00193 mnwdth = *pivmin * 2.; 00194 00195 r__ = *twist; 00196 if (r__ < 1 || r__ > *n) { 00197 r__ = *n; 00198 } 00199 00200 /* Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ]. */ 00201 /* The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while */ 00202 /* Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 ) */ 00203 /* for an unconverged interval is set to the index of the next unconverged */ 00204 /* interval, and is -1 or 0 for a converged interval. Thus a linked */ 00205 /* list of unconverged intervals is set up. */ 00206 00207 i1 = *ifirst; 00208 /* The number of unconverged intervals */ 00209 nint = 0; 00210 /* The last unconverged interval found */ 00211 prev = 0; 00212 rgap = wgap[i1 - *offset]; 00213 i__1 = *ilast; 00214 for (i__ = i1; i__ <= i__1; ++i__) { 00215 k = i__ << 1; 00216 ii = i__ - *offset; 00217 left = w[ii] - werr[ii]; 00218 right = w[ii] + werr[ii]; 00219 lgap = rgap; 00220 rgap = wgap[ii]; 00221 gap = minMACRO(lgap,rgap); 00222 /* Make sure that [LEFT,RIGHT] contains the desired eigenvalue */ 00223 /* Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT */ 00224 00225 /* Do while( NEGCNT(LEFT).GT.I-1 ) */ 00226 00227 back = werr[ii]; 00228 L20: 00229 negcnt = template_lapack_laneg(n, &d__[1], &lld[1], &left, pivmin, &r__); 00230 if (negcnt > i__ - 1) { 00231 left -= back; 00232 back *= 2.; 00233 goto L20; 00234 } 00235 00236 /* Do while( NEGCNT(RIGHT).LT.I ) */ 00237 /* Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT */ 00238 00239 back = werr[ii]; 00240 L50: 00241 negcnt = template_lapack_laneg(n, &d__[1], &lld[1], &right, pivmin, &r__); 00242 if (negcnt < i__) { 00243 right += back; 00244 back *= 2.; 00245 goto L50; 00246 } 00247 width = (d__1 = left - right, absMACRO(d__1)) * .5; 00248 /* Computing MAX */ 00249 d__1 = absMACRO(left), d__2 = absMACRO(right); 00250 tmp = maxMACRO(d__1,d__2); 00251 /* Computing MAX */ 00252 d__1 = *rtol1 * gap, d__2 = *rtol2 * tmp; 00253 cvrgd = maxMACRO(d__1,d__2); 00254 if (width <= cvrgd || width <= mnwdth) { 00255 /* This interval has already converged and does not need refinement. */ 00256 /* (Note that the gaps might change through refining the */ 00257 /* eigenvalues, however, they can only get bigger.) */ 00258 /* Remove it from the list. */ 00259 iwork[k - 1] = -1; 00260 /* Make sure that I1 always points to the first unconverged interval */ 00261 if (i__ == i1 && i__ < *ilast) { 00262 i1 = i__ + 1; 00263 } 00264 if (prev >= i1 && i__ <= *ilast) { 00265 iwork[(prev << 1) - 1] = i__ + 1; 00266 } 00267 } else { 00268 /* unconverged interval found */ 00269 prev = i__; 00270 ++nint; 00271 iwork[k - 1] = i__ + 1; 00272 iwork[k] = negcnt; 00273 } 00274 work[k - 1] = left; 00275 work[k] = right; 00276 /* L75: */ 00277 } 00278 00279 /* Do while( NINT.GT.0 ), i.e. there are still unconverged intervals */ 00280 /* and while (ITER.LT.MAXITR) */ 00281 00282 iter = 0; 00283 L80: 00284 prev = i1 - 1; 00285 i__ = i1; 00286 olnint = nint; 00287 i__1 = olnint; 00288 for (ip = 1; ip <= i__1; ++ip) { 00289 k = i__ << 1; 00290 ii = i__ - *offset; 00291 rgap = wgap[ii]; 00292 lgap = rgap; 00293 if (ii > 1) { 00294 lgap = wgap[ii - 1]; 00295 } 00296 gap = minMACRO(lgap,rgap); 00297 next = iwork[k - 1]; 00298 left = work[k - 1]; 00299 right = work[k]; 00300 mid = (left + right) * .5; 00301 /* semiwidth of interval */ 00302 width = right - mid; 00303 /* Computing MAX */ 00304 d__1 = absMACRO(left), d__2 = absMACRO(right); 00305 tmp = maxMACRO(d__1,d__2); 00306 /* Computing MAX */ 00307 d__1 = *rtol1 * gap, d__2 = *rtol2 * tmp; 00308 cvrgd = maxMACRO(d__1,d__2); 00309 if (width <= cvrgd || width <= mnwdth || iter == maxitr) { 00310 /* reduce number of unconverged intervals */ 00311 --nint; 00312 /* Mark interval as converged. */ 00313 iwork[k - 1] = 0; 00314 if (i1 == i__) { 00315 i1 = next; 00316 } else { 00317 /* Prev holds the last unconverged interval previously examined */ 00318 if (prev >= i1) { 00319 iwork[(prev << 1) - 1] = next; 00320 } 00321 } 00322 i__ = next; 00323 goto L100; 00324 } 00325 prev = i__; 00326 00327 /* Perform one bisection step */ 00328 00329 negcnt = template_lapack_laneg(n, &d__[1], &lld[1], &mid, pivmin, &r__); 00330 if (negcnt <= i__ - 1) { 00331 work[k - 1] = mid; 00332 } else { 00333 work[k] = mid; 00334 } 00335 i__ = next; 00336 L100: 00337 ; 00338 } 00339 ++iter; 00340 /* do another loop if there are still unconverged intervals */ 00341 /* However, in the last iteration, all intervals are accepted */ 00342 /* since this is the best we can do. */ 00343 if (nint > 0 && iter <= maxitr) { 00344 goto L80; 00345 } 00346 00347 00348 /* At this point, all the intervals have converged */ 00349 i__1 = *ilast; 00350 for (i__ = *ifirst; i__ <= i__1; ++i__) { 00351 k = i__ << 1; 00352 ii = i__ - *offset; 00353 /* All intervals marked by '0' have been refined. */ 00354 if (iwork[k - 1] == 0) { 00355 w[ii] = (work[k - 1] + work[k]) * .5; 00356 werr[ii] = work[k] - w[ii]; 00357 } 00358 /* L110: */ 00359 } 00360 00361 i__1 = *ilast; 00362 for (i__ = *ifirst + 1; i__ <= i__1; ++i__) { 00363 k = i__ << 1; 00364 ii = i__ - *offset; 00365 /* Computing MAX */ 00366 d__1 = 0., d__2 = w[ii] - werr[ii] - w[ii - 1] - werr[ii - 1]; 00367 wgap[ii - 1] = maxMACRO(d__1,d__2); 00368 /* L111: */ 00369 } 00370 return 0; 00371 00372 /* End of DLARRB */ 00373 00374 } /* dlarrb_ */ 00375 00376 #endif