ergo
template_lapack_geqr2.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_GEQR2_HEADER
38 #define TEMPLATE_LAPACK_GEQR2_HEADER
39 
40 
41 template<class Treal>
42 int template_lapack_geqr2(const integer *m, const integer *n, Treal *a, const integer *
43  lda, 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  DGEQR2 computes a QR factorization of a real m by n matrix A:
55  A = Q * R.
56 
57  Arguments
58  =========
59 
60  M (input) INTEGER
61  The number of rows of the matrix A. M >= 0.
62 
63  N (input) INTEGER
64  The number of columns of the matrix A. N >= 0.
65 
66  A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
67  On entry, the m by n matrix A.
68  On exit, the elements on and above the diagonal of the array
69  contain the min(m,n) by n upper trapezoidal matrix R (R is
70  upper triangular if m >= n); the elements below the diagonal,
71  with the array TAU, represent the orthogonal matrix Q as a
72  product of elementary reflectors (see Further Details).
73 
74  LDA (input) INTEGER
75  The leading dimension of the array A. LDA >= max(1,M).
76 
77  TAU (output) DOUBLE PRECISION array, dimension (min(M,N))
78  The scalar factors of the elementary reflectors (see Further
79  Details).
80 
81  WORK (workspace) DOUBLE PRECISION array, dimension (N)
82 
83  INFO (output) INTEGER
84  = 0: successful exit
85  < 0: if INFO = -i, the i-th argument had an illegal value
86 
87  Further Details
88  ===============
89 
90  The matrix Q is represented as a product of elementary reflectors
91 
92  Q = H(1) H(2) . . . H(k), where k = min(m,n).
93 
94  Each H(i) has the form
95 
96  H(i) = I - tau * v * v'
97 
98  where tau is a real scalar, and v is a real vector with
99  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
100  and tau in TAU(i).
101 
102  =====================================================================
103 
104 
105  Test the input arguments
106 
107  Parameter adjustments */
108  /* Table of constant values */
109  integer c__1 = 1;
110 
111  /* System generated locals */
112  integer a_dim1, a_offset, i__1, i__2, i__3;
113  /* Local variables */
114  integer i__, k;
115  Treal aii;
116 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
117 
118 
119  a_dim1 = *lda;
120  a_offset = 1 + a_dim1 * 1;
121  a -= a_offset;
122  --tau;
123  --work;
124 
125  /* Function Body */
126  *info = 0;
127  if (*m < 0) {
128  *info = -1;
129  } else if (*n < 0) {
130  *info = -2;
131  } else if (*lda < maxMACRO(1,*m)) {
132  *info = -4;
133  }
134  if (*info != 0) {
135  i__1 = -(*info);
136  template_blas_erbla("GEQR2 ", &i__1);
137  return 0;
138  }
139 
140  k = minMACRO(*m,*n);
141 
142  i__1 = k;
143  for (i__ = 1; i__ <= i__1; ++i__) {
144 
145 /* Generate elementary reflector H(i) to annihilate A(i+1:m,i)
146 
147  Computing MIN */
148  i__2 = i__ + 1;
149  i__3 = *m - i__ + 1;
150  template_lapack_larfg(&i__3, &a_ref(i__, i__), &a_ref(minMACRO(i__2,*m), i__), &c__1, &
151  tau[i__]);
152  if (i__ < *n) {
153 
154 /* Apply H(i) to A(i:m,i+1:n) from the left */
155 
156  aii = a_ref(i__, i__);
157  a_ref(i__, i__) = 1.;
158  i__2 = *m - i__ + 1;
159  i__3 = *n - i__;
160  template_lapack_larf("Left", &i__2, &i__3, &a_ref(i__, i__), &c__1, &tau[i__], &
161  a_ref(i__, i__ + 1), lda, &work[1]);
162  a_ref(i__, i__) = aii;
163  }
164 /* L10: */
165  }
166  return 0;
167 
168 /* End of DGEQR2 */
169 
170 } /* dgeqr2_ */
171 
172 #undef a_ref
173 
174 
175 #endif
template_lapack_larfg
int template_lapack_larfg(const integer *n, Treal *alpha, Treal *x, const integer *incx, Treal *tau)
Definition: template_lapack_larfg.h:43
minMACRO
#define minMACRO(a, b)
Definition: template_blas_common.h:46
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
a_ref
#define a_ref(a_1, a_2)
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
template_lapack_geqr2
int template_lapack_geqr2(const integer *m, const integer *n, Treal *a, const integer *lda, Treal *tau, Treal *work, integer *info)
Definition: template_lapack_geqr2.h:42
maxMACRO
#define maxMACRO(a, b)
Definition: template_blas_common.h:45