ergo
template_lapack_tptri.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_TPTRI_HEADER
38 #define TEMPLATE_LAPACK_TPTRI_HEADER
39 
40 #include "template_lapack_common.h"
41 
42 template<class Treal>
43 int template_lapack_tptri(const char *uplo, const char *diag, const integer *n, Treal *
44  ap, 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  September 30, 1994
50 
51 
52  Purpose
53  =======
54 
55  DTPTRI computes the inverse of a real upper or lower triangular
56  matrix A stored in packed format.
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  AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
73  On entry, the upper or lower triangular matrix A, stored
74  columnwise in a linear array. The j-th column of A is stored
75  in the array AP as follows:
76  if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
77  if UPLO = 'L', AP(i + (j-1)*((2*n-j)/2) = A(i,j) for j<=i<=n.
78  See below for further details.
79  On exit, the (triangular) inverse of the original matrix, in
80  the same packed storage format.
81 
82  INFO (output) INTEGER
83  = 0: successful exit
84  < 0: if INFO = -i, the i-th argument had an illegal value
85  > 0: if INFO = i, A(i,i) is exactly zero. The triangular
86  matrix is singular and its inverse can not be computed.
87 
88  Further Details
89  ===============
90 
91  A triangular matrix A can be transferred to packed storage using one
92  of the following program segments:
93 
94  UPLO = 'U': UPLO = 'L':
95 
96  JC = 1 JC = 1
97  DO 2 J = 1, N DO 2 J = 1, N
98  DO 1 I = 1, J DO 1 I = J, N
99  AP(JC+I-1) = A(I,J) AP(JC+I-J) = A(I,J)
100  1 CONTINUE 1 CONTINUE
101  JC = JC + J JC = JC + N - J + 1
102  2 CONTINUE 2 CONTINUE
103 
104  =====================================================================
105 
106 
107  Test the input parameters.
108 
109  Parameter adjustments */
110  /* Table of constant values */
111  integer c__1 = 1;
112 
113  /* System generated locals */
114  integer i__1, i__2;
115  /* Local variables */
116  integer j;
117  logical upper;
118  integer jc, jj;
119  integer jclast;
120  logical nounit;
121  Treal ajj;
122 
123 
124  --ap;
125 
126  /* Initialization added by Elias to get rid of compiler warnings. */
127  jclast = 0;
128  /* Function Body */
129  *info = 0;
130  upper = template_blas_lsame(uplo, "U");
131  nounit = template_blas_lsame(diag, "N");
132  if (! upper && ! template_blas_lsame(uplo, "L")) {
133  *info = -1;
134  } else if (! nounit && ! template_blas_lsame(diag, "U")) {
135  *info = -2;
136  } else if (*n < 0) {
137  *info = -3;
138  }
139  if (*info != 0) {
140  i__1 = -(*info);
141  template_blas_erbla("TPTRI ", &i__1);
142  return 0;
143  }
144 
145 /* Check for singularity if non-unit. */
146 
147  if (nounit) {
148  if (upper) {
149  jj = 0;
150  i__1 = *n;
151  for (*info = 1; *info <= i__1; ++(*info)) {
152  jj += *info;
153  if (ap[jj] == 0.) {
154  return 0;
155  }
156 /* L10: */
157  }
158  } else {
159  jj = 1;
160  i__1 = *n;
161  for (*info = 1; *info <= i__1; ++(*info)) {
162  if (ap[jj] == 0.) {
163  return 0;
164  }
165  jj = jj + *n - *info + 1;
166 /* L20: */
167  }
168  }
169  *info = 0;
170  }
171 
172  if (upper) {
173 
174 /* Compute inverse of upper triangular matrix. */
175 
176  jc = 1;
177  i__1 = *n;
178  for (j = 1; j <= i__1; ++j) {
179  if (nounit) {
180  ap[jc + j - 1] = 1. / ap[jc + j - 1];
181  ajj = -ap[jc + j - 1];
182  } else {
183  ajj = -1.;
184  }
185 
186 /* Compute elements 1:j-1 of j-th column. */
187 
188  i__2 = j - 1;
189  template_blas_tpmv("Upper", "No transpose", diag, &i__2, &ap[1], &ap[jc], &
190  c__1);
191  i__2 = j - 1;
192  template_blas_scal(&i__2, &ajj, &ap[jc], &c__1);
193  jc += j;
194 /* L30: */
195  }
196 
197  } else {
198 
199 /* Compute inverse of lower triangular matrix. */
200 
201  jc = *n * (*n + 1) / 2;
202  for (j = *n; j >= 1; --j) {
203  if (nounit) {
204  ap[jc] = 1. / ap[jc];
205  ajj = -ap[jc];
206  } else {
207  ajj = -1.;
208  }
209  if (j < *n) {
210 
211 /* Compute elements j+1:n of j-th column. */
212 
213  i__1 = *n - j;
214  template_blas_tpmv("Lower", "No transpose", diag, &i__1, &ap[jclast], &ap[
215  jc + 1], &c__1);
216  i__1 = *n - j;
217  template_blas_scal(&i__1, &ajj, &ap[jc + 1], &c__1);
218  }
219  jclast = jc;
220  jc = jc - *n + j - 2;
221 /* L40: */
222  }
223  }
224 
225  return 0;
226 
227 /* End of DTPTRI */
228 
229 } /* dtptri_ */
230 
231 #endif
template_blas_tpmv
int template_blas_tpmv(const char *uplo, const char *trans, const char *diag, const integer *n, const Treal *ap, Treal *x, const integer *incx)
Definition: template_blas_tpmv.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
template_lapack_tptri
int template_lapack_tptri(const char *uplo, const char *diag, const integer *n, Treal *ap, integer *info)
Definition: template_lapack_tptri.h:43
logical
bool logical
Definition: template_blas_common.h:41
template_lapack_common.h
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