ergo
template_lapack_potrf.h
Go to the documentation of this file.
1 /* Ergo, version 3.4, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2014 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_POTRF_HEADER
36 #define TEMPLATE_LAPACK_POTRF_HEADER
37 
38 #include "template_lapack_potf2.h"
39 
40 template<class Treal>
41 int template_lapack_potrf(const char *uplo, const integer *n, Treal *a, const integer *
42  lda, integer *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  DPOTRF computes the Cholesky factorization of a real symmetric
54  positive definite matrix A.
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  This is the block version of the algorithm, calling Level 3 BLAS.
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  A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
74  On entry, the symmetric matrix A. If UPLO = 'U', the leading
75  N-by-N upper triangular part of A contains the upper
76  triangular part of the matrix A, and the strictly lower
77  triangular part of A is not referenced. If UPLO = 'L', the
78  leading N-by-N lower triangular part of A contains the lower
79  triangular part of the matrix A, and the strictly upper
80  triangular part of A is not referenced.
81 
82  On exit, if INFO = 0, the factor U or L from the Cholesky
83  factorization A = U**T*U or A = L*L**T.
84 
85  LDA (input) INTEGER
86  The leading dimension of the array A. LDA >= max(1,N).
87 
88  INFO (output) INTEGER
89  = 0: successful exit
90  < 0: if INFO = -i, the i-th argument had an illegal value
91  > 0: if INFO = i, the leading minor of order i is not
92  positive definite, and the factorization could not be
93  completed.
94 
95  =====================================================================
96 
97 
98  Test the input parameters.
99 
100  Parameter adjustments */
101  /* Table of constant values */
102  integer c__1 = 1;
103  integer c_n1 = -1;
104  Treal c_b13 = -1.;
105  Treal c_b14 = 1.;
106 
107  /* System generated locals */
108  integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
109  /* Local variables */
110  integer j;
111  logical upper;
112  integer jb, nb;
113 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
114 
115 
116  a_dim1 = *lda;
117  a_offset = 1 + a_dim1 * 1;
118  a -= a_offset;
119 
120  /* Function Body */
121  *info = 0;
122  upper = template_blas_lsame(uplo, "U");
123  if (! upper && ! template_blas_lsame(uplo, "L")) {
124  *info = -1;
125  } else if (*n < 0) {
126  *info = -2;
127  } else if (*lda < maxMACRO(1,*n)) {
128  *info = -4;
129  }
130  if (*info != 0) {
131  i__1 = -(*info);
132  template_blas_erbla("POTRF ", &i__1);
133  return 0;
134  }
135 
136 /* Quick return if possible */
137 
138  if (*n == 0) {
139  return 0;
140  }
141 
142 /* Determine the block size for this environment. */
143 
144  nb = template_lapack_ilaenv(&c__1, "DPOTRF", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (
145  ftnlen)1);
146  if (nb <= 1 || nb >= *n) {
147 
148 /* Use unblocked code. */
149 
150  template_lapack_potf2(uplo, n, &a[a_offset], lda, info);
151  } else {
152 
153 /* Use blocked code. */
154 
155  if (upper) {
156 
157 /* Compute the Cholesky factorization A = U'*U. */
158 
159  i__1 = *n;
160  i__2 = nb;
161  for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
162 
163 /* Update and factorize the current diagonal block and test
164  for non-positive-definiteness.
165 
166  Computing MIN */
167  i__3 = nb, i__4 = *n - j + 1;
168  jb = minMACRO(i__3,i__4);
169  i__3 = j - 1;
170  template_blas_syrk("Upper", "Transpose", &jb, &i__3, &c_b13, &a_ref(1, j),
171  lda, &c_b14, &a_ref(j, j), lda)
172  ;
173  template_lapack_potf2("Upper", &jb, &a_ref(j, j), lda, info);
174  if (*info != 0) {
175  goto L30;
176  }
177  if (j + jb <= *n) {
178 
179 /* Compute the current block row. */
180 
181  i__3 = *n - j - jb + 1;
182  i__4 = j - 1;
183  template_blas_gemm("Transpose", "No transpose", &jb, &i__3, &i__4, &
184  c_b13, &a_ref(1, j), lda, &a_ref(1, j + jb), lda,
185  &c_b14, &a_ref(j, j + jb), lda);
186  i__3 = *n - j - jb + 1;
187  template_blas_trsm("Left", "Upper", "Transpose", "Non-unit", &jb, &
188  i__3, &c_b14, &a_ref(j, j), lda, &a_ref(j, j + jb)
189  , lda)
190  ;
191  }
192 /* L10: */
193  }
194 
195  } else {
196 
197 /* Compute the Cholesky factorization A = L*L'. */
198 
199  i__2 = *n;
200  i__1 = nb;
201  for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
202 
203 /* Update and factorize the current diagonal block and test
204  for non-positive-definiteness.
205 
206  Computing MIN */
207  i__3 = nb, i__4 = *n - j + 1;
208  jb = minMACRO(i__3,i__4);
209  i__3 = j - 1;
210  template_blas_syrk("Lower", "No transpose", &jb, &i__3, &c_b13, &a_ref(j,
211  1), lda, &c_b14, &a_ref(j, j), lda);
212  template_lapack_potf2("Lower", &jb, &a_ref(j, j), lda, info);
213  if (*info != 0) {
214  goto L30;
215  }
216  if (j + jb <= *n) {
217 
218 /* Compute the current block column. */
219 
220  i__3 = *n - j - jb + 1;
221  i__4 = j - 1;
222  template_blas_gemm("No transpose", "Transpose", &i__3, &jb, &i__4, &
223  c_b13, &a_ref(j + jb, 1), lda, &a_ref(j, 1), lda,
224  &c_b14, &a_ref(j + jb, j), lda);
225  i__3 = *n - j - jb + 1;
226  template_blas_trsm("Right", "Lower", "Transpose", "Non-unit", &i__3, &
227  jb, &c_b14, &a_ref(j, j), lda, &a_ref(j + jb, j),
228  lda);
229  }
230 /* L20: */
231  }
232  }
233  }
234  goto L40;
235 
236 L30:
237  *info = *info + j - 1;
238 
239 L40:
240  return 0;
241 
242 /* End of DPOTRF */
243 
244 } /* dpotrf_ */
245 
246 #undef a_ref
247 
248 
249 #endif
int template_lapack_potrf(const char *uplo, const integer *n, Treal *a, const integer *lda, integer *info)
Definition: template_lapack_potrf.h:41
int integer
Definition: template_blas_common.h:38
#define a_ref(a_1, a_2)
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:279
#define maxMACRO(a, b)
Definition: template_blas_common.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:40
#define minMACRO(a, b)
Definition: template_blas_common.h:44
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
int template_lapack_potf2(const char *uplo, const integer *n, Treal *a, const integer *lda, integer *info)
Definition: template_lapack_potf2.h:40
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:41
bool logical
Definition: template_blas_common.h:39
int ftnlen
Definition: template_blas_common.h:40
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:41
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:44