ergo
template_lapack_lasq6.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 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_LASQ6_HEADER
38#define TEMPLATE_LAPACK_LASQ6_HEADER
39
40template<class Treal>
41int template_lapack_lasq6(integer *i0, integer *n0, Treal *z__,
42 integer *pp, Treal *dmin__, Treal *dmin1, Treal *dmin2,
43 Treal *dn, Treal *dnm1, Treal *dnm2)
44{
45 /* System generated locals */
46 integer i__1;
47 Treal d__1, d__2;
48
49 /* Local variables */
50 Treal d__;
51 integer j4, j4p2;
52 Treal emin, temp;
53 Treal safmin;
54
55
56/* -- LAPACK routine (version 3.2) -- */
57
58/* -- Contributed by Osni Marques of the Lawrence Berkeley National -- */
59/* -- Laboratory and Beresford Parlett of the Univ. of California at -- */
60/* -- Berkeley -- */
61/* -- November 2008 -- */
62
63/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
64/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
65
66/* .. Scalar Arguments .. */
67/* .. */
68/* .. Array Arguments .. */
69/* .. */
70
71/* Purpose */
72/* ======= */
73
74/* DLASQ6 computes one dqd (shift equal to zero) transform in */
75/* ping-pong form, with protection against underflow and overflow. */
76
77/* Arguments */
78/* ========= */
79
80/* I0 (input) INTEGER */
81/* First index. */
82
83/* N0 (input) INTEGER */
84/* Last index. */
85
86/* Z (input) DOUBLE PRECISION array, dimension ( 4*N ) */
87/* Z holds the qd array. EMIN is stored in Z(4*N0) to avoid */
88/* an extra argument. */
89
90/* PP (input) INTEGER */
91/* PP=0 for ping, PP=1 for pong. */
92
93/* DMIN (output) DOUBLE PRECISION */
94/* Minimum value of d. */
95
96/* DMIN1 (output) DOUBLE PRECISION */
97/* Minimum value of d, excluding D( N0 ). */
98
99/* DMIN2 (output) DOUBLE PRECISION */
100/* Minimum value of d, excluding D( N0 ) and D( N0-1 ). */
101
102/* DN (output) DOUBLE PRECISION */
103/* d(N0), the last value of d. */
104
105/* DNM1 (output) DOUBLE PRECISION */
106/* d(N0-1). */
107
108/* DNM2 (output) DOUBLE PRECISION */
109/* d(N0-2). */
110
111/* ===================================================================== */
112
113/* .. Parameter .. */
114/* .. */
115/* .. Local Scalars .. */
116/* .. */
117/* .. External Function .. */
118/* .. */
119/* .. Intrinsic Functions .. */
120/* .. */
121/* .. Executable Statements .. */
122
123 /* Parameter adjustments */
124 --z__;
125
126 /* Function Body */
127 if (*n0 - *i0 - 1 <= 0) {
128 return 0;
129 }
130
131 safmin = template_lapack_lamch("Safe minimum", (Treal)0);
132 j4 = (*i0 << 2) + *pp - 3;
133 emin = z__[j4 + 4];
134 d__ = z__[j4];
135 *dmin__ = d__;
136
137 if (*pp == 0) {
138 i__1 = (*n0 - 3) << 2;
139 for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
140 z__[j4 - 2] = d__ + z__[j4 - 1];
141 if (z__[j4 - 2] == 0.) {
142 z__[j4] = 0.;
143 d__ = z__[j4 + 1];
144 *dmin__ = d__;
145 emin = 0.;
146 } else if (safmin * z__[j4 + 1] < z__[j4 - 2] && safmin * z__[j4
147 - 2] < z__[j4 + 1]) {
148 temp = z__[j4 + 1] / z__[j4 - 2];
149 z__[j4] = z__[j4 - 1] * temp;
150 d__ *= temp;
151 } else {
152 z__[j4] = z__[j4 + 1] * (z__[j4 - 1] / z__[j4 - 2]);
153 d__ = z__[j4 + 1] * (d__ / z__[j4 - 2]);
154 }
155 *dmin__ = minMACRO(*dmin__,d__);
156/* Computing MIN */
157 d__1 = emin, d__2 = z__[j4];
158 emin = minMACRO(d__1,d__2);
159/* L10: */
160 }
161 } else {
162 i__1 = ( *n0 - 3 ) << 2;
163 for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
164 z__[j4 - 3] = d__ + z__[j4];
165 if (z__[j4 - 3] == 0.) {
166 z__[j4 - 1] = 0.;
167 d__ = z__[j4 + 2];
168 *dmin__ = d__;
169 emin = 0.;
170 } else if (safmin * z__[j4 + 2] < z__[j4 - 3] && safmin * z__[j4
171 - 3] < z__[j4 + 2]) {
172 temp = z__[j4 + 2] / z__[j4 - 3];
173 z__[j4 - 1] = z__[j4] * temp;
174 d__ *= temp;
175 } else {
176 z__[j4 - 1] = z__[j4 + 2] * (z__[j4] / z__[j4 - 3]);
177 d__ = z__[j4 + 2] * (d__ / z__[j4 - 3]);
178 }
179 *dmin__ = minMACRO(*dmin__,d__);
180/* Computing MIN */
181 d__1 = emin, d__2 = z__[j4 - 1];
182 emin = minMACRO(d__1,d__2);
183/* L20: */
184 }
185 }
186
187/* Unroll last two steps. */
188
189 *dnm2 = d__;
190 *dmin2 = *dmin__;
191 j4 = ( ( *n0 - 2 ) << 2) - *pp;
192 j4p2 = j4 + (*pp << 1) - 1;
193 z__[j4 - 2] = *dnm2 + z__[j4p2];
194 if (z__[j4 - 2] == 0.) {
195 z__[j4] = 0.;
196 *dnm1 = z__[j4p2 + 2];
197 *dmin__ = *dnm1;
198 emin = 0.;
199 } else if (safmin * z__[j4p2 + 2] < z__[j4 - 2] && safmin * z__[j4 - 2] <
200 z__[j4p2 + 2]) {
201 temp = z__[j4p2 + 2] / z__[j4 - 2];
202 z__[j4] = z__[j4p2] * temp;
203 *dnm1 = *dnm2 * temp;
204 } else {
205 z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
206 *dnm1 = z__[j4p2 + 2] * (*dnm2 / z__[j4 - 2]);
207 }
208 *dmin__ = minMACRO(*dmin__,*dnm1);
209
210 *dmin1 = *dmin__;
211 j4 += 4;
212 j4p2 = j4 + (*pp << 1) - 1;
213 z__[j4 - 2] = *dnm1 + z__[j4p2];
214 if (z__[j4 - 2] == 0.) {
215 z__[j4] = 0.;
216 *dn = z__[j4p2 + 2];
217 *dmin__ = *dn;
218 emin = 0.;
219 } else if (safmin * z__[j4p2 + 2] < z__[j4 - 2] && safmin * z__[j4 - 2] <
220 z__[j4p2 + 2]) {
221 temp = z__[j4p2 + 2] / z__[j4 - 2];
222 z__[j4] = z__[j4p2] * temp;
223 *dn = *dnm1 * temp;
224 } else {
225 z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
226 *dn = z__[j4p2 + 2] * (*dnm1 / z__[j4 - 2]);
227 }
228 *dmin__ = minMACRO(*dmin__,*dn);
229
230 z__[j4 + 2] = *dn;
231 z__[(*n0 << 2) - *pp] = emin;
232 return 0;
233
234/* End of DLASQ6 */
235
236} /* dlasq6_ */
237
238#endif
int integer
Definition: template_blas_common.h:40
#define minMACRO(a, b)
Definition: template_blas_common.h:46
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
int template_lapack_lasq6(integer *i0, integer *n0, Treal *z__, integer *pp, Treal *dmin__, Treal *dmin1, Treal *dmin2, Treal *dn, Treal *dnm1, Treal *dnm2)
Definition: template_lapack_lasq6.h:41