ergo
template_lapack_pptrf.h
Go to the documentation of this file.
1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2019 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 
40 #include "template_lapack_common.h"
41 
42 template<class Treal>
43 int 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 
215 L30:
216  *info = j;
217 
218 L40:
219  return 0;
220 
221 /* End of DPPTRF */
222 
223 } /* dpptrf_ */
224 
225 #endif
template_blas_sqrt
Treal template_blas_sqrt(Treal x)
template_blas_spr
int template_blas_spr(const char *uplo, integer *n, Treal *alpha, Treal *x, integer *incx, Treal *ap)
Definition: template_blas_spr.h:43
template_blas_dot
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
template_blas_scal
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition: template_blas_scal.h:43
logical
bool logical
Definition: template_blas_common.h:41
template_lapack_common.h
template_blas_tpsv
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
template_blas_erbla
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
template_blas_lsame
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
integer
int integer
Definition: template_blas_common.h:40
template_lapack_pptrf
int template_lapack_pptrf(const char *uplo, const integer *n, Treal *ap, integer *info)
Definition: template_lapack_pptrf.h:43