ergo
template_lapack_trtri.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_TRTRI_HEADER
38 #define TEMPLATE_LAPACK_TRTRI_HEADER
39 
40 
41 template<class Treal>
42 int 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
ftnlen
int ftnlen
Definition: template_blas_common.h:42
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_lapack_ilaenv
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
minMACRO
#define minMACRO(a, b)
Definition: template_blas_common.h:46
template_blas_trsm
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
template_lapack_trtri
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
logical
bool logical
Definition: template_blas_common.h:41
template_blas_trmm
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
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
template_blas_s_cat
void template_blas_s_cat(char *lp, char *rpp[], ftnlen rnp[], ftnlen *np, ftnlen ll)
Definition: template_blas_common.cc:204
integer
int integer
Definition: template_blas_common.h:40
address
char * address
Definition: template_blas_common.h:43
a_ref
#define a_ref(a_1, a_2)
maxMACRO
#define maxMACRO(a, b)
Definition: template_blas_common.h:45