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
source
matrix
template_lapack
lapack
template_lapack_larra.h
Generated on Sun Aug 16 2020 00:00:00 for ergo by
1.8.18