ergo
template_lapack_larfg.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_LARFG_HEADER
36 #define TEMPLATE_LAPACK_LARFG_HEADER
37 
38 #include "template_lapack_lapy2.h"
39 
40 template<class Treal>
41 int template_lapack_larfg(const integer *n, Treal *alpha, Treal *x,
42  const integer *incx, Treal *tau)
43 {
44 /* -- LAPACK auxiliary routine (version 3.0) --
45  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
46  Courant Institute, Argonne National Lab, and Rice University
47  September 30, 1994
48 
49 
50  Purpose
51  =======
52 
53  DLARFG generates a real elementary reflector H of order n, such
54  that
55 
56  H * ( alpha ) = ( beta ), H' * H = I.
57  ( x ) ( 0 )
58 
59  where alpha and beta are scalars, and x is an (n-1)-element real
60  vector. H is represented in the form
61 
62  H = I - tau * ( 1 ) * ( 1 v' ) ,
63  ( v )
64 
65  where tau is a real scalar and v is a real (n-1)-element
66  vector.
67 
68  If the elements of x are all zero, then tau = 0 and H is taken to be
69  the unit matrix.
70 
71  Otherwise 1 <= tau <= 2.
72 
73  Arguments
74  =========
75 
76  N (input) INTEGER
77  The order of the elementary reflector.
78 
79  ALPHA (input/output) DOUBLE PRECISION
80  On entry, the value alpha.
81  On exit, it is overwritten with the value beta.
82 
83  X (input/output) DOUBLE PRECISION array, dimension
84  (1+(N-2)*abs(INCX))
85  On entry, the vector x.
86  On exit, it is overwritten with the vector v.
87 
88  INCX (input) INTEGER
89  The increment between elements of X. INCX > 0.
90 
91  TAU (output) DOUBLE PRECISION
92  The value tau.
93 
94  =====================================================================
95 
96 
97  Parameter adjustments */
98  /* System generated locals */
99  integer i__1;
100  Treal d__1;
101  /* Local variables */
102  Treal beta;
103  integer j;
104  Treal xnorm;
105  Treal safmin, rsafmn;
106  integer knt;
107 
108  --x;
109 
110  /* Function Body */
111  if (*n <= 1) {
112  *tau = 0.;
113  return 0;
114  }
115 
116  i__1 = *n - 1;
117  xnorm = template_blas_nrm2(&i__1, &x[1], incx);
118 
119  if (xnorm == 0.) {
120 
121 /* H = I */
122 
123  *tau = 0.;
124  } else {
125 
126 /* general case */
127 
128  d__1 = template_lapack_lapy2(alpha, &xnorm);
129  beta = -template_lapack_d_sign(&d__1, alpha);
130  safmin = template_lapack_lamch("S", (Treal)0) / template_lapack_lamch("E", (Treal)0);
131  if (absMACRO(beta) < safmin) {
132 
133 /* XNORM, BETA may be inaccurate; scale X and recompute them */
134 
135  rsafmn = 1. / safmin;
136  knt = 0;
137 L10:
138  ++knt;
139  i__1 = *n - 1;
140  template_blas_scal(&i__1, &rsafmn, &x[1], incx);
141  beta *= rsafmn;
142  *alpha *= rsafmn;
143  if (absMACRO(beta) < safmin) {
144  goto L10;
145  }
146 
147 /* New BETA is at most 1, at least SAFMIN */
148 
149  i__1 = *n - 1;
150  xnorm = template_blas_nrm2(&i__1, &x[1], incx);
151  d__1 = template_lapack_lapy2(alpha, &xnorm);
152  beta = -template_lapack_d_sign(&d__1, alpha);
153  *tau = (beta - *alpha) / beta;
154  i__1 = *n - 1;
155  d__1 = 1. / (*alpha - beta);
156  template_blas_scal(&i__1, &d__1, &x[1], incx);
157 
158 /* If ALPHA is subnormal, it may lose relative accuracy */
159 
160  *alpha = beta;
161  i__1 = knt;
162  for (j = 1; j <= i__1; ++j) {
163  *alpha *= safmin;
164 /* L20: */
165  }
166  } else {
167  *tau = (beta - *alpha) / beta;
168  i__1 = *n - 1;
169  d__1 = 1. / (*alpha - beta);
170  template_blas_scal(&i__1, &d__1, &x[1], incx);
171  *alpha = beta;
172  }
173  }
174 
175  return 0;
176 
177 /* End of DLARFG */
178 
179 } /* dlarfg_ */
180 
181 #endif
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition: template_blas_scal.h:41
Treal template_blas_nrm2(const integer *n, const Treal *x, const integer *incx)
Definition: template_blas_nrm2.h:40
#define absMACRO(x)
Definition: template_blas_common.h:45
int integer
Definition: template_blas_common.h:38
Treal template_lapack_lapy2(Treal *x, Treal *y)
Definition: template_lapack_lapy2.h:40
Treal template_lapack_d_sign(const Treal *a, const Treal *b)
Definition: template_lapack_lamch.h:46
int template_lapack_larfg(const integer *n, Treal *alpha, Treal *x, const integer *incx, Treal *tau)
Definition: template_lapack_larfg.h:41
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:199