ergo
template_lapack_org2r.h
Go to the documentation of this file.
1 /* Ergo, version 3.4, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2014 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_ORG2R_HEADER
36 #define TEMPLATE_LAPACK_ORG2R_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_org2r(const integer *m, const integer *n, const integer *k, Treal *
41  a, const integer *lda, const Treal *tau, Treal *work, integer *info)
42 {
43 /* -- LAPACK 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  DORG2R generates an m by n real matrix Q with orthonormal columns,
53  which is defined as the first n columns of a product of k elementary
54  reflectors of order m
55 
56  Q = H(1) H(2) . . . H(k)
57 
58  as returned by DGEQRF.
59 
60  Arguments
61  =========
62 
63  M (input) INTEGER
64  The number of rows of the matrix Q. M >= 0.
65 
66  N (input) INTEGER
67  The number of columns of the matrix Q. M >= N >= 0.
68 
69  K (input) INTEGER
70  The number of elementary reflectors whose product defines the
71  matrix Q. N >= K >= 0.
72 
73  A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
74  On entry, the i-th column must contain the vector which
75  defines the elementary reflector H(i), for i = 1,2,...,k, as
76  returned by DGEQRF in the first k columns of its array
77  argument A.
78  On exit, the m-by-n matrix Q.
79 
80  LDA (input) INTEGER
81  The first dimension of the array A. LDA >= max(1,M).
82 
83  TAU (input) DOUBLE PRECISION array, dimension (K)
84  TAU(i) must contain the scalar factor of the elementary
85  reflector H(i), as returned by DGEQRF.
86 
87  WORK (workspace) DOUBLE PRECISION array, dimension (N)
88 
89  INFO (output) INTEGER
90  = 0: successful exit
91  < 0: if INFO = -i, the i-th argument has an illegal value
92 
93  =====================================================================
94 
95 
96  Test the input arguments
97 
98  Parameter adjustments */
99  /* Table of constant values */
100  integer c__1 = 1;
101 
102  /* System generated locals */
103  integer a_dim1, a_offset, i__1, i__2;
104  Treal d__1;
105  /* Local variables */
106  integer i__, j, l;
107 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
108 
109 
110  a_dim1 = *lda;
111  a_offset = 1 + a_dim1 * 1;
112  a -= a_offset;
113  --tau;
114  --work;
115 
116  /* Function Body */
117  *info = 0;
118  if (*m < 0) {
119  *info = -1;
120  } else if (*n < 0 || *n > *m) {
121  *info = -2;
122  } else if (*k < 0 || *k > *n) {
123  *info = -3;
124  } else if (*lda < maxMACRO(1,*m)) {
125  *info = -5;
126  }
127  if (*info != 0) {
128  i__1 = -(*info);
129  template_blas_erbla("ORG2R ", &i__1);
130  return 0;
131  }
132 
133 /* Quick return if possible */
134 
135  if (*n <= 0) {
136  return 0;
137  }
138 
139 /* Initialise columns k+1:n to columns of the unit matrix */
140 
141  i__1 = *n;
142  for (j = *k + 1; j <= i__1; ++j) {
143  i__2 = *m;
144  for (l = 1; l <= i__2; ++l) {
145  a_ref(l, j) = 0.;
146 /* L10: */
147  }
148  a_ref(j, j) = 1.;
149 /* L20: */
150  }
151 
152  for (i__ = *k; i__ >= 1; --i__) {
153 
154 /* Apply H(i) to A(i:m,i:n) from the left */
155 
156  if (i__ < *n) {
157  a_ref(i__, i__) = 1.;
158  i__1 = *m - i__ + 1;
159  i__2 = *n - i__;
160  template_lapack_larf("Left", &i__1, &i__2, &a_ref(i__, i__), &c__1, &tau[i__], &
161  a_ref(i__, i__ + 1), lda, &work[1]);
162  }
163  if (i__ < *m) {
164  i__1 = *m - i__;
165  d__1 = -tau[i__];
166  template_blas_scal(&i__1, &d__1, &a_ref(i__ + 1, i__), &c__1);
167  }
168  a_ref(i__, i__) = 1. - tau[i__];
169 
170 /* Set A(1:i-1,i) to zero */
171 
172  i__1 = i__ - 1;
173  for (l = 1; l <= i__1; ++l) {
174  a_ref(l, i__) = 0.;
175 /* L30: */
176  }
177 /* L40: */
178  }
179  return 0;
180 
181 /* End of DORG2R */
182 
183 } /* dorg2r_ */
184 
185 #undef a_ref
186 
187 
188 #endif
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition: template_blas_scal.h:41
int integer
Definition: template_blas_common.h:38
#define maxMACRO(a, b)
Definition: template_blas_common.h:43
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
int template_lapack_larf(const char *side, const integer *m, const integer *n, const Treal *v, const integer *incv, const Treal *tau, Treal *c__, const integer *ldc, Treal *work)
Definition: template_lapack_larf.h:40
#define a_ref(a_1, a_2)
int template_lapack_org2r(const integer *m, const integer *n, const integer *k, Treal *a, const integer *lda, const Treal *tau, Treal *work, integer *info)
Definition: template_lapack_org2r.h:40