ergo
template_lapack_potf2.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_POTF2_HEADER
36 #define TEMPLATE_LAPACK_POTF2_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_potf2(const char *uplo, const integer *n, Treal *a, const integer *
41  lda, integer *info)
42 {
43 /* -- LAPACK routine (version 3.0) --
44  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
45  Courant Institute, Argonne National Lab, and Rice University
46  February 29, 1992
47 
48 
49  Purpose
50  =======
51 
52  DPOTF2 computes the Cholesky factorization of a real symmetric
53  positive definite matrix A.
54 
55  The factorization has the form
56  A = U' * U , if UPLO = 'U', or
57  A = L * L', if UPLO = 'L',
58  where U is an upper triangular matrix and L is lower triangular.
59 
60  This is the unblocked version of the algorithm, calling Level 2 BLAS.
61 
62  Arguments
63  =========
64 
65  UPLO (input) CHARACTER*1
66  Specifies whether the upper or lower triangular part of the
67  symmetric matrix A is stored.
68  = 'U': Upper triangular
69  = 'L': Lower triangular
70 
71  N (input) INTEGER
72  The order of the matrix A. N >= 0.
73 
74  A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
75  On entry, the symmetric matrix A. If UPLO = 'U', the leading
76  n by n upper triangular part of A contains the upper
77  triangular part of the matrix A, and the strictly lower
78  triangular part of A is not referenced. If UPLO = 'L', the
79  leading n by n lower triangular part of A contains the lower
80  triangular part of the matrix A, and the strictly upper
81  triangular part of A is not referenced.
82 
83  On exit, if INFO = 0, the factor U or L from the Cholesky
84  factorization A = U'*U or A = L*L'.
85 
86  LDA (input) INTEGER
87  The leading dimension of the array A. LDA >= max(1,N).
88 
89  INFO (output) INTEGER
90  = 0: successful exit
91  < 0: if INFO = -k, the k-th argument had an illegal value
92  > 0: if INFO = k, the leading minor of order k is not
93  positive definite, and the factorization could not be
94  completed.
95 
96  =====================================================================
97 
98 
99  Test the input parameters.
100 
101  Parameter adjustments */
102  /* Table of constant values */
103  integer c__1 = 1;
104  Treal c_b10 = -1.;
105  Treal c_b12 = 1.;
106 
107  /* System generated locals */
108  integer a_dim1, a_offset, i__1, i__2, i__3;
109  Treal d__1;
110  /* Local variables */
111  integer j;
112  logical upper;
113  Treal ajj;
114 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
115 
116 
117  a_dim1 = *lda;
118  a_offset = 1 + a_dim1 * 1;
119  a -= a_offset;
120 
121  /* Function Body */
122  *info = 0;
123  upper = template_blas_lsame(uplo, "U");
124  if (! upper && ! template_blas_lsame(uplo, "L")) {
125  *info = -1;
126  } else if (*n < 0) {
127  *info = -2;
128  } else if (*lda < maxMACRO(1,*n)) {
129  *info = -4;
130  }
131  if (*info != 0) {
132  i__1 = -(*info);
133  template_blas_erbla("POTF2 ", &i__1);
134  return 0;
135  }
136 
137 /* Quick return if possible */
138 
139  if (*n == 0) {
140  return 0;
141  }
142 
143  if (upper) {
144 
145 /* Compute the Cholesky factorization A = U'*U. */
146 
147  i__1 = *n;
148  for (j = 1; j <= i__1; ++j) {
149 
150 /* Compute U(J,J) and test for non-positive-definiteness. */
151 
152  i__2 = j - 1;
153  ajj = a_ref(j, j) - template_blas_dot(&i__2, &a_ref(1, j), &c__1, &a_ref(1, j)
154  , &c__1);
155  if (ajj <= 0.) {
156  a_ref(j, j) = ajj;
157  goto L30;
158  }
159  ajj = template_blas_sqrt(ajj);
160  a_ref(j, j) = ajj;
161 
162 /* Compute elements J+1:N of row J. */
163 
164  if (j < *n) {
165  i__2 = j - 1;
166  i__3 = *n - j;
167  template_blas_gemv("Transpose", &i__2, &i__3, &c_b10, &a_ref(1, j + 1),
168  lda, &a_ref(1, j), &c__1, &c_b12, &a_ref(j, j + 1),
169  lda);
170  i__2 = *n - j;
171  d__1 = 1. / ajj;
172  template_blas_scal(&i__2, &d__1, &a_ref(j, j + 1), lda);
173  }
174 /* L10: */
175  }
176  } else {
177 
178 /* Compute the Cholesky factorization A = L*L'. */
179 
180  i__1 = *n;
181  for (j = 1; j <= i__1; ++j) {
182 
183 /* Compute L(J,J) and test for non-positive-definiteness. */
184 
185  i__2 = j - 1;
186  ajj = a_ref(j, j) - template_blas_dot(&i__2, &a_ref(j, 1), lda, &a_ref(j, 1),
187  lda);
188  if (ajj <= 0.) {
189  a_ref(j, j) = ajj;
190  goto L30;
191  }
192  ajj = template_blas_sqrt(ajj);
193  a_ref(j, j) = ajj;
194 
195 /* Compute elements J+1:N of column J. */
196 
197  if (j < *n) {
198  i__2 = *n - j;
199  i__3 = j - 1;
200  template_blas_gemv("No transpose", &i__2, &i__3, &c_b10, &a_ref(j + 1, 1),
201  lda, &a_ref(j, 1), lda, &c_b12, &a_ref(j + 1, j), &
202  c__1);
203  i__2 = *n - j;
204  d__1 = 1. / ajj;
205  template_blas_scal(&i__2, &d__1, &a_ref(j + 1, j), &c__1);
206  }
207 /* L20: */
208  }
209  }
210  goto L40;
211 
212 L30:
213  *info = j;
214 
215 L40:
216  return 0;
217 
218 /* End of DPOTF2 */
219 
220 } /* dpotf2_ */
221 
222 #undef a_ref
223 
224 
225 #endif