ergo
template_lapack_pptrf.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_PPTRF_HEADER
38#define TEMPLATE_LAPACK_PPTRF_HEADER
39
41
42template<class Treal>
43int template_lapack_pptrf(const char *uplo, const integer *n, Treal *ap, integer *
44 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 DPPTRF computes the Cholesky factorization of a real symmetric
56 positive definite matrix A stored in packed format.
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 Arguments
64 =========
65
66 UPLO (input) CHARACTER*1
67 = 'U': Upper triangle of A is stored;
68 = 'L': Lower triangle of A is stored.
69
70 N (input) INTEGER
71 The order of the matrix A. N >= 0.
72
73 AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
74 On entry, the upper or lower triangle of the symmetric matrix
75 A, packed columnwise in a linear array. The j-th column of A
76 is stored in the array AP as follows:
77 if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
78 if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
79 See below for further details.
80
81 On exit, if INFO = 0, the triangular factor U or L from the
82 Cholesky factorization A = U**T*U or A = L*L**T, in the same
83 storage format as A.
84
85 INFO (output) INTEGER
86 = 0: successful exit
87 < 0: if INFO = -i, the i-th argument had an illegal value
88 > 0: if INFO = i, the leading minor of order i is not
89 positive definite, and the factorization could not be
90 completed.
91
92 Further Details
93 ======= =======
94
95 The packed storage scheme is illustrated by the following example
96 when N = 4, UPLO = 'U':
97
98 Two-dimensional storage of the symmetric matrix A:
99
100 a11 a12 a13 a14
101 a22 a23 a24
102 a33 a34 (aij = aji)
103 a44
104
105 Packed storage of the upper triangle of A:
106
107 AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
108
109 =====================================================================
110
111
112 Test the input parameters.
113
114 Parameter adjustments */
115 /* Table of constant values */
116 integer c__1 = 1;
117 Treal c_b16 = -1.;
118
119 /* System generated locals */
120 integer i__1, i__2;
121 Treal d__1;
122 /* Local variables */
123 integer j;
124 logical upper;
125 integer jc, jj;
126 Treal ajj;
127
128
129 --ap;
130
131 /* Function Body */
132 *info = 0;
133 upper = template_blas_lsame(uplo, "U");
134 if (! upper && ! template_blas_lsame(uplo, "L")) {
135 *info = -1;
136 } else if (*n < 0) {
137 *info = -2;
138 }
139 if (*info != 0) {
140 i__1 = -(*info);
141 template_blas_erbla("DPPTRF", &i__1);
142 return 0;
143 }
144
145/* Quick return if possible */
146
147 if (*n == 0) {
148 return 0;
149 }
150
151 if (upper) {
152
153/* Compute the Cholesky factorization A = U'*U. */
154
155 jj = 0;
156 i__1 = *n;
157 for (j = 1; j <= i__1; ++j) {
158 jc = jj + 1;
159 jj += j;
160
161/* Compute elements 1:J-1 of column J. */
162
163 if (j > 1) {
164 i__2 = j - 1;
165 template_blas_tpsv("Upper", "Transpose", "Non-unit", &i__2, &ap[1], &ap[
166 jc], &c__1);
167 }
168
169/* Compute U(J,J) and test for non-positive-definiteness. */
170
171 i__2 = j - 1;
172 ajj = ap[jj] - template_blas_dot(&i__2, &ap[jc], &c__1, &ap[jc], &c__1);
173 if (ajj <= 0.) {
174 ap[jj] = ajj;
175 goto L30;
176 }
177 ap[jj] = template_blas_sqrt(ajj);
178/* L10: */
179 }
180 } else {
181
182/* Compute the Cholesky factorization A = L*L'. */
183
184 jj = 1;
185 i__1 = *n;
186 for (j = 1; j <= i__1; ++j) {
187
188/* Compute L(J,J) and test for non-positive-definiteness. */
189
190 ajj = ap[jj];
191 if (ajj <= 0.) {
192 ap[jj] = ajj;
193 goto L30;
194 }
195 ajj = template_blas_sqrt(ajj);
196 ap[jj] = ajj;
197
198/* Compute elements J+1:N of column J and update the trailing
199 submatrix. */
200
201 if (j < *n) {
202 i__2 = *n - j;
203 d__1 = 1. / ajj;
204 template_blas_scal(&i__2, &d__1, &ap[jj + 1], &c__1);
205 i__2 = *n - j;
206 template_blas_spr("Lower", &i__2, &c_b16, &ap[jj + 1], &c__1, &ap[jj + *n
207 - j + 1]);
208 jj = jj + *n - j + 1;
209 }
210/* L20: */
211 }
212 }
213 goto L40;
214
215L30:
216 *info = j;
217
218L40:
219 return 0;
220
221/* End of DPPTRF */
222
223} /* dpptrf_ */
224
225#endif
Treal template_blas_sqrt(Treal x)
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
bool logical
Definition: template_blas_common.h:41
Treal template_blas_dot(const integer *n, const Treal *dx, const integer *incx, const Treal *dy, const integer *incy)
Definition: template_blas_dot.h:43
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition: template_blas_scal.h:43
int template_blas_spr(const char *uplo, integer *n, Treal *alpha, Treal *x, integer *incx, Treal *ap)
Definition: template_blas_spr.h:43
int template_blas_tpsv(const char *uplo, const char *trans, const char *diag, const integer *n, const Treal *ap, Treal *x, const integer *incx)
Definition: template_blas_tpsv.h:42
int template_lapack_pptrf(const char *uplo, const integer *n, Treal *ap, integer *info)
Definition: template_lapack_pptrf.h:43