ergo
template_lapack_larra.h
Go to the documentation of this file.
1 /* Ergo, version 3.2, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2012 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_LARRA_HEADER
36 #define TEMPLATE_LAPACK_LARRA_HEADER
37 
38 template<class Treal>
39 int template_lapack_larra(const integer *n, Treal *d__, Treal *e,
40  Treal *e2, Treal *spltol, Treal *tnrm, integer *nsplit,
41  integer *isplit, integer *info)
42 {
43  /* System generated locals */
44  integer i__1;
45  Treal d__1, d__2;
46 
47 
48  /* Local variables */
49  integer i__;
50  Treal tmp1, eabs;
51 
52 
53 /* -- LAPACK auxiliary routine (version 3.2) -- */
54 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
55 /* November 2006 */
56 
57 /* .. Scalar Arguments .. */
58 /* .. */
59 /* .. Array Arguments .. */
60 /* .. */
61 
62 /* Purpose */
63 /* ======= */
64 
65 /* Compute the splitting points with threshold SPLTOL. */
66 /* DLARRA sets any "small" off-diagonal elements to zero. */
67 
68 /* Arguments */
69 /* ========= */
70 
71 /* N (input) INTEGER */
72 /* The order of the matrix. N > 0. */
73 
74 /* D (input) DOUBLE PRECISION array, dimension (N) */
75 /* On entry, the N diagonal elements of the tridiagonal */
76 /* matrix T. */
77 
78 /* E (input/output) DOUBLE PRECISION array, dimension (N) */
79 /* On entry, the first (N-1) entries contain the subdiagonal */
80 /* elements of the tridiagonal matrix T; E(N) need not be set. */
81 /* On exit, the entries E( ISPLIT( I ) ), 1 <= I <= NSPLIT, */
82 /* are set to zero, the other entries of E are untouched. */
83 
84 /* E2 (input/output) DOUBLE PRECISION array, dimension (N) */
85 /* On entry, the first (N-1) entries contain the SQUARES of the */
86 /* subdiagonal elements of the tridiagonal matrix T; */
87 /* E2(N) need not be set. */
88 /* On exit, the entries E2( ISPLIT( I ) ), */
89 /* 1 <= I <= NSPLIT, have been set to zero */
90 
91 /* SPLTOL (input) DOUBLE PRECISION */
92 /* The threshold for splitting. Two criteria can be used: */
93 /* SPLTOL<0 : criterion based on absolute off-diagonal value */
94 /* SPLTOL>0 : criterion that preserves relative accuracy */
95 
96 /* TNRM (input) DOUBLE PRECISION */
97 /* The norm of the matrix. */
98 
99 /* NSPLIT (output) INTEGER */
100 /* The number of blocks T splits into. 1 <= NSPLIT <= N. */
101 
102 /* ISPLIT (output) INTEGER array, dimension (N) */
103 /* The splitting points, at which T breaks up into blocks. */
104 /* The first block consists of rows/columns 1 to ISPLIT(1), */
105 /* the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), */
106 /* etc., and the NSPLIT-th consists of rows/columns */
107 /* ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. */
108 
109 
110 /* INFO (output) INTEGER */
111 /* = 0: successful exit */
112 
113 /* Further Details */
114 /* =============== */
115 
116 /* Based on contributions by */
117 /* Beresford Parlett, University of California, Berkeley, USA */
118 /* Jim Demmel, University of California, Berkeley, USA */
119 /* Inderjit Dhillon, University of Texas, Austin, USA */
120 /* Osni Marques, LBNL/NERSC, USA */
121 /* Christof Voemel, University of California, Berkeley, USA */
122 
123 /* ===================================================================== */
124 
125 /* .. Parameters .. */
126 /* .. */
127 /* .. Local Scalars .. */
128 /* .. */
129 /* .. Intrinsic Functions .. */
130 /* .. */
131 /* .. Executable Statements .. */
132 
133  /* Parameter adjustments */
134  --isplit;
135  --e2;
136  --e;
137  --d__;
138 
139  /* Function Body */
140  *info = 0;
141 /* Compute splitting points */
142  *nsplit = 1;
143  if (*spltol < 0.) {
144 /* Criterion based on absolute off-diagonal value */
145  tmp1 = absMACRO(*spltol) * *tnrm;
146  i__1 = *n - 1;
147  for (i__ = 1; i__ <= i__1; ++i__) {
148  eabs = (d__1 = e[i__], absMACRO(d__1));
149  if (eabs <= tmp1) {
150  e[i__] = 0.;
151  e2[i__] = 0.;
152  isplit[*nsplit] = i__;
153  ++(*nsplit);
154  }
155 /* L9: */
156  }
157  } else {
158 /* Criterion that guarantees relative accuracy */
159  i__1 = *n - 1;
160  for (i__ = 1; i__ <= i__1; ++i__) {
161  eabs = (d__1 = e[i__], absMACRO(d__1));
162  if (eabs <= *spltol * template_blas_sqrt((d__1 = d__[i__], absMACRO(d__1))) * template_blas_sqrt((
163  d__2 = d__[i__ + 1], absMACRO(d__2)))) {
164  e[i__] = 0.;
165  e2[i__] = 0.;
166  isplit[*nsplit] = i__;
167  ++(*nsplit);
168  }
169 /* L10: */
170  }
171  }
172  isplit[*nsplit] = *n;
173  return 0;
174 
175 /* End of DLARRA */
176 
177 } /* dlarra_ */
178 
179 #endif