template_lapack_pptrf.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_PPTRF_HEADER
00036 #define TEMPLATE_LAPACK_PPTRF_HEADER
00037 
00038 #include "template_lapack_common.h"
00039 
00040 template<class Treal>
00041 int template_lapack_pptrf(const char *uplo, const integer *n, Treal *ap, integer *
00042         info)
00043 {
00044 /*  -- LAPACK routine (version 3.0) --   
00045        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
00046        Courant Institute, Argonne National Lab, and Rice University   
00047        March 31, 1993   
00048 
00049 
00050     Purpose   
00051     =======   
00052 
00053     DPPTRF computes the Cholesky factorization of a real symmetric   
00054     positive definite matrix A stored in packed format.   
00055 
00056     The factorization has the form   
00057        A = U**T * U,  if UPLO = 'U', or   
00058        A = L  * L**T,  if UPLO = 'L',   
00059     where U is an upper triangular matrix and L is lower triangular.   
00060 
00061     Arguments   
00062     =========   
00063 
00064     UPLO    (input) CHARACTER*1   
00065             = 'U':  Upper triangle of A is stored;   
00066             = 'L':  Lower triangle of A is stored.   
00067 
00068     N       (input) INTEGER   
00069             The order of the matrix A.  N >= 0.   
00070 
00071     AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)   
00072             On entry, the upper or lower triangle of the symmetric matrix   
00073             A, packed columnwise in a linear array.  The j-th column of A   
00074             is stored in the array AP as follows:   
00075             if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;   
00076             if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.   
00077             See below for further details.   
00078 
00079             On exit, if INFO = 0, the triangular factor U or L from the   
00080             Cholesky factorization A = U**T*U or A = L*L**T, in the same   
00081             storage format as A.   
00082 
00083     INFO    (output) INTEGER   
00084             = 0:  successful exit   
00085             < 0:  if INFO = -i, the i-th argument had an illegal value   
00086             > 0:  if INFO = i, the leading minor of order i is not   
00087                   positive definite, and the factorization could not be   
00088                   completed.   
00089 
00090     Further Details   
00091     ======= =======   
00092 
00093     The packed storage scheme is illustrated by the following example   
00094     when N = 4, UPLO = 'U':   
00095 
00096     Two-dimensional storage of the symmetric matrix A:   
00097 
00098        a11 a12 a13 a14   
00099            a22 a23 a24   
00100                a33 a34     (aij = aji)   
00101                    a44   
00102 
00103     Packed storage of the upper triangle of A:   
00104 
00105     AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]   
00106 
00107     =====================================================================   
00108 
00109 
00110        Test the input parameters.   
00111 
00112        Parameter adjustments */
00113     /* Table of constant values */
00114      integer c__1 = 1;
00115      Treal c_b16 = -1.;
00116     
00117     /* System generated locals */
00118     integer i__1, i__2;
00119     Treal d__1;
00120     /* Local variables */
00121      integer j;
00122      logical upper;
00123      integer jc, jj;
00124      Treal ajj;
00125 
00126 
00127     --ap;
00128 
00129     /* Function Body */
00130     *info = 0;
00131     upper = template_blas_lsame(uplo, "U");
00132     if (! upper && ! template_blas_lsame(uplo, "L")) {
00133         *info = -1;
00134     } else if (*n < 0) {
00135         *info = -2;
00136     }
00137     if (*info != 0) {
00138         i__1 = -(*info);
00139         template_blas_erbla("DPPTRF", &i__1);
00140         return 0;
00141     }
00142 
00143 /*     Quick return if possible */
00144 
00145     if (*n == 0) {
00146         return 0;
00147     }
00148 
00149     if (upper) {
00150 
00151 /*        Compute the Cholesky factorization A = U'*U. */
00152 
00153         jj = 0;
00154         i__1 = *n;
00155         for (j = 1; j <= i__1; ++j) {
00156             jc = jj + 1;
00157             jj += j;
00158 
00159 /*           Compute elements 1:J-1 of column J. */
00160 
00161             if (j > 1) {
00162                 i__2 = j - 1;
00163                 template_blas_tpsv("Upper", "Transpose", "Non-unit", &i__2, &ap[1], &ap[
00164                         jc], &c__1);
00165             }
00166 
00167 /*           Compute U(J,J) and test for non-positive-definiteness. */
00168 
00169             i__2 = j - 1;
00170             ajj = ap[jj] - template_blas_dot(&i__2, &ap[jc], &c__1, &ap[jc], &c__1);
00171             if (ajj <= 0.) {
00172                 ap[jj] = ajj;
00173                 goto L30;
00174             }
00175             ap[jj] = template_blas_sqrt(ajj);
00176 /* L10: */
00177         }
00178     } else {
00179 
00180 /*        Compute the Cholesky factorization A = L*L'. */
00181 
00182         jj = 1;
00183         i__1 = *n;
00184         for (j = 1; j <= i__1; ++j) {
00185 
00186 /*           Compute L(J,J) and test for non-positive-definiteness. */
00187 
00188             ajj = ap[jj];
00189             if (ajj <= 0.) {
00190                 ap[jj] = ajj;
00191                 goto L30;
00192             }
00193             ajj = template_blas_sqrt(ajj);
00194             ap[jj] = ajj;
00195 
00196 /*           Compute elements J+1:N of column J and update the trailing   
00197              submatrix. */
00198 
00199             if (j < *n) {
00200                 i__2 = *n - j;
00201                 d__1 = 1. / ajj;
00202                 template_blas_scal(&i__2, &d__1, &ap[jj + 1], &c__1);
00203                 i__2 = *n - j;
00204                 template_blas_spr("Lower", &i__2, &c_b16, &ap[jj + 1], &c__1, &ap[jj + *n 
00205                         - j + 1]);
00206                 jj = jj + *n - j + 1;
00207             }
00208 /* L20: */
00209         }
00210     }
00211     goto L40;
00212 
00213 L30:
00214     *info = j;
00215 
00216 L40:
00217     return 0;
00218 
00219 /*     End of DPPTRF */
00220 
00221 } /* dpptrf_ */
00222 
00223 #endif

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