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