ergo
template_lapack_pptrf.h
Go to the documentation of this file.
1 /* Ergo, version 3.2, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28  /* This file belongs to the template_lapack part of the Ergo source
29  * code. The source files in the template_lapack directory are modified
30  * versions of files originally distributed as CLAPACK, see the
31  * Copyright/license notice in the file template_lapack/COPYING.
32  */
33 
34 
35 #ifndef TEMPLATE_LAPACK_PPTRF_HEADER
36 #define TEMPLATE_LAPACK_PPTRF_HEADER
37 
38 #include "template_lapack_common.h"
39 
40 template<class Treal>
41 int template_lapack_pptrf(const char *uplo, const integer *n, Treal *ap, integer *
42  info)
43 {
44 /* -- LAPACK routine (version 3.0) --
45  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
46  Courant Institute, Argonne National Lab, and Rice University
47  March 31, 1993
48 
49 
50  Purpose
51  =======
52 
53  DPPTRF computes the Cholesky factorization of a real symmetric
54  positive definite matrix A stored in packed format.
55 
56  The factorization has the form
57  A = U**T * U, if UPLO = 'U', or
58  A = L * L**T, if UPLO = 'L',
59  where U is an upper triangular matrix and L is lower triangular.
60 
61  Arguments
62  =========
63 
64  UPLO (input) CHARACTER*1
65  = 'U': Upper triangle of A is stored;
66  = 'L': Lower triangle of A is stored.
67 
68  N (input) INTEGER
69  The order of the matrix A. N >= 0.
70 
71  AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
72  On entry, the upper or lower triangle of the symmetric matrix
73  A, packed columnwise in a linear array. The j-th column of A
74  is stored in the array AP as follows:
75  if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
76  if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
77  See below for further details.
78 
79  On exit, if INFO = 0, the triangular factor U or L from the
80  Cholesky factorization A = U**T*U or A = L*L**T, in the same
81  storage format as A.
82 
83  INFO (output) INTEGER
84  = 0: successful exit
85  < 0: if INFO = -i, the i-th argument had an illegal value
86  > 0: if INFO = i, the leading minor of order i is not
87  positive definite, and the factorization could not be
88  completed.
89 
90  Further Details
91  ======= =======
92 
93  The packed storage scheme is illustrated by the following example
94  when N = 4, UPLO = 'U':
95 
96  Two-dimensional storage of the symmetric matrix A:
97 
98  a11 a12 a13 a14
99  a22 a23 a24
100  a33 a34 (aij = aji)
101  a44
102 
103  Packed storage of the upper triangle of A:
104 
105  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
106 
107  =====================================================================
108 
109 
110  Test the input parameters.
111 
112  Parameter adjustments */
113  /* Table of constant values */
114  integer c__1 = 1;
115  Treal c_b16 = -1.;
116 
117  /* System generated locals */
118  integer i__1, i__2;
119  Treal d__1;
120  /* Local variables */
121  integer j;
122  logical upper;
123  integer jc, jj;
124  Treal ajj;
125 
126 
127  --ap;
128 
129  /* Function Body */
130  *info = 0;
131  upper = template_blas_lsame(uplo, "U");
132  if (! upper && ! template_blas_lsame(uplo, "L")) {
133  *info = -1;
134  } else if (*n < 0) {
135  *info = -2;
136  }
137  if (*info != 0) {
138  i__1 = -(*info);
139  template_blas_erbla("DPPTRF", &i__1);
140  return 0;
141  }
142 
143 /* Quick return if possible */
144 
145  if (*n == 0) {
146  return 0;
147  }
148 
149  if (upper) {
150 
151 /* Compute the Cholesky factorization A = U'*U. */
152 
153  jj = 0;
154  i__1 = *n;
155  for (j = 1; j <= i__1; ++j) {
156  jc = jj + 1;
157  jj += j;
158 
159 /* Compute elements 1:J-1 of column J. */
160 
161  if (j > 1) {
162  i__2 = j - 1;
163  template_blas_tpsv("Upper", "Transpose", "Non-unit", &i__2, &ap[1], &ap[
164  jc], &c__1);
165  }
166 
167 /* Compute U(J,J) and test for non-positive-definiteness. */
168 
169  i__2 = j - 1;
170  ajj = ap[jj] - template_blas_dot(&i__2, &ap[jc], &c__1, &ap[jc], &c__1);
171  if (ajj <= 0.) {
172  ap[jj] = ajj;
173  goto L30;
174  }
175  ap[jj] = template_blas_sqrt(ajj);
176 /* L10: */
177  }
178  } else {
179 
180 /* Compute the Cholesky factorization A = L*L'. */
181 
182  jj = 1;
183  i__1 = *n;
184  for (j = 1; j <= i__1; ++j) {
185 
186 /* Compute L(J,J) and test for non-positive-definiteness. */
187 
188  ajj = ap[jj];
189  if (ajj <= 0.) {
190  ap[jj] = ajj;
191  goto L30;
192  }
193  ajj = template_blas_sqrt(ajj);
194  ap[jj] = ajj;
195 
196 /* Compute elements J+1:N of column J and update the trailing
197  submatrix. */
198 
199  if (j < *n) {
200  i__2 = *n - j;
201  d__1 = 1. / ajj;
202  template_blas_scal(&i__2, &d__1, &ap[jj + 1], &c__1);
203  i__2 = *n - j;
204  template_blas_spr("Lower", &i__2, &c_b16, &ap[jj + 1], &c__1, &ap[jj + *n
205  - j + 1]);
206  jj = jj + *n - j + 1;
207  }
208 /* L20: */
209  }
210  }
211  goto L40;
212 
213 L30:
214  *info = j;
215 
216 L40:
217  return 0;
218 
219 /* End of DPPTRF */
220 
221 } /* dpptrf_ */
222 
223 #endif