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