ergo
template_lapack_trti2.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_TRTI2_HEADER
38#define TEMPLATE_LAPACK_TRTI2_HEADER
39
40
41template<class Treal>
42int template_lapack_trti2(const char *uplo, const char *diag, const integer *n, Treal *
43 a, const integer *lda, integer *info)
44{
45/* -- LAPACK routine (version 3.0) --
46 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
47 Courant Institute, Argonne National Lab, and Rice University
48 February 29, 1992
49
50
51 Purpose
52 =======
53
54 DTRTI2 computes the inverse of a real upper or lower triangular
55 matrix.
56
57 This is the Level 2 BLAS version of the algorithm.
58
59 Arguments
60 =========
61
62 UPLO (input) CHARACTER*1
63 Specifies whether the matrix A is upper or lower triangular.
64 = 'U': Upper triangular
65 = 'L': Lower triangular
66
67 DIAG (input) CHARACTER*1
68 Specifies whether or not the matrix A is unit triangular.
69 = 'N': Non-unit triangular
70 = 'U': Unit triangular
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 triangular matrix A. If UPLO = 'U', the
77 leading n by n upper triangular part of the array A contains
78 the upper triangular matrix, 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 the array A contains
81 the lower triangular matrix, and the strictly upper
82 triangular part of A is not referenced. If DIAG = 'U', the
83 diagonal elements of A are also not referenced and are
84 assumed to be 1.
85
86 On exit, the (triangular) inverse of the original matrix, in
87 the same storage format.
88
89 LDA (input) INTEGER
90 The leading dimension of the array A. LDA >= max(1,N).
91
92 INFO (output) INTEGER
93 = 0: successful exit
94 < 0: if INFO = -k, the k-th argument had an illegal value
95
96 =====================================================================
97
98
99 Test the input parameters.
100
101 Parameter adjustments */
102 /* Table of constant values */
103 integer c__1 = 1;
104
105 /* System generated locals */
106 integer a_dim1, a_offset, i__1, i__2;
107 /* Local variables */
108 integer j;
109 logical upper;
110 logical nounit;
111 Treal ajj;
112#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
113
114
115 a_dim1 = *lda;
116 a_offset = 1 + a_dim1 * 1;
117 a -= a_offset;
118
119 /* Function Body */
120 *info = 0;
121 upper = template_blas_lsame(uplo, "U");
122 nounit = template_blas_lsame(diag, "N");
123 if (! upper && ! template_blas_lsame(uplo, "L")) {
124 *info = -1;
125 } else if (! nounit && ! template_blas_lsame(diag, "U")) {
126 *info = -2;
127 } else if (*n < 0) {
128 *info = -3;
129 } else if (*lda < maxMACRO(1,*n)) {
130 *info = -5;
131 }
132 if (*info != 0) {
133 i__1 = -(*info);
134 template_blas_erbla("TRTI2 ", &i__1);
135 return 0;
136 }
137
138 if (upper) {
139
140/* Compute inverse of upper triangular matrix. */
141
142 i__1 = *n;
143 for (j = 1; j <= i__1; ++j) {
144 if (nounit) {
145 a_ref(j, j) = 1. / a_ref(j, j);
146 ajj = -a_ref(j, j);
147 } else {
148 ajj = -1.;
149 }
150
151/* Compute elements 1:j-1 of j-th column. */
152
153 i__2 = j - 1;
154 template_blas_trmv("Upper", "No transpose", diag, &i__2, &a[a_offset], lda, &
155 a_ref(1, j), &c__1);
156 i__2 = j - 1;
157 template_blas_scal(&i__2, &ajj, &a_ref(1, j), &c__1);
158/* L10: */
159 }
160 } else {
161
162/* Compute inverse of lower triangular matrix. */
163
164 for (j = *n; j >= 1; --j) {
165 if (nounit) {
166 a_ref(j, j) = 1. / a_ref(j, j);
167 ajj = -a_ref(j, j);
168 } else {
169 ajj = -1.;
170 }
171 if (j < *n) {
172
173/* Compute elements j+1:n of j-th column. */
174
175 i__1 = *n - j;
176 template_blas_trmv("Lower", "No transpose", diag, &i__1, &a_ref(j + 1, j
177 + 1), lda, &a_ref(j + 1, j), &c__1);
178 i__1 = *n - j;
179 template_blas_scal(&i__1, &ajj, &a_ref(j + 1, j), &c__1);
180 }
181/* L20: */
182 }
183 }
184
185 return 0;
186
187/* End of DTRTI2 */
188
189} /* dtrti2_ */
190
191#undef a_ref
192
193
194#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
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
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_trmv(const char *uplo, const char *trans, const char *diag, const integer *n, const Treal *a, const integer *lda, Treal *x, const integer *incx)
Definition: template_blas_trmv.h:42
#define a_ref(a_1, a_2)
int template_lapack_trti2(const char *uplo, const char *diag, const integer *n, Treal *a, const integer *lda, integer *info)
Definition: template_lapack_trti2.h:42