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_LAGTF_HEADER 00036 #define TEMPLATE_LAPACK_LAGTF_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_lagtf(const integer *n, Treal *a, const Treal *lambda, 00041 Treal *b, Treal *c__, const Treal *tol, Treal *d__, 00042 integer *in, integer *info) 00043 { 00044 /* -- LAPACK routine (version 3.0) -- 00045 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00046 Courant Institute, Argonne National Lab, and Rice University 00047 June 30, 1999 00048 00049 00050 Purpose 00051 ======= 00052 00053 DLAGTF factorizes the matrix (T - lambda*I), where T is an n by n 00054 tridiagonal matrix and lambda is a scalar, as 00055 00056 T - lambda*I = PLU, 00057 00058 where P is a permutation matrix, L is a unit lower tridiagonal matrix 00059 with at most one non-zero sub-diagonal elements per column and U is 00060 an upper triangular matrix with at most two non-zero super-diagonal 00061 elements per column. 00062 00063 The factorization is obtained by Gaussian elimination with partial 00064 pivoting and implicit row scaling. 00065 00066 The parameter LAMBDA is included in the routine so that DLAGTF may 00067 be used, in conjunction with DLAGTS, to obtain eigenvectors of T by 00068 inverse iteration. 00069 00070 Arguments 00071 ========= 00072 00073 N (input) INTEGER 00074 The order of the matrix T. 00075 00076 A (input/output) DOUBLE PRECISION array, dimension (N) 00077 On entry, A must contain the diagonal elements of T. 00078 00079 On exit, A is overwritten by the n diagonal elements of the 00080 upper triangular matrix U of the factorization of T. 00081 00082 LAMBDA (input) DOUBLE PRECISION 00083 On entry, the scalar lambda. 00084 00085 B (input/output) DOUBLE PRECISION array, dimension (N-1) 00086 On entry, B must contain the (n-1) super-diagonal elements of 00087 T. 00088 00089 On exit, B is overwritten by the (n-1) super-diagonal 00090 elements of the matrix U of the factorization of T. 00091 00092 C (input/output) DOUBLE PRECISION array, dimension (N-1) 00093 On entry, C must contain the (n-1) sub-diagonal elements of 00094 T. 00095 00096 On exit, C is overwritten by the (n-1) sub-diagonal elements 00097 of the matrix L of the factorization of T. 00098 00099 TOL (input) DOUBLE PRECISION 00100 On entry, a relative tolerance used to indicate whether or 00101 not the matrix (T - lambda*I) is nearly singular. TOL should 00102 normally be chose as approximately the largest relative error 00103 in the elements of T. For example, if the elements of T are 00104 correct to about 4 significant figures, then TOL should be 00105 set to about 5*10**(-4). If TOL is supplied as less than eps, 00106 where eps is the relative machine precision, then the value 00107 eps is used in place of TOL. 00108 00109 D (output) DOUBLE PRECISION array, dimension (N-2) 00110 On exit, D is overwritten by the (n-2) second super-diagonal 00111 elements of the matrix U of the factorization of T. 00112 00113 IN (output) INTEGER array, dimension (N) 00114 On exit, IN contains details of the permutation matrix P. If 00115 an interchange occurred at the kth step of the elimination, 00116 then IN(k) = 1, otherwise IN(k) = 0. The element IN(n) 00117 returns the smallest positive integer j such that 00118 00119 abs( u(j,j) ).le. norm( (T - lambda*I)(j) )*TOL, 00120 00121 where norm( A(j) ) denotes the sum of the absolute values of 00122 the jth row of the matrix A. If no such j exists then IN(n) 00123 is returned as zero. If IN(n) is returned as positive, then a 00124 diagonal element of U is small, indicating that 00125 (T - lambda*I) is singular or nearly singular, 00126 00127 INFO (output) INTEGER 00128 = 0 : successful exit 00129 .lt. 0: if INFO = -k, the kth argument had an illegal value 00130 00131 ===================================================================== 00132 00133 00134 Parameter adjustments */ 00135 /* System generated locals */ 00136 integer i__1; 00137 Treal d__1, d__2; 00138 /* Local variables */ 00139 Treal temp, mult; 00140 integer k; 00141 Treal scale1, scale2; 00142 Treal tl; 00143 Treal eps, piv1, piv2; 00144 00145 --in; 00146 --d__; 00147 --c__; 00148 --b; 00149 --a; 00150 00151 /* Function Body */ 00152 *info = 0; 00153 if (*n < 0) { 00154 *info = -1; 00155 i__1 = -(*info); 00156 template_blas_erbla("LAGTF ", &i__1); 00157 return 0; 00158 } 00159 00160 if (*n == 0) { 00161 return 0; 00162 } 00163 00164 a[1] -= *lambda; 00165 in[*n] = 0; 00166 if (*n == 1) { 00167 if (a[1] == 0.) { 00168 in[1] = 1; 00169 } 00170 return 0; 00171 } 00172 00173 eps = template_lapack_lamch("Epsilon", (Treal)0); 00174 00175 tl = maxMACRO(*tol,eps); 00176 scale1 = absMACRO(a[1]) + absMACRO(b[1]); 00177 i__1 = *n - 1; 00178 for (k = 1; k <= i__1; ++k) { 00179 a[k + 1] -= *lambda; 00180 scale2 = (d__1 = c__[k], absMACRO(d__1)) + (d__2 = a[k + 1], absMACRO(d__2)); 00181 if (k < *n - 1) { 00182 scale2 += (d__1 = b[k + 1], absMACRO(d__1)); 00183 } 00184 if (a[k] == 0.) { 00185 piv1 = 0.; 00186 } else { 00187 piv1 = (d__1 = a[k], absMACRO(d__1)) / scale1; 00188 } 00189 if (c__[k] == 0.) { 00190 in[k] = 0; 00191 piv2 = 0.; 00192 scale1 = scale2; 00193 if (k < *n - 1) { 00194 d__[k] = 0.; 00195 } 00196 } else { 00197 piv2 = (d__1 = c__[k], absMACRO(d__1)) / scale2; 00198 if (piv2 <= piv1) { 00199 in[k] = 0; 00200 scale1 = scale2; 00201 c__[k] /= a[k]; 00202 a[k + 1] -= c__[k] * b[k]; 00203 if (k < *n - 1) { 00204 d__[k] = 0.; 00205 } 00206 } else { 00207 in[k] = 1; 00208 mult = a[k] / c__[k]; 00209 a[k] = c__[k]; 00210 temp = a[k + 1]; 00211 a[k + 1] = b[k] - mult * temp; 00212 if (k < *n - 1) { 00213 d__[k] = b[k + 1]; 00214 b[k + 1] = -mult * d__[k]; 00215 } 00216 b[k] = temp; 00217 c__[k] = mult; 00218 } 00219 } 00220 if (maxMACRO(piv1,piv2) <= tl && in[*n] == 0) { 00221 in[*n] = k; 00222 } 00223 /* L10: */ 00224 } 00225 if ((d__1 = a[*n], absMACRO(d__1)) <= scale1 * tl && in[*n] == 0) { 00226 in[*n] = *n; 00227 } 00228 00229 return 0; 00230 00231 /* End of DLAGTF */ 00232 00233 } /* dlagtf_ */ 00234 00235 #endif