ergo
template_lapack_potrf.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
30 /* This file belongs to the template_lapack part of the Ergo source
31 * code. The source files in the template_lapack directory are modified
32 * versions of files originally distributed as CLAPACK, see the
33 * Copyright/license notice in the file template_lapack/COPYING.
34 */
35
36
37#ifndef TEMPLATE_LAPACK_POTRF_HEADER
38#define TEMPLATE_LAPACK_POTRF_HEADER
39
41
42template<class Treal>
43int template_lapack_potrf(const char *uplo, const integer *n, Treal *a, const integer *
44 lda, integer *info)
45{
46/* -- LAPACK routine (version 3.0) --
47 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
48 Courant Institute, Argonne National Lab, and Rice University
49 March 31, 1993
50
51
52 Purpose
53 =======
54
55 DPOTRF computes the Cholesky factorization of a real symmetric
56 positive definite matrix A.
57
58 The factorization has the form
59 A = U**T * U, if UPLO = 'U', or
60 A = L * L**T, if UPLO = 'L',
61 where U is an upper triangular matrix and L is lower triangular.
62
63 This is the block version of the algorithm, calling Level 3 BLAS.
64
65 Arguments
66 =========
67
68 UPLO (input) CHARACTER*1
69 = 'U': Upper triangle of A is stored;
70 = 'L': Lower triangle of A is stored.
71
72 N (input) INTEGER
73 The order of the matrix A. N >= 0.
74
75 A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
76 On entry, the symmetric matrix A. If UPLO = 'U', the leading
77 N-by-N upper triangular part of A contains the upper
78 triangular part of the matrix A, and the strictly lower
79 triangular part of A is not referenced. If UPLO = 'L', the
80 leading N-by-N lower triangular part of A contains the lower
81 triangular part of the matrix A, and the strictly upper
82 triangular part of A is not referenced.
83
84 On exit, if INFO = 0, the factor U or L from the Cholesky
85 factorization A = U**T*U or A = L*L**T.
86
87 LDA (input) INTEGER
88 The leading dimension of the array A. LDA >= max(1,N).
89
90 INFO (output) INTEGER
91 = 0: successful exit
92 < 0: if INFO = -i, the i-th argument had an illegal value
93 > 0: if INFO = i, the leading minor of order i is not
94 positive definite, and the factorization could not be
95 completed.
96
97 =====================================================================
98
99
100 Test the input parameters.
101
102 Parameter adjustments */
103 /* Table of constant values */
104 integer c__1 = 1;
105 integer c_n1 = -1;
106 Treal c_b13 = -1.;
107 Treal c_b14 = 1.;
108
109 /* System generated locals */
110 integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
111 /* Local variables */
112 integer j;
113 logical upper;
114 integer jb, nb;
115#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
116
117
118 a_dim1 = *lda;
119 a_offset = 1 + a_dim1 * 1;
120 a -= a_offset;
121
122 /* Function Body */
123 *info = 0;
124 upper = template_blas_lsame(uplo, "U");
125 if (! upper && ! template_blas_lsame(uplo, "L")) {
126 *info = -1;
127 } else if (*n < 0) {
128 *info = -2;
129 } else if (*lda < maxMACRO(1,*n)) {
130 *info = -4;
131 }
132 if (*info != 0) {
133 i__1 = -(*info);
134 template_blas_erbla("POTRF ", &i__1);
135 return 0;
136 }
137
138/* Quick return if possible */
139
140 if (*n == 0) {
141 return 0;
142 }
143
144/* Determine the block size for this environment. */
145
146 nb = template_lapack_ilaenv(&c__1, "DPOTRF", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (
147 ftnlen)1);
148 if (nb <= 1 || nb >= *n) {
149
150/* Use unblocked code. */
151
152 template_lapack_potf2(uplo, n, &a[a_offset], lda, info);
153 } else {
154
155/* Use blocked code. */
156
157 if (upper) {
158
159/* Compute the Cholesky factorization A = U'*U. */
160
161 i__1 = *n;
162 i__2 = nb;
163 for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
164
165/* Update and factorize the current diagonal block and test
166 for non-positive-definiteness.
167
168 Computing MIN */
169 i__3 = nb, i__4 = *n - j + 1;
170 jb = minMACRO(i__3,i__4);
171 i__3 = j - 1;
172 template_blas_syrk("Upper", "Transpose", &jb, &i__3, &c_b13, &a_ref(1, j),
173 lda, &c_b14, &a_ref(j, j), lda)
174 ;
175 template_lapack_potf2("Upper", &jb, &a_ref(j, j), lda, info);
176 if (*info != 0) {
177 goto L30;
178 }
179 if (j + jb <= *n) {
180
181/* Compute the current block row. */
182
183 i__3 = *n - j - jb + 1;
184 i__4 = j - 1;
185 template_blas_gemm("Transpose", "No transpose", &jb, &i__3, &i__4, &
186 c_b13, &a_ref(1, j), lda, &a_ref(1, j + jb), lda,
187 &c_b14, &a_ref(j, j + jb), lda);
188 i__3 = *n - j - jb + 1;
189 template_blas_trsm("Left", "Upper", "Transpose", "Non-unit", &jb, &
190 i__3, &c_b14, &a_ref(j, j), lda, &a_ref(j, j + jb)
191 , lda)
192 ;
193 }
194/* L10: */
195 }
196
197 } else {
198
199/* Compute the Cholesky factorization A = L*L'. */
200
201 i__2 = *n;
202 i__1 = nb;
203 for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
204
205/* Update and factorize the current diagonal block and test
206 for non-positive-definiteness.
207
208 Computing MIN */
209 i__3 = nb, i__4 = *n - j + 1;
210 jb = minMACRO(i__3,i__4);
211 i__3 = j - 1;
212 template_blas_syrk("Lower", "No transpose", &jb, &i__3, &c_b13, &a_ref(j,
213 1), lda, &c_b14, &a_ref(j, j), lda);
214 template_lapack_potf2("Lower", &jb, &a_ref(j, j), lda, info);
215 if (*info != 0) {
216 goto L30;
217 }
218 if (j + jb <= *n) {
219
220/* Compute the current block column. */
221
222 i__3 = *n - j - jb + 1;
223 i__4 = j - 1;
224 template_blas_gemm("No transpose", "Transpose", &i__3, &jb, &i__4, &
225 c_b13, &a_ref(j + jb, 1), lda, &a_ref(j, 1), lda,
226 &c_b14, &a_ref(j + jb, j), lda);
227 i__3 = *n - j - jb + 1;
228 template_blas_trsm("Right", "Lower", "Transpose", "Non-unit", &i__3, &
229 jb, &c_b14, &a_ref(j, j), lda, &a_ref(j + jb, j),
230 lda);
231 }
232/* L20: */
233 }
234 }
235 }
236 goto L40;
237
238L30:
239 *info = *info + j - 1;
240
241L40:
242 return 0;
243
244/* End of DPOTRF */
245
246} /* dpotrf_ */
247
248#undef a_ref
249
250
251#endif
int template_blas_erbla(const char *srname, integer *info)
Definition template_blas_common.cc:146
logical template_blas_lsame(const char *ca, const char *cb)
Definition template_blas_common.cc:46
int integer
Definition template_blas_common.h:40
int ftnlen
Definition template_blas_common.h:42
#define minMACRO(a, b)
Definition template_blas_common.h:46
#define maxMACRO(a, b)
Definition template_blas_common.h:45
bool logical
Definition template_blas_common.h:41
int template_blas_gemm(const char *transa, const char *transb, const integer *m, const integer *n, const integer *k, const Treal *alpha, const Treal *a, const integer *lda, const Treal *b, const integer *ldb, const Treal *beta, Treal *c__, const integer *ldc)
Definition template_blas_gemm.h:43
int template_blas_syrk(const char *uplo, const char *trans, const integer *n, const integer *k, const Treal *alpha, const Treal *a, const integer *lda, const Treal *beta, Treal *c__, const integer *ldc)
Definition template_blas_syrk.h:42
int template_blas_trsm(const char *side, const char *uplo, const char *transa, const char *diag, const integer *m, const integer *n, const Treal *alpha, const Treal *a, const integer *lda, Treal *b, const integer *ldb)
Definition template_blas_trsm.h:43
integer template_lapack_ilaenv(const integer *ispec, const char *name__, const char *opts, const integer *n1, const integer *n2, const integer *n3, const integer *n4, ftnlen name_len, ftnlen opts_len)
Definition template_lapack_common.cc:281
int template_lapack_potf2(const char *uplo, const integer *n, Treal *a, const integer *lda, integer *info)
Definition template_lapack_potf2.h:42
#define a_ref(a_1, a_2)
int template_lapack_potrf(const char *uplo, const integer *n, Treal *a, const integer *lda, integer *info)
Definition template_lapack_potrf.h:43