ergo
|
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