template_lapack_getrf.h

Go to the documentation of this file.
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

Generated on Wed Nov 21 09:32:01 2012 for ergo by  doxygen 1.4.7