ergo
template_lapack_larra.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_LARRA_HEADER
38 #define TEMPLATE_LAPACK_LARRA_HEADER
39 
40 template<class Treal>
41 int template_lapack_larra(const integer *n, Treal *d__, Treal *e,
42  Treal *e2, Treal *spltol, Treal *tnrm, integer *nsplit,
43  integer *isplit, integer *info)
44 {
45  /* System generated locals */
46  integer i__1;
47  Treal d__1, d__2;
48 
49 
50  /* Local variables */
51  integer i__;
52  Treal tmp1, eabs;
53 
54 
55 /* -- LAPACK auxiliary routine (version 3.2) -- */
56 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
57 /* November 2006 */
58 
59 /* .. Scalar Arguments .. */
60 /* .. */
61 /* .. Array Arguments .. */
62 /* .. */
63 
64 /* Purpose */
65 /* ======= */
66 
67 /* Compute the splitting points with threshold SPLTOL. */
68 /* DLARRA sets any "small" off-diagonal elements to zero. */
69 
70 /* Arguments */
71 /* ========= */
72 
73 /* N (input) INTEGER */
74 /* The order of the matrix. N > 0. */
75 
76 /* D (input) DOUBLE PRECISION array, dimension (N) */
77 /* On entry, the N diagonal elements of the tridiagonal */
78 /* matrix T. */
79 
80 /* E (input/output) DOUBLE PRECISION array, dimension (N) */
81 /* On entry, the first (N-1) entries contain the subdiagonal */
82 /* elements of the tridiagonal matrix T; E(N) need not be set. */
83 /* On exit, the entries E( ISPLIT( I ) ), 1 <= I <= NSPLIT, */
84 /* are set to zero, the other entries of E are untouched. */
85 
86 /* E2 (input/output) DOUBLE PRECISION array, dimension (N) */
87 /* On entry, the first (N-1) entries contain the SQUARES of the */
88 /* subdiagonal elements of the tridiagonal matrix T; */
89 /* E2(N) need not be set. */
90 /* On exit, the entries E2( ISPLIT( I ) ), */
91 /* 1 <= I <= NSPLIT, have been set to zero */
92 
93 /* SPLTOL (input) DOUBLE PRECISION */
94 /* The threshold for splitting. Two criteria can be used: */
95 /* SPLTOL<0 : criterion based on absolute off-diagonal value */
96 /* SPLTOL>0 : criterion that preserves relative accuracy */
97 
98 /* TNRM (input) DOUBLE PRECISION */
99 /* The norm of the matrix. */
100 
101 /* NSPLIT (output) INTEGER */
102 /* The number of blocks T splits into. 1 <= NSPLIT <= N. */
103 
104 /* ISPLIT (output) INTEGER array, dimension (N) */
105 /* The splitting points, at which T breaks up into blocks. */
106 /* The first block consists of rows/columns 1 to ISPLIT(1), */
107 /* the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), */
108 /* etc., and the NSPLIT-th consists of rows/columns */
109 /* ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. */
110 
111 
112 /* INFO (output) INTEGER */
113 /* = 0: successful exit */
114 
115 /* Further Details */
116 /* =============== */
117 
118 /* Based on contributions by */
119 /* Beresford Parlett, University of California, Berkeley, USA */
120 /* Jim Demmel, University of California, Berkeley, USA */
121 /* Inderjit Dhillon, University of Texas, Austin, USA */
122 /* Osni Marques, LBNL/NERSC, USA */
123 /* Christof Voemel, University of California, Berkeley, USA */
124 
125 /* ===================================================================== */
126 
127 /* .. Parameters .. */
128 /* .. */
129 /* .. Local Scalars .. */
130 /* .. */
131 /* .. Intrinsic Functions .. */
132 /* .. */
133 /* .. Executable Statements .. */
134 
135  /* Parameter adjustments */
136  --isplit;
137  --e2;
138  --e;
139  --d__;
140 
141  /* Function Body */
142  *info = 0;
143 /* Compute splitting points */
144  *nsplit = 1;
145  if (*spltol < 0.) {
146 /* Criterion based on absolute off-diagonal value */
147  tmp1 = absMACRO(*spltol) * *tnrm;
148  i__1 = *n - 1;
149  for (i__ = 1; i__ <= i__1; ++i__) {
150  eabs = (d__1 = e[i__], absMACRO(d__1));
151  if (eabs <= tmp1) {
152  e[i__] = 0.;
153  e2[i__] = 0.;
154  isplit[*nsplit] = i__;
155  ++(*nsplit);
156  }
157 /* L9: */
158  }
159  } else {
160 /* Criterion that guarantees relative accuracy */
161  i__1 = *n - 1;
162  for (i__ = 1; i__ <= i__1; ++i__) {
163  eabs = (d__1 = e[i__], absMACRO(d__1));
164  if (eabs <= *spltol * template_blas_sqrt((d__1 = d__[i__], absMACRO(d__1))) * template_blas_sqrt((
165  d__2 = d__[i__ + 1], absMACRO(d__2)))) {
166  e[i__] = 0.;
167  e2[i__] = 0.;
168  isplit[*nsplit] = i__;
169  ++(*nsplit);
170  }
171 /* L10: */
172  }
173  }
174  isplit[*nsplit] = *n;
175  return 0;
176 
177 /* End of DLARRA */
178 
179 } /* dlarra_ */
180 
181 #endif
template_blas_sqrt
Treal template_blas_sqrt(Treal x)
absMACRO
#define absMACRO(x)
Definition: template_blas_common.h:47
template_lapack_larra
int template_lapack_larra(const integer *n, Treal *d__, Treal *e, Treal *e2, Treal *spltol, Treal *tnrm, integer *nsplit, integer *isplit, integer *info)
Definition: template_lapack_larra.h:41
integer
int integer
Definition: template_blas_common.h:40