ergo
template_lapack_lacon.h
Go to the documentation of this file.
1 /* Ergo, version 3.3, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2013 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28  /* This file belongs to the template_lapack part of the Ergo source
29  * code. The source files in the template_lapack directory are modified
30  * versions of files originally distributed as CLAPACK, see the
31  * Copyright/license notice in the file template_lapack/COPYING.
32  */
33 
34 
35 #ifndef TEMPLATE_LAPACK_LACON_HEADER
36 #define TEMPLATE_LAPACK_LACON_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_lacon(const integer *n, Treal *v, Treal *x,
41  integer *isgn, Treal *est, integer *kase)
42 {
43 /* -- LAPACK auxiliary routine (version 3.0) --
44  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
45  Courant Institute, Argonne National Lab, and Rice University
46  February 29, 1992
47 
48 
49  Purpose
50  =======
51 
52  DLACON estimates the 1-norm of a square, real matrix A.
53  Reverse communication is used for evaluating matrix-vector products.
54 
55  Arguments
56  =========
57 
58  N (input) INTEGER
59  The order of the matrix. N >= 1.
60 
61  V (workspace) DOUBLE PRECISION array, dimension (N)
62  On the final return, V = A*W, where EST = norm(V)/norm(W)
63  (W is not returned).
64 
65  X (input/output) DOUBLE PRECISION array, dimension (N)
66  On an intermediate return, X should be overwritten by
67  A * X, if KASE=1,
68  A' * X, if KASE=2,
69  and DLACON must be re-called with all the other parameters
70  unchanged.
71 
72  ISGN (workspace) INTEGER array, dimension (N)
73 
74  EST (output) DOUBLE PRECISION
75  An estimate (a lower bound) for norm(A).
76 
77  KASE (input/output) INTEGER
78  On the initial call to DLACON, KASE should be 0.
79  On an intermediate return, KASE will be 1 or 2, indicating
80  whether X should be overwritten by A * X or A' * X.
81  On the final return from DLACON, KASE will again be 0.
82 
83  Further Details
84  ======= =======
85 
86  Contributed by Nick Higham, University of Manchester.
87  Originally named SONEST, dated March 16, 1988.
88 
89  Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of
90  a real or complex matrix, with applications to condition estimation",
91  ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
92 
93  =====================================================================
94 
95 
96  Parameter adjustments */
97  /* Table of constant values */
98  integer c__1 = 1;
99  Treal c_b11 = 1.;
100 
101  /* System generated locals */
102  integer i__1;
103  Treal d__1;
104  /* Builtin functions */
105  double d_sign(Treal *, Treal *);
106  integer i_dnnt(Treal *);
107  /* Local variables */
108  integer iter;
109  Treal temp;
110  integer jump, i__, j;
111  integer jlast;
112  Treal altsgn, estold;
113 
114 
115  --isgn;
116  --x;
117  --v;
118 
119  /* Function Body */
120  if (*kase == 0) {
121  i__1 = *n;
122  for (i__ = 1; i__ <= i__1; ++i__) {
123  x[i__] = 1. / (Treal) (*n);
124 /* L10: */
125  }
126  *kase = 1;
127  jump = 1;
128  return 0;
129  }
130 
131  switch (jump) {
132  case 1: goto L20;
133  case 2: goto L40;
134  case 3: goto L70;
135  case 4: goto L110;
136  case 5: goto L140;
137  }
138 
139 /* ................ ENTRY (JUMP = 1)
140  FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X. */
141 
142 L20:
143  if (*n == 1) {
144  v[1] = x[1];
145  *est = absMACRO(v[1]);
146 /* ... QUIT */
147  goto L150;
148  }
149  *est = template_blas_asum(n, &x[1], &c__1);
150 
151  i__1 = *n;
152  for (i__ = 1; i__ <= i__1; ++i__) {
153  x[i__] = d_sign(&c_b11, &x[i__]);
154  isgn[i__] = i_dnnt(&x[i__]);
155 /* L30: */
156  }
157  *kase = 2;
158  jump = 2;
159  return 0;
160 
161 /* ................ ENTRY (JUMP = 2)
162  FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. */
163 
164 L40:
165  j = template_blas_idamax(n, &x[1], &c__1);
166  iter = 2;
167 
168 /* MAIN LOOP - ITERATIONS 2,3,...,ITMAX. */
169 
170 L50:
171  i__1 = *n;
172  for (i__ = 1; i__ <= i__1; ++i__) {
173  x[i__] = 0.;
174 /* L60: */
175  }
176  x[j] = 1.;
177  *kase = 1;
178  jump = 3;
179  return 0;
180 
181 /* ................ ENTRY (JUMP = 3)
182  X HAS BEEN OVERWRITTEN BY A*X. */
183 
184 L70:
185  template_blas_copy(n, &x[1], &c__1, &v[1], &c__1);
186  estold = *est;
187  *est = template_blas_asum(n, &v[1], &c__1);
188  i__1 = *n;
189  for (i__ = 1; i__ <= i__1; ++i__) {
190  d__1 = d_sign(&c_b11, &x[i__]);
191  if (i_dnnt(&d__1) != isgn[i__]) {
192  goto L90;
193  }
194 /* L80: */
195  }
196 /* REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED. */
197  goto L120;
198 
199 L90:
200 /* TEST FOR CYCLING. */
201  if (*est <= estold) {
202  goto L120;
203  }
204 
205  i__1 = *n;
206  for (i__ = 1; i__ <= i__1; ++i__) {
207  x[i__] = d_sign(&c_b11, &x[i__]);
208  isgn[i__] = i_dnnt(&x[i__]);
209 /* L100: */
210  }
211  *kase = 2;
212  jump = 4;
213  return 0;
214 
215 /* ................ ENTRY (JUMP = 4)
216  X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. */
217 
218 L110:
219  jlast = j;
220  j = template_blas_idamax(n, &x[1], &c__1);
221  if (x[jlast] != (d__1 = x[j], absMACRO(d__1)) && iter < 5) {
222  ++iter;
223  goto L50;
224  }
225 
226 /* ITERATION COMPLETE. FINAL STAGE. */
227 
228 L120:
229  altsgn = 1.;
230  i__1 = *n;
231  for (i__ = 1; i__ <= i__1; ++i__) {
232  x[i__] = altsgn * ((Treal) (i__ - 1) / (Treal) (*n - 1) +
233  1.);
234  altsgn = -altsgn;
235 /* L130: */
236  }
237  *kase = 1;
238  jump = 5;
239  return 0;
240 
241 /* ................ ENTRY (JUMP = 5)
242  X HAS BEEN OVERWRITTEN BY A*X. */
243 
244 L140:
245  temp = template_blas_asum(n, &x[1], &c__1) / (Treal) (*n * 3) * 2.;
246  if (temp > *est) {
247  template_blas_copy(n, &x[1], &c__1, &v[1], &c__1);
248  *est = temp;
249  }
250 
251 L150:
252  *kase = 0;
253  return 0;
254 
255 /* End of DLACON */
256 
257 } /* dlacon_ */
258 
259 #endif
#define absMACRO(x)
Definition: template_blas_common.h:45
integer template_blas_idamax(const integer *n, const Treal *dx, const integer *incx)
Definition: template_blas_idamax.h:40
int integer
Definition: template_blas_common.h:38
int template_lapack_lacon(const integer *n, Treal *v, Treal *x, integer *isgn, Treal *est, integer *kase)
Definition: template_lapack_lacon.h:40
Treal template_blas_asum(const integer *n, const Treal *dx, const integer *incx)
Definition: template_blas_asum.h:40
int template_blas_copy(const integer *n, const Treal *dx, const integer *incx, Treal *dy, const integer *incy)
Definition: template_blas_copy.h:40