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