ergo
template_lapack_lagtf.h
Go to the documentation of this file.
1 /* Ergo, version 3.3, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2013 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28  /* This file belongs to the template_lapack part of the Ergo source
29  * code. The source files in the template_lapack directory are modified
30  * versions of files originally distributed as CLAPACK, see the
31  * Copyright/license notice in the file template_lapack/COPYING.
32  */
33 
34 
35 #ifndef TEMPLATE_LAPACK_LAGTF_HEADER
36 #define TEMPLATE_LAPACK_LAGTF_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_lagtf(const integer *n, Treal *a, const Treal *lambda,
41  Treal *b, Treal *c__, const Treal *tol, Treal *d__,
42  integer *in, integer *info)
43 {
44 /* -- LAPACK routine (version 3.0) --
45  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
46  Courant Institute, Argonne National Lab, and Rice University
47  June 30, 1999
48 
49 
50  Purpose
51  =======
52 
53  DLAGTF factorizes the matrix (T - lambda*I), where T is an n by n
54  tridiagonal matrix and lambda is a scalar, as
55 
56  T - lambda*I = PLU,
57 
58  where P is a permutation matrix, L is a unit lower tridiagonal matrix
59  with at most one non-zero sub-diagonal elements per column and U is
60  an upper triangular matrix with at most two non-zero super-diagonal
61  elements per column.
62 
63  The factorization is obtained by Gaussian elimination with partial
64  pivoting and implicit row scaling.
65 
66  The parameter LAMBDA is included in the routine so that DLAGTF may
67  be used, in conjunction with DLAGTS, to obtain eigenvectors of T by
68  inverse iteration.
69 
70  Arguments
71  =========
72 
73  N (input) INTEGER
74  The order of the matrix T.
75 
76  A (input/output) DOUBLE PRECISION array, dimension (N)
77  On entry, A must contain the diagonal elements of T.
78 
79  On exit, A is overwritten by the n diagonal elements of the
80  upper triangular matrix U of the factorization of T.
81 
82  LAMBDA (input) DOUBLE PRECISION
83  On entry, the scalar lambda.
84 
85  B (input/output) DOUBLE PRECISION array, dimension (N-1)
86  On entry, B must contain the (n-1) super-diagonal elements of
87  T.
88 
89  On exit, B is overwritten by the (n-1) super-diagonal
90  elements of the matrix U of the factorization of T.
91 
92  C (input/output) DOUBLE PRECISION array, dimension (N-1)
93  On entry, C must contain the (n-1) sub-diagonal elements of
94  T.
95 
96  On exit, C is overwritten by the (n-1) sub-diagonal elements
97  of the matrix L of the factorization of T.
98 
99  TOL (input) DOUBLE PRECISION
100  On entry, a relative tolerance used to indicate whether or
101  not the matrix (T - lambda*I) is nearly singular. TOL should
102  normally be chose as approximately the largest relative error
103  in the elements of T. For example, if the elements of T are
104  correct to about 4 significant figures, then TOL should be
105  set to about 5*10**(-4). If TOL is supplied as less than eps,
106  where eps is the relative machine precision, then the value
107  eps is used in place of TOL.
108 
109  D (output) DOUBLE PRECISION array, dimension (N-2)
110  On exit, D is overwritten by the (n-2) second super-diagonal
111  elements of the matrix U of the factorization of T.
112 
113  IN (output) INTEGER array, dimension (N)
114  On exit, IN contains details of the permutation matrix P. If
115  an interchange occurred at the kth step of the elimination,
116  then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)
117  returns the smallest positive integer j such that
118 
119  abs( u(j,j) ).le. norm( (T - lambda*I)(j) )*TOL,
120 
121  where norm( A(j) ) denotes the sum of the absolute values of
122  the jth row of the matrix A. If no such j exists then IN(n)
123  is returned as zero. If IN(n) is returned as positive, then a
124  diagonal element of U is small, indicating that
125  (T - lambda*I) is singular or nearly singular,
126 
127  INFO (output) INTEGER
128  = 0 : successful exit
129  .lt. 0: if INFO = -k, the kth argument had an illegal value
130 
131  =====================================================================
132 
133 
134  Parameter adjustments */
135  /* System generated locals */
136  integer i__1;
137  Treal d__1, d__2;
138  /* Local variables */
139  Treal temp, mult;
140  integer k;
141  Treal scale1, scale2;
142  Treal tl;
143  Treal eps, piv1, piv2;
144 
145  --in;
146  --d__;
147  --c__;
148  --b;
149  --a;
150 
151  /* Function Body */
152  *info = 0;
153  if (*n < 0) {
154  *info = -1;
155  i__1 = -(*info);
156  template_blas_erbla("LAGTF ", &i__1);
157  return 0;
158  }
159 
160  if (*n == 0) {
161  return 0;
162  }
163 
164  a[1] -= *lambda;
165  in[*n] = 0;
166  if (*n == 1) {
167  if (a[1] == 0.) {
168  in[1] = 1;
169  }
170  return 0;
171  }
172 
173  eps = template_lapack_lamch("Epsilon", (Treal)0);
174 
175  tl = maxMACRO(*tol,eps);
176  scale1 = absMACRO(a[1]) + absMACRO(b[1]);
177  i__1 = *n - 1;
178  for (k = 1; k <= i__1; ++k) {
179  a[k + 1] -= *lambda;
180  scale2 = (d__1 = c__[k], absMACRO(d__1)) + (d__2 = a[k + 1], absMACRO(d__2));
181  if (k < *n - 1) {
182  scale2 += (d__1 = b[k + 1], absMACRO(d__1));
183  }
184  if (a[k] == 0.) {
185  piv1 = 0.;
186  } else {
187  piv1 = (d__1 = a[k], absMACRO(d__1)) / scale1;
188  }
189  if (c__[k] == 0.) {
190  in[k] = 0;
191  piv2 = 0.;
192  scale1 = scale2;
193  if (k < *n - 1) {
194  d__[k] = 0.;
195  }
196  } else {
197  piv2 = (d__1 = c__[k], absMACRO(d__1)) / scale2;
198  if (piv2 <= piv1) {
199  in[k] = 0;
200  scale1 = scale2;
201  c__[k] /= a[k];
202  a[k + 1] -= c__[k] * b[k];
203  if (k < *n - 1) {
204  d__[k] = 0.;
205  }
206  } else {
207  in[k] = 1;
208  mult = a[k] / c__[k];
209  a[k] = c__[k];
210  temp = a[k + 1];
211  a[k + 1] = b[k] - mult * temp;
212  if (k < *n - 1) {
213  d__[k] = b[k + 1];
214  b[k + 1] = -mult * d__[k];
215  }
216  b[k] = temp;
217  c__[k] = mult;
218  }
219  }
220  if (maxMACRO(piv1,piv2) <= tl && in[*n] == 0) {
221  in[*n] = k;
222  }
223 /* L10: */
224  }
225  if ((d__1 = a[*n], absMACRO(d__1)) <= scale1 * tl && in[*n] == 0) {
226  in[*n] = *n;
227  }
228 
229  return 0;
230 
231 /* End of DLAGTF */
232 
233 } /* dlagtf_ */
234 
235 #endif
#define absMACRO(x)
Definition: template_blas_common.h:45
int integer
Definition: template_blas_common.h:38
#define maxMACRO(a, b)
Definition: template_blas_common.h:43
int template_lapack_lagtf(const integer *n, Treal *a, const Treal *lambda, Treal *b, Treal *c__, const Treal *tol, Treal *d__, integer *in, integer *info)
Definition: template_lapack_lagtf.h:40
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:199