ergo
template_lapack_getf2.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_GETF2_HEADER
00036 #define TEMPLATE_LAPACK_GETF2_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_getf2(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        June 30, 1992   
00047 
00048 
00049     Purpose   
00050     =======   
00051 
00052     DGETF2 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 2 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 = -k, the k-th argument had an illegal value   
00087             > 0: if INFO = k, U(k,k) 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      Treal c_b6 = -1.;
00101     
00102     /* System generated locals */
00103     integer a_dim1, a_offset, i__1, i__2, i__3;
00104     Treal d__1;
00105     /* Local variables */
00106      integer j;
00107      integer jp;
00108 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00109 
00110 
00111     a_dim1 = *lda;
00112     a_offset = 1 + a_dim1 * 1;
00113     a -= a_offset;
00114     --ipiv;
00115 
00116     /* Function Body */
00117     *info = 0;
00118     if (*m < 0) {
00119         *info = -1;
00120     } else if (*n < 0) {
00121         *info = -2;
00122     } else if (*lda < maxMACRO(1,*m)) {
00123         *info = -4;
00124     }
00125     if (*info != 0) {
00126         i__1 = -(*info);
00127         template_blas_erbla("GETF2 ", &i__1);
00128         return 0;
00129     }
00130 
00131 /*     Quick return if possible */
00132 
00133     if (*m == 0 || *n == 0) {
00134         return 0;
00135     }
00136 
00137     i__1 = minMACRO(*m,*n);
00138     for (j = 1; j <= i__1; ++j) {
00139 
00140 /*        Find pivot and test for singularity. */
00141 
00142         i__2 = *m - j + 1;
00143         jp = j - 1 + template_blas_idamax(&i__2, &a_ref(j, j), &c__1);
00144         ipiv[j] = jp;
00145         if (a_ref(jp, j) != 0.) {
00146 
00147 /*           Apply the interchange to columns 1:N. */
00148 
00149             if (jp != j) {
00150               template_blas_swap(n, &a_ref(j, 1), lda, &a_ref(jp, 1), lda);
00151             }
00152 
00153 /*           Compute elements J+1:M of J-th column. */
00154 
00155             if (j < *m) {
00156                 i__2 = *m - j;
00157                 d__1 = 1. / a_ref(j, j);
00158                 template_blas_scal(&i__2, &d__1, &a_ref(j + 1, j), &c__1);
00159             }
00160 
00161         } else if (*info == 0) {
00162 
00163             *info = j;
00164         }
00165 
00166         if (j < minMACRO(*m,*n)) {
00167 
00168 /*           Update trailing submatrix. */
00169 
00170             i__2 = *m - j;
00171             i__3 = *n - j;
00172             template_blas_ger(&i__2, &i__3, &c_b6, &a_ref(j + 1, j), &c__1, &a_ref(j, j + 
00173                     1), lda, &a_ref(j + 1, j + 1), lda);
00174         }
00175 /* L10: */
00176     }
00177     return 0;
00178 
00179 /*     End of DGETF2 */
00180 
00181 } /* dgetf2_ */
00182 
00183 #undef a_ref
00184 
00185 
00186 #endif