template_lapack_potrf.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_POTRF_HEADER
00036 #define TEMPLATE_LAPACK_POTRF_HEADER
00037 
00038 #include "template_lapack_potf2.h"
00039 
00040 template<class Treal>
00041 int template_lapack_potrf(const char *uplo, const integer *n, Treal *a, const integer *
00042         lda, integer *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     DPOTRF computes the Cholesky factorization of a real symmetric   
00054     positive definite matrix A.   
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     This is the block version of the algorithm, calling Level 3 BLAS.   
00062 
00063     Arguments   
00064     =========   
00065 
00066     UPLO    (input) CHARACTER*1   
00067             = 'U':  Upper triangle of A is stored;   
00068             = 'L':  Lower triangle of A is stored.   
00069 
00070     N       (input) INTEGER   
00071             The order of the matrix A.  N >= 0.   
00072 
00073     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)   
00074             On entry, the symmetric matrix A.  If UPLO = 'U', the leading   
00075             N-by-N upper triangular part of A contains the upper   
00076             triangular part of the matrix A, and the strictly lower   
00077             triangular part of A is not referenced.  If UPLO = 'L', the   
00078             leading N-by-N lower triangular part of A contains the lower   
00079             triangular part of the matrix A, and the strictly upper   
00080             triangular part of A is not referenced.   
00081 
00082             On exit, if INFO = 0, the factor U or L from the Cholesky   
00083             factorization A = U**T*U or A = L*L**T.   
00084 
00085     LDA     (input) INTEGER   
00086             The leading dimension of the array A.  LDA >= max(1,N).   
00087 
00088     INFO    (output) INTEGER   
00089             = 0:  successful exit   
00090             < 0:  if INFO = -i, the i-th argument had an illegal value   
00091             > 0:  if INFO = i, the leading minor of order i is not   
00092                   positive definite, and the factorization could not be   
00093                   completed.   
00094 
00095     =====================================================================   
00096 
00097 
00098        Test the input parameters.   
00099 
00100        Parameter adjustments */
00101     /* Table of constant values */
00102      integer c__1 = 1;
00103      integer c_n1 = -1;
00104      Treal c_b13 = -1.;
00105      Treal c_b14 = 1.;
00106     
00107     /* System generated locals */
00108     integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
00109     /* Local variables */
00110      integer j;
00111      logical upper;
00112      integer jb, nb;
00113 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00114 
00115 
00116     a_dim1 = *lda;
00117     a_offset = 1 + a_dim1 * 1;
00118     a -= a_offset;
00119 
00120     /* Function Body */
00121     *info = 0;
00122     upper = template_blas_lsame(uplo, "U");
00123     if (! upper && ! template_blas_lsame(uplo, "L")) {
00124         *info = -1;
00125     } else if (*n < 0) {
00126         *info = -2;
00127     } else if (*lda < maxMACRO(1,*n)) {
00128         *info = -4;
00129     }
00130     if (*info != 0) {
00131         i__1 = -(*info);
00132         template_blas_erbla("POTRF ", &i__1);
00133         return 0;
00134     }
00135 
00136 /*     Quick return if possible */
00137 
00138     if (*n == 0) {
00139         return 0;
00140     }
00141 
00142 /*     Determine the block size for this environment. */
00143 
00144     nb = template_lapack_ilaenv(&c__1, "DPOTRF", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (
00145             ftnlen)1);
00146     if (nb <= 1 || nb >= *n) {
00147 
00148 /*        Use unblocked code. */
00149 
00150         template_lapack_potf2(uplo, n, &a[a_offset], lda, info);
00151     } else {
00152 
00153 /*        Use blocked code. */
00154 
00155         if (upper) {
00156 
00157 /*           Compute the Cholesky factorization A = U'*U. */
00158 
00159             i__1 = *n;
00160             i__2 = nb;
00161             for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
00162 
00163 /*              Update and factorize the current diagonal block and test   
00164                 for non-positive-definiteness.   
00165 
00166    Computing MIN */
00167                 i__3 = nb, i__4 = *n - j + 1;
00168                 jb = minMACRO(i__3,i__4);
00169                 i__3 = j - 1;
00170                 template_blas_syrk("Upper", "Transpose", &jb, &i__3, &c_b13, &a_ref(1, j),
00171                          lda, &c_b14, &a_ref(j, j), lda)
00172                         ;
00173                 template_lapack_potf2("Upper", &jb, &a_ref(j, j), lda, info);
00174                 if (*info != 0) {
00175                     goto L30;
00176                 }
00177                 if (j + jb <= *n) {
00178 
00179 /*                 Compute the current block row. */
00180 
00181                     i__3 = *n - j - jb + 1;
00182                     i__4 = j - 1;
00183                     template_blas_gemm("Transpose", "No transpose", &jb, &i__3, &i__4, &
00184                             c_b13, &a_ref(1, j), lda, &a_ref(1, j + jb), lda, 
00185                             &c_b14, &a_ref(j, j + jb), lda);
00186                     i__3 = *n - j - jb + 1;
00187                     template_blas_trsm("Left", "Upper", "Transpose", "Non-unit", &jb, &
00188                             i__3, &c_b14, &a_ref(j, j), lda, &a_ref(j, j + jb)
00189                             , lda)
00190                             ;
00191                 }
00192 /* L10: */
00193             }
00194 
00195         } else {
00196 
00197 /*           Compute the Cholesky factorization A = L*L'. */
00198 
00199             i__2 = *n;
00200             i__1 = nb;
00201             for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
00202 
00203 /*              Update and factorize the current diagonal block and test   
00204                 for non-positive-definiteness.   
00205 
00206    Computing MIN */
00207                 i__3 = nb, i__4 = *n - j + 1;
00208                 jb = minMACRO(i__3,i__4);
00209                 i__3 = j - 1;
00210                 template_blas_syrk("Lower", "No transpose", &jb, &i__3, &c_b13, &a_ref(j, 
00211                         1), lda, &c_b14, &a_ref(j, j), lda);
00212                 template_lapack_potf2("Lower", &jb, &a_ref(j, j), lda, info);
00213                 if (*info != 0) {
00214                     goto L30;
00215                 }
00216                 if (j + jb <= *n) {
00217 
00218 /*                 Compute the current block column. */
00219 
00220                     i__3 = *n - j - jb + 1;
00221                     i__4 = j - 1;
00222                     template_blas_gemm("No transpose", "Transpose", &i__3, &jb, &i__4, &
00223                             c_b13, &a_ref(j + jb, 1), lda, &a_ref(j, 1), lda, 
00224                             &c_b14, &a_ref(j + jb, j), lda);
00225                     i__3 = *n - j - jb + 1;
00226                     template_blas_trsm("Right", "Lower", "Transpose", "Non-unit", &i__3, &
00227                             jb, &c_b14, &a_ref(j, j), lda, &a_ref(j + jb, j), 
00228                             lda);
00229                 }
00230 /* L20: */
00231             }
00232         }
00233     }
00234     goto L40;
00235 
00236 L30:
00237     *info = *info + j - 1;
00238 
00239 L40:
00240     return 0;
00241 
00242 /*     End of DPOTRF */
00243 
00244 } /* dpotrf_ */
00245 
00246 #undef a_ref
00247 
00248 
00249 #endif

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