ergo
template_lapack_getf2.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 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
41template<class Treal>
42int 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_erbla(const char *srname, integer *info)
Definition template_blas_common.cc:146
int integer
Definition template_blas_common.h:40
#define minMACRO(a, b)
Definition template_blas_common.h:46
#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
integer template_blas_idamax(const integer *n, const Treal *dx, const integer *incx)
Definition template_blas_idamax.h:42
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition template_blas_scal.h:43
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