ergo
template_lapack_trtri.h
Go to the documentation of this file.
1 /* Ergo, version 3.4, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2014 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28  /* This file belongs to the template_lapack part of the Ergo source
29  * code. The source files in the template_lapack directory are modified
30  * versions of files originally distributed as CLAPACK, see the
31  * Copyright/license notice in the file template_lapack/COPYING.
32  */
33 
34 
35 #ifndef TEMPLATE_LAPACK_TRTRI_HEADER
36 #define TEMPLATE_LAPACK_TRTRI_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_trtri(const char *uplo, const char *diag,
41  const integer *n, Treal *a,
42  const integer *lda, integer *info)
43 {
44 /* -- LAPACK routine (version 3.0) --
45  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
46  Courant Institute, Argonne National Lab, and Rice University
47  March 31, 1993
48 
49 
50  Purpose
51  =======
52 
53  DTRTRI computes the inverse of a real upper or lower triangular
54  matrix A.
55 
56  This is the Level 3 BLAS version of the algorithm.
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  A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
73  On entry, the triangular matrix A. If UPLO = 'U', the
74  leading N-by-N upper triangular part of the array A contains
75  the upper triangular matrix, and the strictly lower
76  triangular part of A is not referenced. If UPLO = 'L', the
77  leading N-by-N lower triangular part of the array A contains
78  the lower triangular matrix, and the strictly upper
79  triangular part of A is not referenced. If DIAG = 'U', the
80  diagonal elements of A are also not referenced and are
81  assumed to be 1.
82  On exit, the (triangular) inverse of the original matrix, in
83  the same storage format.
84 
85  LDA (input) INTEGER
86  The leading dimension of the array A. LDA >= max(1,N).
87 
88  INFO (output) INTEGER
89  = 0: successful exit
90  < 0: if INFO = -i, the i-th argument had an illegal value
91  > 0: if INFO = i, A(i,i) is exactly zero. The triangular
92  matrix is singular and its inverse can not be computed.
93 
94  =====================================================================
95 
96 
97  Test the input parameters.
98 
99  Parameter adjustments */
100  /* Table of constant values */
101  integer c__1 = 1;
102  integer c_n1 = -1;
103  integer c__2 = 2;
104  Treal c_b18 = 1.;
105  Treal c_b22 = -1.;
106 
107  /* System generated locals */
108  address a__1[2];
109  integer a_dim1, a_offset, i__1, i__2[2], i__3, i__4, i__5;
110  char ch__1[2];
111  /* Local variables */
112  integer j;
113  logical upper;
114  integer jb, nb, nn;
115  logical nounit;
116 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
117 
118 
119  a_dim1 = *lda;
120  a_offset = 1 + a_dim1 * 1;
121  a -= a_offset;
122 
123  /* Function Body */
124  *info = 0;
125  upper = template_blas_lsame(uplo, "U");
126  nounit = template_blas_lsame(diag, "N");
127  if (! upper && ! template_blas_lsame(uplo, "L")) {
128  *info = -1;
129  } else if (! nounit && ! template_blas_lsame(diag, "U")) {
130  *info = -2;
131  } else if (*n < 0) {
132  *info = -3;
133  } else if (*lda < maxMACRO(1,*n)) {
134  *info = -5;
135  }
136  if (*info != 0) {
137  i__1 = -(*info);
138  template_blas_erbla("TRTRI ", &i__1);
139  return 0;
140  }
141 
142 /* Quick return if possible */
143 
144  if (*n == 0) {
145  return 0;
146  }
147 
148 /* Check for singularity if non-unit. */
149 
150  if (nounit) {
151  i__1 = *n;
152  for (*info = 1; *info <= i__1; ++(*info)) {
153  if (a_ref(*info, *info) == 0.) {
154  return 0;
155  }
156 /* L10: */
157  }
158  *info = 0;
159  }
160 
161 /* Determine the block size for this environment.
162 
163  Writing concatenation */
164  i__2[0] = 1, a__1[0] = (char*)uplo;
165  i__2[1] = 1, a__1[1] = (char*)diag;
166  template_blas_s_cat(ch__1, a__1, i__2, &c__2, (ftnlen)2);
167  nb = template_lapack_ilaenv(&c__1, "DTRTRI", ch__1, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (
168  ftnlen)2);
169  if (nb <= 1 || nb >= *n) {
170 
171 /* Use unblocked code */
172 
173  template_lapack_trti2(uplo, diag, n, &a[a_offset], lda, info);
174  } else {
175 
176 /* Use blocked code */
177 
178  if (upper) {
179 
180 /* Compute inverse of upper triangular matrix */
181 
182  i__1 = *n;
183  i__3 = nb;
184  for (j = 1; i__3 < 0 ? j >= i__1 : j <= i__1; j += i__3) {
185 /* Computing MIN */
186  i__4 = nb, i__5 = *n - j + 1;
187  jb = minMACRO(i__4,i__5);
188 
189 /* Compute rows 1:j-1 of current block column */
190 
191  i__4 = j - 1;
192  template_blas_trmm("Left", "Upper", "No transpose", diag, &i__4, &jb, &
193  c_b18, &a[a_offset], lda, &a_ref(1, j), lda);
194  i__4 = j - 1;
195  template_blas_trsm("Right", "Upper", "No transpose", diag, &i__4, &jb, &
196  c_b22, &a_ref(j, j), lda, &a_ref(1, j), lda);
197 
198 /* Compute inverse of current diagonal block */
199 
200  template_lapack_trti2("Upper", diag, &jb, &a_ref(j, j), lda, info);
201 /* L20: */
202  }
203  } else {
204 
205 /* Compute inverse of lower triangular matrix */
206 
207  nn = (*n - 1) / nb * nb + 1;
208  i__3 = -nb;
209  for (j = nn; i__3 < 0 ? j >= 1 : j <= 1; j += i__3) {
210 /* Computing MIN */
211  i__1 = nb, i__4 = *n - j + 1;
212  jb = minMACRO(i__1,i__4);
213  if (j + jb <= *n) {
214 
215 /* Compute rows j+jb:n of current block column */
216 
217  i__1 = *n - j - jb + 1;
218  template_blas_trmm("Left", "Lower", "No transpose", diag, &i__1, &jb,
219  &c_b18, &a_ref(j + jb, j + jb), lda, &a_ref(j +
220  jb, j), lda);
221  i__1 = *n - j - jb + 1;
222  template_blas_trsm("Right", "Lower", "No transpose", diag, &i__1, &jb,
223  &c_b22, &a_ref(j, j), lda, &a_ref(j + jb, j),
224  lda);
225  }
226 
227 /* Compute inverse of current diagonal block */
228 
229  template_lapack_trti2("Lower", diag, &jb, &a_ref(j, j), lda, info);
230 /* L30: */
231  }
232  }
233  }
234 
235  return 0;
236 
237 /* End of DTRTRI */
238 
239 } /* dtrtri_ */
240 
241 #undef a_ref
242 
243 
244 #endif
#define a_ref(a_1, a_2)
int integer
Definition: template_blas_common.h:38
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:41
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:279
char * address
Definition: template_blas_common.h:41
#define maxMACRO(a, b)
Definition: template_blas_common.h:43
#define minMACRO(a, b)
Definition: template_blas_common.h:44
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
void template_blas_s_cat(char *lp, char *rpp[], ftnlen rnp[], ftnlen *np, ftnlen ll)
Definition: template_blas_common.cc:202
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:40
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:41
bool logical
Definition: template_blas_common.h:39
int ftnlen
Definition: template_blas_common.h:40
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:44
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:40