ergo
template_lapack_trtri.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_TRTRI_HEADER
38#define TEMPLATE_LAPACK_TRTRI_HEADER
39
40
41template<class Treal>
42int template_lapack_trtri(const char *uplo, const char *diag,
43 const integer *n, Treal *a,
44 const integer *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 DTRTRI computes the inverse of a real upper or lower triangular
56 matrix A.
57
58 This is the Level 3 BLAS version of the algorithm.
59
60 Arguments
61 =========
62
63 UPLO (input) CHARACTER*1
64 = 'U': A is upper triangular;
65 = 'L': A is lower triangular.
66
67 DIAG (input) CHARACTER*1
68 = 'N': A is non-unit triangular;
69 = 'U': A is unit triangular.
70
71 N (input) INTEGER
72 The order of the matrix A. N >= 0.
73
74 A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
75 On entry, the triangular matrix A. If UPLO = 'U', the
76 leading N-by-N upper triangular part of the array A contains
77 the upper triangular matrix, and the strictly lower
78 triangular part of A is not referenced. If UPLO = 'L', the
79 leading N-by-N lower triangular part of the array A contains
80 the lower triangular matrix, and the strictly upper
81 triangular part of A is not referenced. If DIAG = 'U', the
82 diagonal elements of A are also not referenced and are
83 assumed to be 1.
84 On exit, the (triangular) inverse of the original matrix, in
85 the same storage format.
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, A(i,i) is exactly zero. The triangular
94 matrix is singular and its inverse can not be computed.
95
96 =====================================================================
97
98
99 Test the input parameters.
100
101 Parameter adjustments */
102 /* Table of constant values */
103 integer c__1 = 1;
104 integer c_n1 = -1;
105 integer c__2 = 2;
106 Treal c_b18 = 1.;
107 Treal c_b22 = -1.;
108
109 /* System generated locals */
110 address a__1[2];
111 integer a_dim1, a_offset, i__1, i__2[2], i__3, i__4, i__5;
112 char ch__1[2];
113 /* Local variables */
114 integer j;
115 logical upper;
116 integer jb, nb, nn;
117 logical nounit;
118#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
119
120
121 a_dim1 = *lda;
122 a_offset = 1 + a_dim1 * 1;
123 a -= a_offset;
124
125 /* Function Body */
126 *info = 0;
127 upper = template_blas_lsame(uplo, "U");
128 nounit = template_blas_lsame(diag, "N");
129 if (! upper && ! template_blas_lsame(uplo, "L")) {
130 *info = -1;
131 } else if (! nounit && ! template_blas_lsame(diag, "U")) {
132 *info = -2;
133 } else if (*n < 0) {
134 *info = -3;
135 } else if (*lda < maxMACRO(1,*n)) {
136 *info = -5;
137 }
138 if (*info != 0) {
139 i__1 = -(*info);
140 template_blas_erbla("TRTRI ", &i__1);
141 return 0;
142 }
143
144/* Quick return if possible */
145
146 if (*n == 0) {
147 return 0;
148 }
149
150/* Check for singularity if non-unit. */
151
152 if (nounit) {
153 i__1 = *n;
154 for (*info = 1; *info <= i__1; ++(*info)) {
155 if (a_ref(*info, *info) == 0.) {
156 return 0;
157 }
158/* L10: */
159 }
160 *info = 0;
161 }
162
163/* Determine the block size for this environment.
164
165 Writing concatenation */
166 i__2[0] = 1, a__1[0] = (char*)uplo;
167 i__2[1] = 1, a__1[1] = (char*)diag;
168 template_blas_s_cat(ch__1, a__1, i__2, &c__2, (ftnlen)2);
169 nb = template_lapack_ilaenv(&c__1, "DTRTRI", ch__1, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (
170 ftnlen)2);
171 if (nb <= 1 || nb >= *n) {
172
173/* Use unblocked code */
174
175 template_lapack_trti2(uplo, diag, n, &a[a_offset], lda, info);
176 } else {
177
178/* Use blocked code */
179
180 if (upper) {
181
182/* Compute inverse of upper triangular matrix */
183
184 i__1 = *n;
185 i__3 = nb;
186 for (j = 1; i__3 < 0 ? j >= i__1 : j <= i__1; j += i__3) {
187/* Computing MIN */
188 i__4 = nb, i__5 = *n - j + 1;
189 jb = minMACRO(i__4,i__5);
190
191/* Compute rows 1:j-1 of current block column */
192
193 i__4 = j - 1;
194 template_blas_trmm("Left", "Upper", "No transpose", diag, &i__4, &jb, &
195 c_b18, &a[a_offset], lda, &a_ref(1, j), lda);
196 i__4 = j - 1;
197 template_blas_trsm("Right", "Upper", "No transpose", diag, &i__4, &jb, &
198 c_b22, &a_ref(j, j), lda, &a_ref(1, j), lda);
199
200/* Compute inverse of current diagonal block */
201
202 template_lapack_trti2("Upper", diag, &jb, &a_ref(j, j), lda, info);
203/* L20: */
204 }
205 } else {
206
207/* Compute inverse of lower triangular matrix */
208
209 nn = (*n - 1) / nb * nb + 1;
210 i__3 = -nb;
211 for (j = nn; i__3 < 0 ? j >= 1 : j <= 1; j += i__3) {
212/* Computing MIN */
213 i__1 = nb, i__4 = *n - j + 1;
214 jb = minMACRO(i__1,i__4);
215 if (j + jb <= *n) {
216
217/* Compute rows j+jb:n of current block column */
218
219 i__1 = *n - j - jb + 1;
220 template_blas_trmm("Left", "Lower", "No transpose", diag, &i__1, &jb,
221 &c_b18, &a_ref(j + jb, j + jb), lda, &a_ref(j +
222 jb, j), lda);
223 i__1 = *n - j - jb + 1;
224 template_blas_trsm("Right", "Lower", "No transpose", diag, &i__1, &jb,
225 &c_b22, &a_ref(j, j), lda, &a_ref(j + jb, j),
226 lda);
227 }
228
229/* Compute inverse of current diagonal block */
230
231 template_lapack_trti2("Lower", diag, &jb, &a_ref(j, j), lda, info);
232/* L30: */
233 }
234 }
235 }
236
237 return 0;
238
239/* End of DTRTRI */
240
241} /* dtrtri_ */
242
243#undef a_ref
244
245
246#endif
void template_blas_s_cat(char *lp, char *rpp[], ftnlen rnp[], ftnlen *np, ftnlen ll)
Definition: template_blas_common.cc:204
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
char * address
Definition: template_blas_common.h:43
int integer
Definition: template_blas_common.h:40
int ftnlen
Definition: template_blas_common.h:42
#define minMACRO(a, b)
Definition: template_blas_common.h:46
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
bool logical
Definition: template_blas_common.h:41
int template_blas_trmm(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_trmm.h:43
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
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
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
#define a_ref(a_1, a_2)
int template_lapack_trtri(const char *uplo, const char *diag, const integer *n, Treal *a, const integer *lda, integer *info)
Definition: template_lapack_trtri.h:42