ergo
template_lapack_tptri.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 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_TPTRI_HEADER
38#define TEMPLATE_LAPACK_TPTRI_HEADER
39
41
42template<class Treal>
43int template_lapack_tptri(const char *uplo, const char *diag, const integer *n, Treal *
44 ap, 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 September 30, 1994
50
51
52 Purpose
53 =======
54
55 DTPTRI computes the inverse of a real upper or lower triangular
56 matrix A stored in packed format.
57
58 Arguments
59 =========
60
61 UPLO (input) CHARACTER*1
62 = 'U': A is upper triangular;
63 = 'L': A is lower triangular.
64
65 DIAG (input) CHARACTER*1
66 = 'N': A is non-unit triangular;
67 = 'U': A is unit triangular.
68
69 N (input) INTEGER
70 The order of the matrix A. N >= 0.
71
72 AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
73 On entry, the upper or lower triangular matrix A, stored
74 columnwise in a linear array. The j-th column of A is stored
75 in the array AP as follows:
76 if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
77 if UPLO = 'L', AP(i + (j-1)*((2*n-j)/2) = A(i,j) for j<=i<=n.
78 See below for further details.
79 On exit, the (triangular) inverse of the original matrix, in
80 the same packed storage format.
81
82 INFO (output) INTEGER
83 = 0: successful exit
84 < 0: if INFO = -i, the i-th argument had an illegal value
85 > 0: if INFO = i, A(i,i) is exactly zero. The triangular
86 matrix is singular and its inverse can not be computed.
87
88 Further Details
89 ===============
90
91 A triangular matrix A can be transferred to packed storage using one
92 of the following program segments:
93
94 UPLO = 'U': UPLO = 'L':
95
96 JC = 1 JC = 1
97 DO 2 J = 1, N DO 2 J = 1, N
98 DO 1 I = 1, J DO 1 I = J, N
99 AP(JC+I-1) = A(I,J) AP(JC+I-J) = A(I,J)
100 1 CONTINUE 1 CONTINUE
101 JC = JC + J JC = JC + N - J + 1
102 2 CONTINUE 2 CONTINUE
103
104 =====================================================================
105
106
107 Test the input parameters.
108
109 Parameter adjustments */
110 /* Table of constant values */
111 integer c__1 = 1;
112
113 /* System generated locals */
114 integer i__1, i__2;
115 /* Local variables */
116 integer j;
117 logical upper;
118 integer jc, jj;
119 integer jclast;
120 logical nounit;
121 Treal ajj;
122
123
124 --ap;
125
126 /* Initialization added by Elias to get rid of compiler warnings. */
127 jclast = 0;
128 /* Function Body */
129 *info = 0;
130 upper = template_blas_lsame(uplo, "U");
131 nounit = template_blas_lsame(diag, "N");
132 if (! upper && ! template_blas_lsame(uplo, "L")) {
133 *info = -1;
134 } else if (! nounit && ! template_blas_lsame(diag, "U")) {
135 *info = -2;
136 } else if (*n < 0) {
137 *info = -3;
138 }
139 if (*info != 0) {
140 i__1 = -(*info);
141 template_blas_erbla("TPTRI ", &i__1);
142 return 0;
143 }
144
145/* Check for singularity if non-unit. */
146
147 if (nounit) {
148 if (upper) {
149 jj = 0;
150 i__1 = *n;
151 for (*info = 1; *info <= i__1; ++(*info)) {
152 jj += *info;
153 if (ap[jj] == 0.) {
154 return 0;
155 }
156/* L10: */
157 }
158 } else {
159 jj = 1;
160 i__1 = *n;
161 for (*info = 1; *info <= i__1; ++(*info)) {
162 if (ap[jj] == 0.) {
163 return 0;
164 }
165 jj = jj + *n - *info + 1;
166/* L20: */
167 }
168 }
169 *info = 0;
170 }
171
172 if (upper) {
173
174/* Compute inverse of upper triangular matrix. */
175
176 jc = 1;
177 i__1 = *n;
178 for (j = 1; j <= i__1; ++j) {
179 if (nounit) {
180 ap[jc + j - 1] = 1. / ap[jc + j - 1];
181 ajj = -ap[jc + j - 1];
182 } else {
183 ajj = -1.;
184 }
185
186/* Compute elements 1:j-1 of j-th column. */
187
188 i__2 = j - 1;
189 template_blas_tpmv("Upper", "No transpose", diag, &i__2, &ap[1], &ap[jc], &
190 c__1);
191 i__2 = j - 1;
192 template_blas_scal(&i__2, &ajj, &ap[jc], &c__1);
193 jc += j;
194/* L30: */
195 }
196
197 } else {
198
199/* Compute inverse of lower triangular matrix. */
200
201 jc = *n * (*n + 1) / 2;
202 for (j = *n; j >= 1; --j) {
203 if (nounit) {
204 ap[jc] = 1. / ap[jc];
205 ajj = -ap[jc];
206 } else {
207 ajj = -1.;
208 }
209 if (j < *n) {
210
211/* Compute elements j+1:n of j-th column. */
212
213 i__1 = *n - j;
214 template_blas_tpmv("Lower", "No transpose", diag, &i__1, &ap[jclast], &ap[
215 jc + 1], &c__1);
216 i__1 = *n - j;
217 template_blas_scal(&i__1, &ajj, &ap[jc + 1], &c__1);
218 }
219 jclast = jc;
220 jc = jc - *n + j - 2;
221/* L40: */
222 }
223 }
224
225 return 0;
226
227/* End of DTPTRI */
228
229} /* dtptri_ */
230
231#endif
int template_blas_erbla(const char *srname, integer *info)
Definition template_blas_common.cc:146
logical template_blas_lsame(const char *ca, const char *cb)
Definition template_blas_common.cc:46
int integer
Definition template_blas_common.h:40
bool logical
Definition template_blas_common.h:41
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition template_blas_scal.h:43
int template_blas_tpmv(const char *uplo, const char *trans, const char *diag, const integer *n, const Treal *ap, Treal *x, const integer *incx)
Definition template_blas_tpmv.h:42
int template_lapack_tptri(const char *uplo, const char *diag, const integer *n, Treal *ap, integer *info)
Definition template_lapack_tptri.h:43