template_lapack_potf2.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_POTF2_HEADER
00036 #define TEMPLATE_LAPACK_POTF2_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_potf2(const char *uplo, const integer *n, Treal *a, const integer *
00041         lda, 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        February 29, 1992   
00047 
00048 
00049     Purpose   
00050     =======   
00051 
00052     DPOTF2 computes the Cholesky factorization of a real symmetric   
00053     positive definite matrix A.   
00054 
00055     The factorization has the form   
00056        A = U' * U ,  if UPLO = 'U', or   
00057        A = L  * L',  if UPLO = 'L',   
00058     where U is an upper triangular matrix and L is lower triangular.   
00059 
00060     This is the unblocked version of the algorithm, calling Level 2 BLAS.   
00061 
00062     Arguments   
00063     =========   
00064 
00065     UPLO    (input) CHARACTER*1   
00066             Specifies whether the upper or lower triangular part of the   
00067             symmetric matrix A is stored.   
00068             = 'U':  Upper triangular   
00069             = 'L':  Lower triangular   
00070 
00071     N       (input) INTEGER   
00072             The order of the matrix A.  N >= 0.   
00073 
00074     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)   
00075             On entry, the symmetric matrix A.  If UPLO = 'U', the leading   
00076             n by n upper triangular part of A contains the upper   
00077             triangular part of the matrix A, and the strictly lower   
00078             triangular part of A is not referenced.  If UPLO = 'L', the   
00079             leading n by n lower triangular part of A contains the lower   
00080             triangular part of the matrix A, and the strictly upper   
00081             triangular part of A is not referenced.   
00082 
00083             On exit, if INFO = 0, the factor U or L from the Cholesky   
00084             factorization A = U'*U  or A = L*L'.   
00085 
00086     LDA     (input) INTEGER   
00087             The leading dimension of the array A.  LDA >= max(1,N).   
00088 
00089     INFO    (output) INTEGER   
00090             = 0: successful exit   
00091             < 0: if INFO = -k, the k-th argument had an illegal value   
00092             > 0: if INFO = k, the leading minor of order k is not   
00093                  positive definite, and the factorization could not be   
00094                  completed.   
00095 
00096     =====================================================================   
00097 
00098 
00099        Test the input parameters.   
00100 
00101        Parameter adjustments */
00102     /* Table of constant values */
00103      integer c__1 = 1;
00104      Treal c_b10 = -1.;
00105      Treal c_b12 = 1.;
00106     
00107     /* System generated locals */
00108     integer a_dim1, a_offset, i__1, i__2, i__3;
00109     Treal d__1;
00110     /* Local variables */
00111      integer j;
00112      logical upper;
00113      Treal ajj;
00114 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00115 
00116 
00117     a_dim1 = *lda;
00118     a_offset = 1 + a_dim1 * 1;
00119     a -= a_offset;
00120 
00121     /* Function Body */
00122     *info = 0;
00123     upper = template_blas_lsame(uplo, "U");
00124     if (! upper && ! template_blas_lsame(uplo, "L")) {
00125         *info = -1;
00126     } else if (*n < 0) {
00127         *info = -2;
00128     } else if (*lda < maxMACRO(1,*n)) {
00129         *info = -4;
00130     }
00131     if (*info != 0) {
00132         i__1 = -(*info);
00133         template_blas_erbla("POTF2 ", &i__1);
00134         return 0;
00135     }
00136 
00137 /*     Quick return if possible */
00138 
00139     if (*n == 0) {
00140         return 0;
00141     }
00142 
00143     if (upper) {
00144 
00145 /*        Compute the Cholesky factorization A = U'*U. */
00146 
00147         i__1 = *n;
00148         for (j = 1; j <= i__1; ++j) {
00149 
00150 /*           Compute U(J,J) and test for non-positive-definiteness. */
00151 
00152             i__2 = j - 1;
00153             ajj = a_ref(j, j) - template_blas_dot(&i__2, &a_ref(1, j), &c__1, &a_ref(1, j)
00154                     , &c__1);
00155             if (ajj <= 0.) {
00156                 a_ref(j, j) = ajj;
00157                 goto L30;
00158             }
00159             ajj = template_blas_sqrt(ajj);
00160             a_ref(j, j) = ajj;
00161 
00162 /*           Compute elements J+1:N of row J. */
00163 
00164             if (j < *n) {
00165                 i__2 = j - 1;
00166                 i__3 = *n - j;
00167                 template_blas_gemv("Transpose", &i__2, &i__3, &c_b10, &a_ref(1, j + 1), 
00168                         lda, &a_ref(1, j), &c__1, &c_b12, &a_ref(j, j + 1), 
00169                         lda);
00170                 i__2 = *n - j;
00171                 d__1 = 1. / ajj;
00172                 template_blas_scal(&i__2, &d__1, &a_ref(j, j + 1), lda);
00173             }
00174 /* L10: */
00175         }
00176     } else {
00177 
00178 /*        Compute the Cholesky factorization A = L*L'. */
00179 
00180         i__1 = *n;
00181         for (j = 1; j <= i__1; ++j) {
00182 
00183 /*           Compute L(J,J) and test for non-positive-definiteness. */
00184 
00185             i__2 = j - 1;
00186             ajj = a_ref(j, j) - template_blas_dot(&i__2, &a_ref(j, 1), lda, &a_ref(j, 1), 
00187                     lda);
00188             if (ajj <= 0.) {
00189                 a_ref(j, j) = ajj;
00190                 goto L30;
00191             }
00192             ajj = template_blas_sqrt(ajj);
00193             a_ref(j, j) = ajj;
00194 
00195 /*           Compute elements J+1:N of column J. */
00196 
00197             if (j < *n) {
00198                 i__2 = *n - j;
00199                 i__3 = j - 1;
00200                 template_blas_gemv("No transpose", &i__2, &i__3, &c_b10, &a_ref(j + 1, 1),
00201                          lda, &a_ref(j, 1), lda, &c_b12, &a_ref(j + 1, j), &
00202                         c__1);
00203                 i__2 = *n - j;
00204                 d__1 = 1. / ajj;
00205                 template_blas_scal(&i__2, &d__1, &a_ref(j + 1, j), &c__1);
00206             }
00207 /* L20: */
00208         }
00209     }
00210     goto L40;
00211 
00212 L30:
00213     *info = j;
00214 
00215 L40:
00216     return 0;
00217 
00218 /*     End of DPOTF2 */
00219 
00220 } /* dpotf2_ */
00221 
00222 #undef a_ref
00223 
00224 
00225 #endif

Generated on Mon Sep 17 14:30:40 2012 for ergo by  doxygen 1.4.7