ergo
template_lapack_orm2r.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4  * and Anastasia Kruchinina.
5  *
6  * This program is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program. If not, see <http://www.gnu.org/licenses/>.
18  *
19  * Primary academic reference:
20  * Ergo: An open-source program for linear-scaling electronic structure
21  * calculations,
22  * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23  * Kruchinina,
24  * SoftwareX 7, 107 (2018),
25  * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26  *
27  * For further information about Ergo, see <http://www.ergoscf.org>.
28  */
29 
30  /* This file belongs to the template_lapack part of the Ergo source
31  * code. The source files in the template_lapack directory are modified
32  * versions of files originally distributed as CLAPACK, see the
33  * Copyright/license notice in the file template_lapack/COPYING.
34  */
35 
36 
37 #ifndef TEMPLATE_LAPACK_ORM2R_HEADER
38 #define TEMPLATE_LAPACK_ORM2R_HEADER
39 
40 
41 template<class Treal>
42 int template_lapack_orm2r(const char *side, const char *trans, const integer *m, const integer *n,
43  const integer *k, Treal *a, const integer *lda, const Treal *tau, Treal *
44  c__, const integer *ldc, Treal *work, integer *info)
45 {
46 /* -- LAPACK routine (version 3.0) --
47  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
48  Courant Institute, Argonne National Lab, and Rice University
49  February 29, 1992
50 
51 
52  Purpose
53  =======
54 
55  DORM2R overwrites the general real m by n matrix C with
56 
57  Q * C if SIDE = 'L' and TRANS = 'N', or
58 
59  Q'* C if SIDE = 'L' and TRANS = 'T', or
60 
61  C * Q if SIDE = 'R' and TRANS = 'N', or
62 
63  C * Q' if SIDE = 'R' and TRANS = 'T',
64 
65  where Q is a real orthogonal matrix defined as the product of k
66  elementary reflectors
67 
68  Q = H(1) H(2) . . . H(k)
69 
70  as returned by DGEQRF. Q is of order m if SIDE = 'L' and of order n
71  if SIDE = 'R'.
72 
73  Arguments
74  =========
75 
76  SIDE (input) CHARACTER*1
77  = 'L': apply Q or Q' from the Left
78  = 'R': apply Q or Q' from the Right
79 
80  TRANS (input) CHARACTER*1
81  = 'N': apply Q (No transpose)
82  = 'T': apply Q' (Transpose)
83 
84  M (input) INTEGER
85  The number of rows of the matrix C. M >= 0.
86 
87  N (input) INTEGER
88  The number of columns of the matrix C. N >= 0.
89 
90  K (input) INTEGER
91  The number of elementary reflectors whose product defines
92  the matrix Q.
93  If SIDE = 'L', M >= K >= 0;
94  if SIDE = 'R', N >= K >= 0.
95 
96  A (input) DOUBLE PRECISION array, dimension (LDA,K)
97  The i-th column must contain the vector which defines the
98  elementary reflector H(i), for i = 1,2,...,k, as returned by
99  DGEQRF in the first k columns of its array argument A.
100  A is modified by the routine but restored on exit.
101 
102  LDA (input) INTEGER
103  The leading dimension of the array A.
104  If SIDE = 'L', LDA >= max(1,M);
105  if SIDE = 'R', LDA >= max(1,N).
106 
107  TAU (input) DOUBLE PRECISION array, dimension (K)
108  TAU(i) must contain the scalar factor of the elementary
109  reflector H(i), as returned by DGEQRF.
110 
111  C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
112  On entry, the m by n matrix C.
113  On exit, C is overwritten by Q*C or Q'*C or C*Q' or C*Q.
114 
115  LDC (input) INTEGER
116  The leading dimension of the array C. LDC >= max(1,M).
117 
118  WORK (workspace) DOUBLE PRECISION array, dimension
119  (N) if SIDE = 'L',
120  (M) if SIDE = 'R'
121 
122  INFO (output) INTEGER
123  = 0: successful exit
124  < 0: if INFO = -i, the i-th argument had an illegal value
125 
126  =====================================================================
127 
128 
129  Test the input arguments
130 
131  Parameter adjustments */
132  /* Table of constant values */
133  integer c__1 = 1;
134 
135  /* System generated locals */
136  integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2;
137  /* Local variables */
138  logical left;
139  integer i__;
140  integer i1, i2, i3, ic, jc, mi, ni, nq;
141  logical notran;
142  Treal aii;
143 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
144 #define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1]
145 
146 
147  a_dim1 = *lda;
148  a_offset = 1 + a_dim1 * 1;
149  a -= a_offset;
150  --tau;
151  c_dim1 = *ldc;
152  c_offset = 1 + c_dim1 * 1;
153  c__ -= c_offset;
154  --work;
155 
156  /* Function Body */
157  *info = 0;
158  left = template_blas_lsame(side, "L");
159  notran = template_blas_lsame(trans, "N");
160 
161 /* NQ is the order of Q */
162 
163  if (left) {
164  nq = *m;
165  } else {
166  nq = *n;
167  }
168  if (! left && ! template_blas_lsame(side, "R")) {
169  *info = -1;
170  } else if (! notran && ! template_blas_lsame(trans, "T")) {
171  *info = -2;
172  } else if (*m < 0) {
173  *info = -3;
174  } else if (*n < 0) {
175  *info = -4;
176  } else if (*k < 0 || *k > nq) {
177  *info = -5;
178  } else if (*lda < maxMACRO(1,nq)) {
179  *info = -7;
180  } else if (*ldc < maxMACRO(1,*m)) {
181  *info = -10;
182  }
183  if (*info != 0) {
184  i__1 = -(*info);
185  template_blas_erbla("ORM2R ", &i__1);
186  return 0;
187  }
188 
189 /* Quick return if possible */
190 
191  if (*m == 0 || *n == 0 || *k == 0) {
192  return 0;
193  }
194 
195  if ( ( left && ! notran ) || ( ! left && notran ) ) {
196  i1 = 1;
197  i2 = *k;
198  i3 = 1;
199  } else {
200  i1 = *k;
201  i2 = 1;
202  i3 = -1;
203  }
204 
205  if (left) {
206  ni = *n;
207  jc = 1;
208  } else {
209  mi = *m;
210  ic = 1;
211  }
212 
213  i__1 = i2;
214  i__2 = i3;
215  for (i__ = i1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
216  if (left) {
217 
218 /* H(i) is applied to C(i:m,1:n) */
219 
220  mi = *m - i__ + 1;
221  ic = i__;
222  } else {
223 
224 /* H(i) is applied to C(1:m,i:n) */
225 
226  ni = *n - i__ + 1;
227  jc = i__;
228  }
229 
230 /* Apply H(i) */
231 
232  aii = a_ref(i__, i__);
233  a_ref(i__, i__) = 1.;
234  template_lapack_larf(side, &mi, &ni, &a_ref(i__, i__), &c__1, &tau[i__], &c___ref(
235  ic, jc), ldc, &work[1]);
236  a_ref(i__, i__) = aii;
237 /* L10: */
238  }
239  return 0;
240 
241 /* End of DORM2R */
242 
243 } /* dorm2r_ */
244 
245 #undef c___ref
246 #undef a_ref
247 
248 
249 #endif
int integer
Definition: template_blas_common.h:40
#define a_ref(a_1, a_2)
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
#define c___ref(a_1, a_2)
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
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:42
bool logical
Definition: template_blas_common.h:41
side
Definition: Matrix.h:75
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:42
Definition: Matrix.h:75
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46