00001 /* Ergo, version 3.7, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2018 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, 00004 * and Anastasia Kruchinina. 00005 * 00006 * This program is free software: you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation, either version 3 of the License, or 00009 * (at your option) any later version. 00010 * 00011 * This program is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00018 * 00019 * Primary academic reference: 00020 * Ergo: An open-source program for linear-scaling electronic structure 00021 * calculations, 00022 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia 00023 * Kruchinina, 00024 * SoftwareX 7, 107 (2018), 00025 * <http://dx.doi.org/10.1016/j.softx.2018.03.005> 00026 * 00027 * For further information about Ergo, see <http://www.ergoscf.org>. 00028 */ 00029 00030 /* This file belongs to the template_lapack part of the Ergo source 00031 * code. The source files in the template_lapack directory are modified 00032 * versions of files originally distributed as CLAPACK, see the 00033 * Copyright/license notice in the file template_lapack/COPYING. 00034 */ 00035 00036 00037 #ifndef TEMPLATE_LAPACK_GETRF_HEADER 00038 #define TEMPLATE_LAPACK_GETRF_HEADER 00039 00040 00041 template<class Treal> 00042 int template_lapack_getrf(const integer *m, const integer *n, Treal *a, const integer * 00043 lda, integer *ipiv, integer *info) 00044 { 00045 /* -- LAPACK routine (version 3.0) -- 00046 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00047 Courant Institute, Argonne National Lab, and Rice University 00048 March 31, 1993 00049 00050 00051 Purpose 00052 ======= 00053 00054 DGETRF computes an LU factorization of a general M-by-N matrix A 00055 using partial pivoting with row interchanges. 00056 00057 The factorization has the form 00058 A = P * L * U 00059 where P is a permutation matrix, L is lower triangular with unit 00060 diagonal elements (lower trapezoidal if m > n), and U is upper 00061 triangular (upper trapezoidal if m < n). 00062 00063 This is the right-looking Level 3 BLAS version of the algorithm. 00064 00065 Arguments 00066 ========= 00067 00068 M (input) INTEGER 00069 The number of rows of the matrix A. M >= 0. 00070 00071 N (input) INTEGER 00072 The number of columns of the matrix A. N >= 0. 00073 00074 A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00075 On entry, the M-by-N matrix to be factored. 00076 On exit, the factors L and U from the factorization 00077 A = P*L*U; the unit diagonal elements of L are not stored. 00078 00079 LDA (input) INTEGER 00080 The leading dimension of the array A. LDA >= max(1,M). 00081 00082 IPIV (output) INTEGER array, dimension (min(M,N)) 00083 The pivot indices; for 1 <= i <= min(M,N), row i of the 00084 matrix was interchanged with row IPIV(i). 00085 00086 INFO (output) INTEGER 00087 = 0: successful exit 00088 < 0: if INFO = -i, the i-th argument had an illegal value 00089 > 0: if INFO = i, U(i,i) is exactly zero. The factorization 00090 has been completed, but the factor U is exactly 00091 singular, and division by zero will occur if it is used 00092 to solve a system of equations. 00093 00094 ===================================================================== 00095 00096 00097 Test the input parameters. 00098 00099 Parameter adjustments */ 00100 /* Table of constant values */ 00101 integer c__1 = 1; 00102 integer c_n1 = -1; 00103 Treal c_b16 = 1.; 00104 Treal c_b19 = -1.; 00105 00106 /* System generated locals */ 00107 integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5; 00108 /* Local variables */ 00109 integer i__, j; 00110 integer iinfo; 00111 integer jb, nb; 00112 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 00113 00114 00115 a_dim1 = *lda; 00116 a_offset = 1 + a_dim1 * 1; 00117 a -= a_offset; 00118 --ipiv; 00119 00120 /* Function Body */ 00121 *info = 0; 00122 if (*m < 0) { 00123 *info = -1; 00124 } else if (*n < 0) { 00125 *info = -2; 00126 } else if (*lda < maxMACRO(1,*m)) { 00127 *info = -4; 00128 } 00129 if (*info != 0) { 00130 i__1 = -(*info); 00131 template_blas_erbla("GETRF ", &i__1); 00132 return 0; 00133 } 00134 00135 /* Quick return if possible */ 00136 00137 if (*m == 0 || *n == 0) { 00138 return 0; 00139 } 00140 00141 /* Determine the block size for this environment. */ 00142 00143 nb = template_lapack_ilaenv(&c__1, "DGETRF", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen) 00144 1); 00145 if (nb <= 1 || nb >= minMACRO(*m,*n)) { 00146 00147 /* Use unblocked code. */ 00148 00149 template_lapack_getf2(m, n, &a[a_offset], lda, &ipiv[1], info); 00150 } else { 00151 00152 /* Use blocked code. */ 00153 00154 i__1 = minMACRO(*m,*n); 00155 i__2 = nb; 00156 for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) { 00157 /* Computing MIN */ 00158 i__3 = minMACRO(*m,*n) - j + 1; 00159 jb = minMACRO(i__3,nb); 00160 00161 /* Factor diagonal and subdiagonal blocks and test for exact 00162 singularity. */ 00163 00164 i__3 = *m - j + 1; 00165 template_lapack_getf2(&i__3, &jb, &a_ref(j, j), lda, &ipiv[j], &iinfo); 00166 00167 /* Adjust INFO and the pivot indices. */ 00168 00169 if (*info == 0 && iinfo > 0) { 00170 *info = iinfo + j - 1; 00171 } 00172 /* Computing MIN */ 00173 i__4 = *m, i__5 = j + jb - 1; 00174 i__3 = minMACRO(i__4,i__5); 00175 for (i__ = j; i__ <= i__3; ++i__) { 00176 ipiv[i__] = j - 1 + ipiv[i__]; 00177 /* L10: */ 00178 } 00179 00180 /* Apply interchanges to columns 1:J-1. */ 00181 00182 i__3 = j - 1; 00183 i__4 = j + jb - 1; 00184 template_lapack_laswp(&i__3, &a[a_offset], lda, &j, &i__4, &ipiv[1], &c__1); 00185 00186 if (j + jb <= *n) { 00187 00188 /* Apply interchanges to columns J+JB:N. */ 00189 00190 i__3 = *n - j - jb + 1; 00191 i__4 = j + jb - 1; 00192 template_lapack_laswp(&i__3, &a_ref(1, j + jb), lda, &j, &i__4, &ipiv[1], & 00193 c__1); 00194 00195 /* Compute block row of U. */ 00196 00197 i__3 = *n - j - jb + 1; 00198 template_blas_trsm("Left", "Lower", "No transpose", "Unit", &jb, &i__3, & 00199 c_b16, &a_ref(j, j), lda, &a_ref(j, j + jb), lda); 00200 if (j + jb <= *m) { 00201 00202 /* Update trailing submatrix. */ 00203 00204 i__3 = *m - j - jb + 1; 00205 i__4 = *n - j - jb + 1; 00206 template_blas_gemm("No transpose", "No transpose", &i__3, &i__4, &jb, 00207 &c_b19, &a_ref(j + jb, j), lda, &a_ref(j, j + jb), 00208 lda, &c_b16, &a_ref(j + jb, j + jb), lda); 00209 } 00210 } 00211 /* L20: */ 00212 } 00213 } 00214 return 0; 00215 00216 /* End of DGETRF */ 00217 00218 } /* dgetrf_ */ 00219 00220 #undef a_ref 00221 00222 00223 #endif