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