ergo
template_lapack_trti2.h
Go to the documentation of this file.
1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2019 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 
41 template<class Treal>
42 int 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
template_lapack_trti2
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
template_blas_scal
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition: template_blas_scal.h:43
logical
bool logical
Definition: template_blas_common.h:41
a_ref
#define a_ref(a_1, a_2)
template_blas_erbla
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
template_blas_lsame
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
integer
int integer
Definition: template_blas_common.h:40
template_blas_trmv
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
maxMACRO
#define maxMACRO(a, b)
Definition: template_blas_common.h:45