ergo
Main Page
Namespaces
Classes
Files
File List
File Members
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
source
matrix
template_lapack
lapack
template_lapack_pptrf.h
Generated on Wed Nov 21 2012 09:32:25 for ergo by
1.8.1.2