ergo
template_lapack_lacon.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_LACON_HEADER
38#define TEMPLATE_LAPACK_LACON_HEADER
39
40
41template<class Treal>
42int template_lapack_lacon(const integer *n, Treal *v, Treal *x,
43 integer *isgn, Treal *est, integer *kase)
44{
45/* -- LAPACK auxiliary routine (version 3.0) --
46 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
47 Courant Institute, Argonne National Lab, and Rice University
48 February 29, 1992
49
50
51 Purpose
52 =======
53
54 DLACON estimates the 1-norm of a square, real matrix A.
55 Reverse communication is used for evaluating matrix-vector products.
56
57 Arguments
58 =========
59
60 N (input) INTEGER
61 The order of the matrix. N >= 1.
62
63 V (workspace) DOUBLE PRECISION array, dimension (N)
64 On the final return, V = A*W, where EST = norm(V)/norm(W)
65 (W is not returned).
66
67 X (input/output) DOUBLE PRECISION array, dimension (N)
68 On an intermediate return, X should be overwritten by
69 A * X, if KASE=1,
70 A' * X, if KASE=2,
71 and DLACON must be re-called with all the other parameters
72 unchanged.
73
74 ISGN (workspace) INTEGER array, dimension (N)
75
76 EST (output) DOUBLE PRECISION
77 An estimate (a lower bound) for norm(A).
78
79 KASE (input/output) INTEGER
80 On the initial call to DLACON, KASE should be 0.
81 On an intermediate return, KASE will be 1 or 2, indicating
82 whether X should be overwritten by A * X or A' * X.
83 On the final return from DLACON, KASE will again be 0.
84
85 Further Details
86 ======= =======
87
88 Contributed by Nick Higham, University of Manchester.
89 Originally named SONEST, dated March 16, 1988.
90
91 Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of
92 a real or complex matrix, with applications to condition estimation",
93 ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
94
95 =====================================================================
96
97
98 Parameter adjustments */
99 /* Table of constant values */
100 integer c__1 = 1;
101 Treal c_b11 = 1.;
102
103 /* System generated locals */
104 integer i__1;
105 Treal d__1;
106 /* Builtin functions */
107 double d_sign(Treal *, Treal *);
108 integer i_dnnt(Treal *);
109 /* Local variables */
110 integer iter;
111 Treal temp;
112 integer jump, i__, j;
113 integer jlast;
114 Treal altsgn, estold;
115
116
117 --isgn;
118 --x;
119 --v;
120
121 /* Function Body */
122 if (*kase == 0) {
123 i__1 = *n;
124 for (i__ = 1; i__ <= i__1; ++i__) {
125 x[i__] = 1. / (Treal) (*n);
126/* L10: */
127 }
128 *kase = 1;
129 jump = 1;
130 return 0;
131 }
132
133 switch (jump) {
134 case 1: goto L20;
135 case 2: goto L40;
136 case 3: goto L70;
137 case 4: goto L110;
138 case 5: goto L140;
139 }
140
141/* ................ ENTRY (JUMP = 1)
142 FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X. */
143
144L20:
145 if (*n == 1) {
146 v[1] = x[1];
147 *est = absMACRO(v[1]);
148/* ... QUIT */
149 goto L150;
150 }
151 *est = template_blas_asum(n, &x[1], &c__1);
152
153 i__1 = *n;
154 for (i__ = 1; i__ <= i__1; ++i__) {
155 x[i__] = d_sign(&c_b11, &x[i__]);
156 isgn[i__] = i_dnnt(&x[i__]);
157/* L30: */
158 }
159 *kase = 2;
160 jump = 2;
161 return 0;
162
163/* ................ ENTRY (JUMP = 2)
164 FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. */
165
166L40:
167 j = template_blas_idamax(n, &x[1], &c__1);
168 iter = 2;
169
170/* MAIN LOOP - ITERATIONS 2,3,...,ITMAX. */
171
172L50:
173 i__1 = *n;
174 for (i__ = 1; i__ <= i__1; ++i__) {
175 x[i__] = 0.;
176/* L60: */
177 }
178 x[j] = 1.;
179 *kase = 1;
180 jump = 3;
181 return 0;
182
183/* ................ ENTRY (JUMP = 3)
184 X HAS BEEN OVERWRITTEN BY A*X. */
185
186L70:
187 template_blas_copy(n, &x[1], &c__1, &v[1], &c__1);
188 estold = *est;
189 *est = template_blas_asum(n, &v[1], &c__1);
190 i__1 = *n;
191 for (i__ = 1; i__ <= i__1; ++i__) {
192 d__1 = d_sign(&c_b11, &x[i__]);
193 if (i_dnnt(&d__1) != isgn[i__]) {
194 goto L90;
195 }
196/* L80: */
197 }
198/* REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED. */
199 goto L120;
200
201L90:
202/* TEST FOR CYCLING. */
203 if (*est <= estold) {
204 goto L120;
205 }
206
207 i__1 = *n;
208 for (i__ = 1; i__ <= i__1; ++i__) {
209 x[i__] = d_sign(&c_b11, &x[i__]);
210 isgn[i__] = i_dnnt(&x[i__]);
211/* L100: */
212 }
213 *kase = 2;
214 jump = 4;
215 return 0;
216
217/* ................ ENTRY (JUMP = 4)
218 X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. */
219
220L110:
221 jlast = j;
222 j = template_blas_idamax(n, &x[1], &c__1);
223 if (x[jlast] != (d__1 = x[j], absMACRO(d__1)) && iter < 5) {
224 ++iter;
225 goto L50;
226 }
227
228/* ITERATION COMPLETE. FINAL STAGE. */
229
230L120:
231 altsgn = 1.;
232 i__1 = *n;
233 for (i__ = 1; i__ <= i__1; ++i__) {
234 x[i__] = altsgn * ((Treal) (i__ - 1) / (Treal) (*n - 1) +
235 1.);
236 altsgn = -altsgn;
237/* L130: */
238 }
239 *kase = 1;
240 jump = 5;
241 return 0;
242
243/* ................ ENTRY (JUMP = 5)
244 X HAS BEEN OVERWRITTEN BY A*X. */
245
246L140:
247 temp = template_blas_asum(n, &x[1], &c__1) / (Treal) (*n * 3) * 2.;
248 if (temp > *est) {
249 template_blas_copy(n, &x[1], &c__1, &v[1], &c__1);
250 *est = temp;
251 }
252
253L150:
254 *kase = 0;
255 return 0;
256
257/* End of DLACON */
258
259} /* dlacon_ */
260
261#endif
Treal template_blas_asum(const integer *n, const Treal *dx, const integer *incx)
Definition: template_blas_asum.h:42
int integer
Definition: template_blas_common.h:40
#define absMACRO(x)
Definition: template_blas_common.h:47
int template_blas_copy(const integer *n, const Treal *dx, const integer *incx, Treal *dy, const integer *incy)
Definition: template_blas_copy.h:42
integer template_blas_idamax(const integer *n, const Treal *dx, const integer *incx)
Definition: template_blas_idamax.h:42
int template_lapack_lacon(const integer *n, Treal *v, Treal *x, integer *isgn, Treal *est, integer *kase)
Definition: template_lapack_lacon.h:42