ergo
template_lapack_org2l.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_ORG2L_HEADER
36 #define TEMPLATE_LAPACK_ORG2L_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_org2l(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  DORG2L generates an m by n real matrix Q with orthonormal columns,
53  which is defined as the last n columns of a product of k elementary
54  reflectors of order m
55 
56  Q = H(k) . . . H(2) H(1)
57 
58  as returned by DGEQLF.
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 (n-k+i)-th column must contain the vector which
75  defines the elementary reflector H(i), for i = 1,2,...,k, as
76  returned by DGEQLF in the last 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 DGEQLF.
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, i__3;
104  Treal d__1;
105  /* Local variables */
106  integer i__, j, l;
107  integer ii;
108 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
109 
110 
111  a_dim1 = *lda;
112  a_offset = 1 + a_dim1 * 1;
113  a -= a_offset;
114  --tau;
115  --work;
116 
117  /* Function Body */
118  *info = 0;
119  if (*m < 0) {
120  *info = -1;
121  } else if (*n < 0 || *n > *m) {
122  *info = -2;
123  } else if (*k < 0 || *k > *n) {
124  *info = -3;
125  } else if (*lda < maxMACRO(1,*m)) {
126  *info = -5;
127  }
128  if (*info != 0) {
129  i__1 = -(*info);
130  template_blas_erbla("ORG2L ", &i__1);
131  return 0;
132  }
133 
134 /* Quick return if possible */
135 
136  if (*n <= 0) {
137  return 0;
138  }
139 
140 /* Initialise columns 1:n-k to columns of the unit matrix */
141 
142  i__1 = *n - *k;
143  for (j = 1; j <= i__1; ++j) {
144  i__2 = *m;
145  for (l = 1; l <= i__2; ++l) {
146  a_ref(l, j) = 0.;
147 /* L10: */
148  }
149  a_ref(*m - *n + j, j) = 1.;
150 /* L20: */
151  }
152 
153  i__1 = *k;
154  for (i__ = 1; i__ <= i__1; ++i__) {
155  ii = *n - *k + i__;
156 
157 /* Apply H(i) to A(1:m-k+i,1:n-k+i) from the left */
158 
159  a_ref(*m - *n + ii, ii) = 1.;
160  i__2 = *m - *n + ii;
161  i__3 = ii - 1;
162  template_lapack_larf("Left", &i__2, &i__3, &a_ref(1, ii), &c__1, &tau[i__], &a[
163  a_offset], lda, &work[1]);
164  i__2 = *m - *n + ii - 1;
165  d__1 = -tau[i__];
166  template_blas_scal(&i__2, &d__1, &a_ref(1, ii), &c__1);
167  a_ref(*m - *n + ii, ii) = 1. - tau[i__];
168 
169 /* Set A(m-k+i+1:m,n-k+i) to zero */
170 
171  i__2 = *m;
172  for (l = *m - *n + ii + 1; l <= i__2; ++l) {
173  a_ref(l, ii) = 0.;
174 /* L30: */
175  }
176 /* L40: */
177  }
178  return 0;
179 
180 /* End of DORG2L */
181 
182 } /* dorg2l_ */
183 
184 #undef a_ref
185 
186 
187 #endif