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