ergo
template_lapack_rscl.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_RSCL_HEADER
36 #define TEMPLATE_LAPACK_RSCL_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_rscl(const integer *n, const Treal *sa, Treal *sx,
41  const integer *incx)
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  September 30, 1994
47 
48 
49  Purpose
50  =======
51 
52  DRSCL multiplies an n-element real vector x by the real scalar 1/a.
53  This is done without overflow or underflow as long as
54  the final result x/a does not overflow or underflow.
55 
56  Arguments
57  =========
58 
59  N (input) INTEGER
60  The number of components of the vector x.
61 
62  SA (input) DOUBLE PRECISION
63  The scalar a which is used to divide each component of x.
64  SA must be >= 0, or the subroutine will divide by zero.
65 
66  SX (input/output) DOUBLE PRECISION array, dimension
67  (1+(N-1)*abs(INCX))
68  The n-element vector x.
69 
70  INCX (input) INTEGER
71  The increment between successive values of the vector SX.
72  > 0: SX(1) = X(1) and SX(1+(i-1)*INCX) = x(i), 1< i<= n
73 
74  =====================================================================
75 
76 
77  Quick return if possible
78 
79  Parameter adjustments */
80  Treal cden;
81  logical done;
82  Treal cnum, cden1, cnum1;
83  Treal bignum, smlnum, mul;
84 
85  --sx;
86 
87  /* Function Body */
88  if (*n <= 0) {
89  return 0;
90  }
91 
92 /* Get machine parameters */
93 
94  smlnum = template_lapack_lamch("S", (Treal)0);
95  bignum = 1. / smlnum;
96  template_lapack_labad(&smlnum, &bignum);
97 
98 /* Initialize the denominator to SA and the numerator to 1. */
99 
100  cden = *sa;
101  cnum = 1.;
102 
103 L10:
104  cden1 = cden * smlnum;
105  cnum1 = cnum / bignum;
106  if (absMACRO(cden1) > absMACRO(cnum) && cnum != 0.) {
107 
108 /* Pre-multiply X by SMLNUM if CDEN is large compared to CNUM. */
109 
110  mul = smlnum;
111  done = FALSE_;
112  cden = cden1;
113  } else if (absMACRO(cnum1) > absMACRO(cden)) {
114 
115 /* Pre-multiply X by BIGNUM if CDEN is small compared to CNUM. */
116 
117  mul = bignum;
118  done = FALSE_;
119  cnum = cnum1;
120  } else {
121 
122 /* Multiply X by CNUM / CDEN and return. */
123 
124  mul = cnum / cden;
125  done = TRUE_;
126  }
127 
128 /* Scale the vector X by MUL */
129 
130  dscal_(n, &mul, &sx[1], incx);
131 
132  if (! done) {
133  goto L10;
134  }
135 
136  return 0;
137 
138 /* End of DRSCL */
139 
140 } /* drscl_ */
141 
142 #endif