ergo
template_lapack_getrf.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_GETRF_HEADER
38#define TEMPLATE_LAPACK_GETRF_HEADER
39
40
41template<class Treal>
42int template_lapack_getrf(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 March 31, 1993
49
50
51 Purpose
52 =======
53
54 DGETRF 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 3 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 = -i, the i-th argument had an illegal value
89 > 0: if INFO = i, U(i,i) 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 integer c_n1 = -1;
103 Treal c_b16 = 1.;
104 Treal c_b19 = -1.;
105
106 /* System generated locals */
107 integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
108 /* Local variables */
109 integer i__, j;
110 integer iinfo;
111 integer jb, nb;
112#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
113
114
115 a_dim1 = *lda;
116 a_offset = 1 + a_dim1 * 1;
117 a -= a_offset;
118 --ipiv;
119
120 /* Function Body */
121 *info = 0;
122 if (*m < 0) {
123 *info = -1;
124 } else if (*n < 0) {
125 *info = -2;
126 } else if (*lda < maxMACRO(1,*m)) {
127 *info = -4;
128 }
129 if (*info != 0) {
130 i__1 = -(*info);
131 template_blas_erbla("GETRF ", &i__1);
132 return 0;
133 }
134
135/* Quick return if possible */
136
137 if (*m == 0 || *n == 0) {
138 return 0;
139 }
140
141/* Determine the block size for this environment. */
142
143 nb = template_lapack_ilaenv(&c__1, "DGETRF", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)
144 1);
145 if (nb <= 1 || nb >= minMACRO(*m,*n)) {
146
147/* Use unblocked code. */
148
149 template_lapack_getf2(m, n, &a[a_offset], lda, &ipiv[1], info);
150 } else {
151
152/* Use blocked code. */
153
154 i__1 = minMACRO(*m,*n);
155 i__2 = nb;
156 for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
157/* Computing MIN */
158 i__3 = minMACRO(*m,*n) - j + 1;
159 jb = minMACRO(i__3,nb);
160
161/* Factor diagonal and subdiagonal blocks and test for exact
162 singularity. */
163
164 i__3 = *m - j + 1;
165 template_lapack_getf2(&i__3, &jb, &a_ref(j, j), lda, &ipiv[j], &iinfo);
166
167/* Adjust INFO and the pivot indices. */
168
169 if (*info == 0 && iinfo > 0) {
170 *info = iinfo + j - 1;
171 }
172/* Computing MIN */
173 i__4 = *m, i__5 = j + jb - 1;
174 i__3 = minMACRO(i__4,i__5);
175 for (i__ = j; i__ <= i__3; ++i__) {
176 ipiv[i__] = j - 1 + ipiv[i__];
177/* L10: */
178 }
179
180/* Apply interchanges to columns 1:J-1. */
181
182 i__3 = j - 1;
183 i__4 = j + jb - 1;
184 template_lapack_laswp(&i__3, &a[a_offset], lda, &j, &i__4, &ipiv[1], &c__1);
185
186 if (j + jb <= *n) {
187
188/* Apply interchanges to columns J+JB:N. */
189
190 i__3 = *n - j - jb + 1;
191 i__4 = j + jb - 1;
192 template_lapack_laswp(&i__3, &a_ref(1, j + jb), lda, &j, &i__4, &ipiv[1], &
193 c__1);
194
195/* Compute block row of U. */
196
197 i__3 = *n - j - jb + 1;
198 template_blas_trsm("Left", "Lower", "No transpose", "Unit", &jb, &i__3, &
199 c_b16, &a_ref(j, j), lda, &a_ref(j, j + jb), lda);
200 if (j + jb <= *m) {
201
202/* Update trailing submatrix. */
203
204 i__3 = *m - j - jb + 1;
205 i__4 = *n - j - jb + 1;
206 template_blas_gemm("No transpose", "No transpose", &i__3, &i__4, &jb,
207 &c_b19, &a_ref(j + jb, j), lda, &a_ref(j, j + jb),
208 lda, &c_b16, &a_ref(j + jb, j + jb), lda);
209 }
210 }
211/* L20: */
212 }
213 }
214 return 0;
215
216/* End of DGETRF */
217
218} /* dgetrf_ */
219
220#undef a_ref
221
222
223#endif
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
int integer
Definition: template_blas_common.h:40
int ftnlen
Definition: template_blas_common.h:42
#define minMACRO(a, b)
Definition: template_blas_common.h:46
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
int template_blas_gemm(const char *transa, const char *transb, const integer *m, const integer *n, const integer *k, const Treal *alpha, const Treal *a, const integer *lda, const Treal *b, const integer *ldb, const Treal *beta, Treal *c__, const integer *ldc)
Definition: template_blas_gemm.h:43
int template_blas_trsm(const char *side, const char *uplo, const char *transa, const char *diag, const integer *m, const integer *n, const Treal *alpha, const Treal *a, const integer *lda, Treal *b, const integer *ldb)
Definition: template_blas_trsm.h:43
integer template_lapack_ilaenv(const integer *ispec, const char *name__, const char *opts, const integer *n1, const integer *n2, const integer *n3, const integer *n4, ftnlen name_len, ftnlen opts_len)
Definition: template_lapack_common.cc:281
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
#define a_ref(a_1, a_2)
int template_lapack_getrf(const integer *m, const integer *n, Treal *a, const integer *lda, integer *ipiv, integer *info)
Definition: template_lapack_getrf.h:42
int template_lapack_laswp(const integer *n, Treal *a, const integer *lda, const integer *k1, const integer *k2, const integer *ipiv, const integer *incx)
Definition: template_lapack_laswp.h:42