ergo
template_lapack_getf2.h
Go to the documentation of this file.
1 /* Ergo, version 3.2, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2012 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_GETF2_HEADER
36 #define TEMPLATE_LAPACK_GETF2_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_getf2(const integer *m, const integer *n, Treal *a, const integer *
41  lda, integer *ipiv, 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  June 30, 1992
47 
48 
49  Purpose
50  =======
51 
52  DGETF2 computes an LU factorization of a general m-by-n matrix A
53  using partial pivoting with row interchanges.
54 
55  The factorization has the form
56  A = P * L * U
57  where P is a permutation matrix, L is lower triangular with unit
58  diagonal elements (lower trapezoidal if m > n), and U is upper
59  triangular (upper trapezoidal if m < n).
60 
61  This is the right-looking Level 2 BLAS version of the algorithm.
62 
63  Arguments
64  =========
65 
66  M (input) INTEGER
67  The number of rows of the matrix A. M >= 0.
68 
69  N (input) INTEGER
70  The number of columns of the matrix A. N >= 0.
71 
72  A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
73  On entry, the m by n matrix to be factored.
74  On exit, the factors L and U from the factorization
75  A = P*L*U; the unit diagonal elements of L are not stored.
76 
77  LDA (input) INTEGER
78  The leading dimension of the array A. LDA >= max(1,M).
79 
80  IPIV (output) INTEGER array, dimension (min(M,N))
81  The pivot indices; for 1 <= i <= min(M,N), row i of the
82  matrix was interchanged with row IPIV(i).
83 
84  INFO (output) INTEGER
85  = 0: successful exit
86  < 0: if INFO = -k, the k-th argument had an illegal value
87  > 0: if INFO = k, U(k,k) is exactly zero. The factorization
88  has been completed, but the factor U is exactly
89  singular, and division by zero will occur if it is used
90  to solve a system of equations.
91 
92  =====================================================================
93 
94 
95  Test the input parameters.
96 
97  Parameter adjustments */
98  /* Table of constant values */
99  integer c__1 = 1;
100  Treal c_b6 = -1.;
101 
102  /* System generated locals */
103  integer a_dim1, a_offset, i__1, i__2, i__3;
104  Treal d__1;
105  /* Local variables */
106  integer j;
107  integer jp;
108 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
109 
110 
111  a_dim1 = *lda;
112  a_offset = 1 + a_dim1 * 1;
113  a -= a_offset;
114  --ipiv;
115 
116  /* Function Body */
117  *info = 0;
118  if (*m < 0) {
119  *info = -1;
120  } else if (*n < 0) {
121  *info = -2;
122  } else if (*lda < maxMACRO(1,*m)) {
123  *info = -4;
124  }
125  if (*info != 0) {
126  i__1 = -(*info);
127  template_blas_erbla("GETF2 ", &i__1);
128  return 0;
129  }
130 
131 /* Quick return if possible */
132 
133  if (*m == 0 || *n == 0) {
134  return 0;
135  }
136 
137  i__1 = minMACRO(*m,*n);
138  for (j = 1; j <= i__1; ++j) {
139 
140 /* Find pivot and test for singularity. */
141 
142  i__2 = *m - j + 1;
143  jp = j - 1 + template_blas_idamax(&i__2, &a_ref(j, j), &c__1);
144  ipiv[j] = jp;
145  if (a_ref(jp, j) != 0.) {
146 
147 /* Apply the interchange to columns 1:N. */
148 
149  if (jp != j) {
150  template_blas_swap(n, &a_ref(j, 1), lda, &a_ref(jp, 1), lda);
151  }
152 
153 /* Compute elements J+1:M of J-th column. */
154 
155  if (j < *m) {
156  i__2 = *m - j;
157  d__1 = 1. / a_ref(j, j);
158  template_blas_scal(&i__2, &d__1, &a_ref(j + 1, j), &c__1);
159  }
160 
161  } else if (*info == 0) {
162 
163  *info = j;
164  }
165 
166  if (j < minMACRO(*m,*n)) {
167 
168 /* Update trailing submatrix. */
169 
170  i__2 = *m - j;
171  i__3 = *n - j;
172  template_blas_ger(&i__2, &i__3, &c_b6, &a_ref(j + 1, j), &c__1, &a_ref(j, j +
173  1), lda, &a_ref(j + 1, j + 1), lda);
174  }
175 /* L10: */
176  }
177  return 0;
178 
179 /* End of DGETF2 */
180 
181 } /* dgetf2_ */
182 
183 #undef a_ref
184 
185 
186 #endif