ergo
template_lapack_lacon.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_LACON_HEADER
38 #define TEMPLATE_LAPACK_LACON_HEADER
39 
40 
41 template<class Treal>
42 int 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 
144 L20:
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 
166 L40:
167  j = template_blas_idamax(n, &x[1], &c__1);
168  iter = 2;
169 
170 /* MAIN LOOP - ITERATIONS 2,3,...,ITMAX. */
171 
172 L50:
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 
186 L70:
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 
201 L90:
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 
220 L110:
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 
230 L120:
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 
246 L140:
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 
253 L150:
254  *kase = 0;
255  return 0;
256 
257 /* End of DLACON */
258 
259 } /* dlacon_ */
260 
261 #endif
template_blas_idamax
integer template_blas_idamax(const integer *n, const Treal *dx, const integer *incx)
Definition: template_blas_idamax.h:42
template_lapack_lacon
int template_lapack_lacon(const integer *n, Treal *v, Treal *x, integer *isgn, Treal *est, integer *kase)
Definition: template_lapack_lacon.h:42
absMACRO
#define absMACRO(x)
Definition: template_blas_common.h:47
template_blas_asum
Treal template_blas_asum(const integer *n, const Treal *dx, const integer *incx)
Definition: template_blas_asum.h:42
template_blas_copy
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
int integer
Definition: template_blas_common.h:40