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