ergo
|
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek. 00004 * 00005 * This program is free software: you can redistribute it and/or modify 00006 * it under the terms of the GNU General Public License as published by 00007 * the Free Software Foundation, either version 3 of the License, or 00008 * (at your option) any later version. 00009 * 00010 * This program is distributed in the hope that it will be useful, 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00013 * GNU General Public License for more details. 00014 * 00015 * You should have received a copy of the GNU General Public License 00016 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00017 * 00018 * Primary academic reference: 00019 * KohnâSham Density Functional Theory Electronic Structure Calculations 00020 * with Linearly Scaling Computational Time and Memory Usage, 00021 * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek, 00022 * J. Chem. Theory Comput. 7, 340 (2011), 00023 * <http://dx.doi.org/10.1021/ct100611z> 00024 * 00025 * For further information about Ergo, see <http://www.ergoscf.org>. 00026 */ 00027 00028 /* This file belongs to the template_lapack part of the Ergo source 00029 * code. The source files in the template_lapack directory are modified 00030 * versions of files originally distributed as CLAPACK, see the 00031 * Copyright/license notice in the file template_lapack/COPYING. 00032 */ 00033 00034 00035 #ifndef TEMPLATE_LAPACK_LAGTS_HEADER 00036 #define TEMPLATE_LAPACK_LAGTS_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_lagts(const integer *job, const integer *n, const Treal *a, 00041 const Treal *b, const Treal *c__, const Treal *d__, const integer *in, 00042 Treal *y, Treal *tol, integer *info) 00043 { 00044 /* -- LAPACK auxiliary routine (version 3.0) -- 00045 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00046 Courant Institute, Argonne National Lab, and Rice University 00047 October 31, 1992 00048 00049 00050 Purpose 00051 ======= 00052 00053 DLAGTS may be used to solve one of the systems of equations 00054 00055 (T - lambda*I)*x = y or (T - lambda*I)'*x = y, 00056 00057 where T is an n by n tridiagonal matrix, for x, following the 00058 factorization of (T - lambda*I) as 00059 00060 (T - lambda*I) = P*L*U , 00061 00062 by routine DLAGTF. The choice of equation to be solved is 00063 controlled by the argument JOB, and in each case there is an option 00064 to perturb zero or very small diagonal elements of U, this option 00065 being intended for use in applications such as inverse iteration. 00066 00067 Arguments 00068 ========= 00069 00070 JOB (input) INTEGER 00071 Specifies the job to be performed by DLAGTS as follows: 00072 = 1: The equations (T - lambda*I)x = y are to be solved, 00073 but diagonal elements of U are not to be perturbed. 00074 = -1: The equations (T - lambda*I)x = y are to be solved 00075 and, if overflow would otherwise occur, the diagonal 00076 elements of U are to be perturbed. See argument TOL 00077 below. 00078 = 2: The equations (T - lambda*I)'x = y are to be solved, 00079 but diagonal elements of U are not to be perturbed. 00080 = -2: The equations (T - lambda*I)'x = y are to be solved 00081 and, if overflow would otherwise occur, the diagonal 00082 elements of U are to be perturbed. See argument TOL 00083 below. 00084 00085 N (input) INTEGER 00086 The order of the matrix T. 00087 00088 A (input) DOUBLE PRECISION array, dimension (N) 00089 On entry, A must contain the diagonal elements of U as 00090 returned from DLAGTF. 00091 00092 B (input) DOUBLE PRECISION array, dimension (N-1) 00093 On entry, B must contain the first super-diagonal elements of 00094 U as returned from DLAGTF. 00095 00096 C (input) DOUBLE PRECISION array, dimension (N-1) 00097 On entry, C must contain the sub-diagonal elements of L as 00098 returned from DLAGTF. 00099 00100 D (input) DOUBLE PRECISION array, dimension (N-2) 00101 On entry, D must contain the second super-diagonal elements 00102 of U as returned from DLAGTF. 00103 00104 IN (input) INTEGER array, dimension (N) 00105 On entry, IN must contain details of the matrix P as returned 00106 from DLAGTF. 00107 00108 Y (input/output) DOUBLE PRECISION array, dimension (N) 00109 On entry, the right hand side vector y. 00110 On exit, Y is overwritten by the solution vector x. 00111 00112 TOL (input/output) DOUBLE PRECISION 00113 On entry, with JOB .lt. 0, TOL should be the minimum 00114 perturbation to be made to very small diagonal elements of U. 00115 TOL should normally be chosen as about eps*norm(U), where eps 00116 is the relative machine precision, but if TOL is supplied as 00117 non-positive, then it is reset to eps*max( abs( u(i,j) ) ). 00118 If JOB .gt. 0 then TOL is not referenced. 00119 00120 On exit, TOL is changed as described above, only if TOL is 00121 non-positive on entry. Otherwise TOL is unchanged. 00122 00123 INFO (output) INTEGER 00124 = 0 : successful exit 00125 .lt. 0: if INFO = -i, the i-th argument had an illegal value 00126 .gt. 0: overflow would occur when computing the INFO(th) 00127 element of the solution vector x. This can only occur 00128 when JOB is supplied as positive and either means 00129 that a diagonal element of U is very small, or that 00130 the elements of the right-hand side vector y are very 00131 large. 00132 00133 ===================================================================== 00134 00135 00136 Parameter adjustments */ 00137 /* System generated locals */ 00138 integer i__1; 00139 Treal d__1, d__2, d__3, d__4, d__5; 00140 /* Local variables */ 00141 Treal temp, pert; 00142 integer k; 00143 Treal absak, sfmin, ak; 00144 Treal bignum, eps; 00145 00146 --y; 00147 --in; 00148 --d__; 00149 --c__; 00150 --b; 00151 --a; 00152 00153 /* Function Body */ 00154 *info = 0; 00155 if (absMACRO(*job) > 2 || *job == 0) { 00156 *info = -1; 00157 } else if (*n < 0) { 00158 *info = -2; 00159 } 00160 if (*info != 0) { 00161 i__1 = -(*info); 00162 template_blas_erbla("LAGTS ", &i__1); 00163 return 0; 00164 } 00165 00166 if (*n == 0) { 00167 return 0; 00168 } 00169 00170 eps = template_lapack_lamch("Epsilon", (Treal)0); 00171 sfmin = template_lapack_lamch("Safe minimum", (Treal)0); 00172 bignum = 1. / sfmin; 00173 00174 if (*job < 0) { 00175 if (*tol <= 0.) { 00176 *tol = absMACRO(a[1]); 00177 if (*n > 1) { 00178 /* Computing MAX */ 00179 d__1 = *tol, d__2 = absMACRO(a[2]), d__1 = maxMACRO(d__1,d__2), d__2 = 00180 absMACRO(b[1]); 00181 *tol = maxMACRO(d__1,d__2); 00182 } 00183 i__1 = *n; 00184 for (k = 3; k <= i__1; ++k) { 00185 /* Computing MAX */ 00186 d__4 = *tol, d__5 = (d__1 = a[k], absMACRO(d__1)), d__4 = maxMACRO(d__4, 00187 d__5), d__5 = (d__2 = b[k - 1], absMACRO(d__2)), d__4 = 00188 maxMACRO(d__4,d__5), d__5 = (d__3 = d__[k - 2], absMACRO(d__3)); 00189 *tol = maxMACRO(d__4,d__5); 00190 /* L10: */ 00191 } 00192 *tol *= eps; 00193 if (*tol == 0.) { 00194 *tol = eps; 00195 } 00196 } 00197 } 00198 00199 if (absMACRO(*job) == 1) { 00200 i__1 = *n; 00201 for (k = 2; k <= i__1; ++k) { 00202 if (in[k - 1] == 0) { 00203 y[k] -= c__[k - 1] * y[k - 1]; 00204 } else { 00205 temp = y[k - 1]; 00206 y[k - 1] = y[k]; 00207 y[k] = temp - c__[k - 1] * y[k]; 00208 } 00209 /* L20: */ 00210 } 00211 if (*job == 1) { 00212 for (k = *n; k >= 1; --k) { 00213 if (k <= *n - 2) { 00214 temp = y[k] - b[k] * y[k + 1] - d__[k] * y[k + 2]; 00215 } else if (k == *n - 1) { 00216 temp = y[k] - b[k] * y[k + 1]; 00217 } else { 00218 temp = y[k]; 00219 } 00220 ak = a[k]; 00221 absak = absMACRO(ak); 00222 if (absak < 1.) { 00223 if (absak < sfmin) { 00224 if (absak == 0. || absMACRO(temp) * sfmin > absak) { 00225 *info = k; 00226 return 0; 00227 } else { 00228 temp *= bignum; 00229 ak *= bignum; 00230 } 00231 } else if (absMACRO(temp) > absak * bignum) { 00232 *info = k; 00233 return 0; 00234 } 00235 } 00236 y[k] = temp / ak; 00237 /* L30: */ 00238 } 00239 } else { 00240 for (k = *n; k >= 1; --k) { 00241 if (k <= *n - 2) { 00242 temp = y[k] - b[k] * y[k + 1] - d__[k] * y[k + 2]; 00243 } else if (k == *n - 1) { 00244 temp = y[k] - b[k] * y[k + 1]; 00245 } else { 00246 temp = y[k]; 00247 } 00248 ak = a[k]; 00249 pert = template_lapack_d_sign(tol, &ak); 00250 L40: 00251 absak = absMACRO(ak); 00252 if (absak < 1.) { 00253 if (absak < sfmin) { 00254 if (absak == 0. || absMACRO(temp) * sfmin > absak) { 00255 ak += pert; 00256 pert *= 2; 00257 goto L40; 00258 } else { 00259 temp *= bignum; 00260 ak *= bignum; 00261 } 00262 } else if (absMACRO(temp) > absak * bignum) { 00263 ak += pert; 00264 pert *= 2; 00265 goto L40; 00266 } 00267 } 00268 y[k] = temp / ak; 00269 /* L50: */ 00270 } 00271 } 00272 } else { 00273 00274 /* Come to here if JOB = 2 or -2 */ 00275 00276 if (*job == 2) { 00277 i__1 = *n; 00278 for (k = 1; k <= i__1; ++k) { 00279 if (k >= 3) { 00280 temp = y[k] - b[k - 1] * y[k - 1] - d__[k - 2] * y[k - 2]; 00281 } else if (k == 2) { 00282 temp = y[k] - b[k - 1] * y[k - 1]; 00283 } else { 00284 temp = y[k]; 00285 } 00286 ak = a[k]; 00287 absak = absMACRO(ak); 00288 if (absak < 1.) { 00289 if (absak < sfmin) { 00290 if (absak == 0. || absMACRO(temp) * sfmin > absak) { 00291 *info = k; 00292 return 0; 00293 } else { 00294 temp *= bignum; 00295 ak *= bignum; 00296 } 00297 } else if (absMACRO(temp) > absak * bignum) { 00298 *info = k; 00299 return 0; 00300 } 00301 } 00302 y[k] = temp / ak; 00303 /* L60: */ 00304 } 00305 } else { 00306 i__1 = *n; 00307 for (k = 1; k <= i__1; ++k) { 00308 if (k >= 3) { 00309 temp = y[k] - b[k - 1] * y[k - 1] - d__[k - 2] * y[k - 2]; 00310 } else if (k == 2) { 00311 temp = y[k] - b[k - 1] * y[k - 1]; 00312 } else { 00313 temp = y[k]; 00314 } 00315 ak = a[k]; 00316 pert = template_lapack_d_sign(tol, &ak); 00317 L70: 00318 absak = absMACRO(ak); 00319 if (absak < 1.) { 00320 if (absak < sfmin) { 00321 if (absak == 0. || absMACRO(temp) * sfmin > absak) { 00322 ak += pert; 00323 pert *= 2; 00324 goto L70; 00325 } else { 00326 temp *= bignum; 00327 ak *= bignum; 00328 } 00329 } else if (absMACRO(temp) > absak * bignum) { 00330 ak += pert; 00331 pert *= 2; 00332 goto L70; 00333 } 00334 } 00335 y[k] = temp / ak; 00336 /* L80: */ 00337 } 00338 } 00339 00340 for (k = *n; k >= 2; --k) { 00341 if (in[k - 1] == 0) { 00342 y[k - 1] -= c__[k - 1] * y[k]; 00343 } else { 00344 temp = y[k - 1]; 00345 y[k - 1] = y[k]; 00346 y[k] = temp - c__[k - 1] * y[k]; 00347 } 00348 /* L90: */ 00349 } 00350 } 00351 00352 /* End of DLAGTS */ 00353 00354 return 0; 00355 } /* dlagts_ */ 00356 00357 #endif