ergo
template_lapack_getrf.h
Go to the documentation of this file.
1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2019 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 
41 template<class Treal>
42 int 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
ftnlen
int ftnlen
Definition: template_blas_common.h:42
template_lapack_ilaenv
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
minMACRO
#define minMACRO(a, b)
Definition: template_blas_common.h:46
template_blas_trsm
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
template_blas_gemm
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
template_lapack_getf2
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
template_lapack_laswp
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
a_ref
#define a_ref(a_1, a_2)
template_blas_erbla
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
integer
int integer
Definition: template_blas_common.h:40
template_lapack_getrf
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
maxMACRO
#define maxMACRO(a, b)
Definition: template_blas_common.h:45