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