ergo
template_lapack_orm2r.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_ORM2R_HEADER
36 #define TEMPLATE_LAPACK_ORM2R_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_orm2r(const char *side, const char *trans, const integer *m, const integer *n,
41  const integer *k, Treal *a, const integer *lda, const Treal *tau, Treal *
42  c__, const integer *ldc, Treal *work, integer *info)
43 {
44 /* -- LAPACK routine (version 3.0) --
45  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
46  Courant Institute, Argonne National Lab, and Rice University
47  February 29, 1992
48 
49 
50  Purpose
51  =======
52 
53  DORM2R overwrites the general real m by n matrix C with
54 
55  Q * C if SIDE = 'L' and TRANS = 'N', or
56 
57  Q'* C if SIDE = 'L' and TRANS = 'T', or
58 
59  C * Q if SIDE = 'R' and TRANS = 'N', or
60 
61  C * Q' if SIDE = 'R' and TRANS = 'T',
62 
63  where Q is a real orthogonal matrix defined as the product of k
64  elementary reflectors
65 
66  Q = H(1) H(2) . . . H(k)
67 
68  as returned by DGEQRF. Q is of order m if SIDE = 'L' and of order n
69  if SIDE = 'R'.
70 
71  Arguments
72  =========
73 
74  SIDE (input) CHARACTER*1
75  = 'L': apply Q or Q' from the Left
76  = 'R': apply Q or Q' from the Right
77 
78  TRANS (input) CHARACTER*1
79  = 'N': apply Q (No transpose)
80  = 'T': apply Q' (Transpose)
81 
82  M (input) INTEGER
83  The number of rows of the matrix C. M >= 0.
84 
85  N (input) INTEGER
86  The number of columns of the matrix C. N >= 0.
87 
88  K (input) INTEGER
89  The number of elementary reflectors whose product defines
90  the matrix Q.
91  If SIDE = 'L', M >= K >= 0;
92  if SIDE = 'R', N >= K >= 0.
93 
94  A (input) DOUBLE PRECISION array, dimension (LDA,K)
95  The i-th column must contain the vector which defines the
96  elementary reflector H(i), for i = 1,2,...,k, as returned by
97  DGEQRF in the first k columns of its array argument A.
98  A is modified by the routine but restored on exit.
99 
100  LDA (input) INTEGER
101  The leading dimension of the array A.
102  If SIDE = 'L', LDA >= max(1,M);
103  if SIDE = 'R', LDA >= max(1,N).
104 
105  TAU (input) DOUBLE PRECISION array, dimension (K)
106  TAU(i) must contain the scalar factor of the elementary
107  reflector H(i), as returned by DGEQRF.
108 
109  C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
110  On entry, the m by n matrix C.
111  On exit, C is overwritten by Q*C or Q'*C or C*Q' or C*Q.
112 
113  LDC (input) INTEGER
114  The leading dimension of the array C. LDC >= max(1,M).
115 
116  WORK (workspace) DOUBLE PRECISION array, dimension
117  (N) if SIDE = 'L',
118  (M) if SIDE = 'R'
119 
120  INFO (output) INTEGER
121  = 0: successful exit
122  < 0: if INFO = -i, the i-th argument had an illegal value
123 
124  =====================================================================
125 
126 
127  Test the input arguments
128 
129  Parameter adjustments */
130  /* Table of constant values */
131  integer c__1 = 1;
132 
133  /* System generated locals */
134  integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2;
135  /* Local variables */
136  logical left;
137  integer i__;
138  integer i1, i2, i3, ic, jc, mi, ni, nq;
139  logical notran;
140  Treal aii;
141 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
142 #define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1]
143 
144 
145  a_dim1 = *lda;
146  a_offset = 1 + a_dim1 * 1;
147  a -= a_offset;
148  --tau;
149  c_dim1 = *ldc;
150  c_offset = 1 + c_dim1 * 1;
151  c__ -= c_offset;
152  --work;
153 
154  /* Function Body */
155  *info = 0;
156  left = template_blas_lsame(side, "L");
157  notran = template_blas_lsame(trans, "N");
158 
159 /* NQ is the order of Q */
160 
161  if (left) {
162  nq = *m;
163  } else {
164  nq = *n;
165  }
166  if (! left && ! template_blas_lsame(side, "R")) {
167  *info = -1;
168  } else if (! notran && ! template_blas_lsame(trans, "T")) {
169  *info = -2;
170  } else if (*m < 0) {
171  *info = -3;
172  } else if (*n < 0) {
173  *info = -4;
174  } else if (*k < 0 || *k > nq) {
175  *info = -5;
176  } else if (*lda < maxMACRO(1,nq)) {
177  *info = -7;
178  } else if (*ldc < maxMACRO(1,*m)) {
179  *info = -10;
180  }
181  if (*info != 0) {
182  i__1 = -(*info);
183  template_blas_erbla("ORM2R ", &i__1);
184  return 0;
185  }
186 
187 /* Quick return if possible */
188 
189  if (*m == 0 || *n == 0 || *k == 0) {
190  return 0;
191  }
192 
193  if ( ( left && ! notran ) || ( ! left && notran ) ) {
194  i1 = 1;
195  i2 = *k;
196  i3 = 1;
197  } else {
198  i1 = *k;
199  i2 = 1;
200  i3 = -1;
201  }
202 
203  if (left) {
204  ni = *n;
205  jc = 1;
206  } else {
207  mi = *m;
208  ic = 1;
209  }
210 
211  i__1 = i2;
212  i__2 = i3;
213  for (i__ = i1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
214  if (left) {
215 
216 /* H(i) is applied to C(i:m,1:n) */
217 
218  mi = *m - i__ + 1;
219  ic = i__;
220  } else {
221 
222 /* H(i) is applied to C(1:m,i:n) */
223 
224  ni = *n - i__ + 1;
225  jc = i__;
226  }
227 
228 /* Apply H(i) */
229 
230  aii = a_ref(i__, i__);
231  a_ref(i__, i__) = 1.;
232  template_lapack_larf(side, &mi, &ni, &a_ref(i__, i__), &c__1, &tau[i__], &c___ref(
233  ic, jc), ldc, &work[1]);
234  a_ref(i__, i__) = aii;
235 /* L10: */
236  }
237  return 0;
238 
239 /* End of DORM2R */
240 
241 } /* dorm2r_ */
242 
243 #undef c___ref
244 #undef a_ref
245 
246 
247 #endif
int integer
Definition: template_blas_common.h:38
#define a_ref(a_1, a_2)
#define maxMACRO(a, b)
Definition: template_blas_common.h:43
#define c___ref(a_1, a_2)
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
int template_lapack_orm2r(const char *side, const char *trans, const integer *m, const integer *n, const integer *k, Treal *a, const integer *lda, const Treal *tau, Treal *c__, const integer *ldc, Treal *work, integer *info)
Definition: template_lapack_orm2r.h:40
bool logical
Definition: template_blas_common.h:39
side
Definition: Matrix.h:73
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
Definition: Matrix.h:73
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:44