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