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