00001 /* Ergo, version 3.7, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2018 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, 00004 * and Anastasia Kruchinina. 00005 * 00006 * This program is free software: you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation, either version 3 of the License, or 00009 * (at your option) any later version. 00010 * 00011 * This program is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00018 * 00019 * Primary academic reference: 00020 * Ergo: An open-source program for linear-scaling electronic structure 00021 * calculations, 00022 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia 00023 * Kruchinina, 00024 * SoftwareX 7, 107 (2018), 00025 * <http://dx.doi.org/10.1016/j.softx.2018.03.005> 00026 * 00027 * For further information about Ergo, see <http://www.ergoscf.org>. 00028 */ 00029 00030 /* This file belongs to the template_lapack part of the Ergo source 00031 * code. The source files in the template_lapack directory are modified 00032 * versions of files originally distributed as CLAPACK, see the 00033 * Copyright/license notice in the file template_lapack/COPYING. 00034 */ 00035 00036 00037 #ifndef TEMPLATE_LAPACK_TRTI2_HEADER 00038 #define TEMPLATE_LAPACK_TRTI2_HEADER 00039 00040 00041 template<class Treal> 00042 int template_lapack_trti2(const char *uplo, const char *diag, const integer *n, Treal * 00043 a, const integer *lda, integer *info) 00044 { 00045 /* -- LAPACK routine (version 3.0) -- 00046 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00047 Courant Institute, Argonne National Lab, and Rice University 00048 February 29, 1992 00049 00050 00051 Purpose 00052 ======= 00053 00054 DTRTI2 computes the inverse of a real upper or lower triangular 00055 matrix. 00056 00057 This is the Level 2 BLAS version of the algorithm. 00058 00059 Arguments 00060 ========= 00061 00062 UPLO (input) CHARACTER*1 00063 Specifies whether the matrix A is upper or lower triangular. 00064 = 'U': Upper triangular 00065 = 'L': Lower triangular 00066 00067 DIAG (input) CHARACTER*1 00068 Specifies whether or not the matrix A is unit triangular. 00069 = 'N': Non-unit triangular 00070 = 'U': Unit triangular 00071 00072 N (input) INTEGER 00073 The order of the matrix A. N >= 0. 00074 00075 A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00076 On entry, the triangular matrix A. If UPLO = 'U', the 00077 leading n by n upper triangular part of the array A contains 00078 the upper triangular matrix, and the strictly lower 00079 triangular part of A is not referenced. If UPLO = 'L', the 00080 leading n by n lower triangular part of the array A contains 00081 the lower triangular matrix, and the strictly upper 00082 triangular part of A is not referenced. If DIAG = 'U', the 00083 diagonal elements of A are also not referenced and are 00084 assumed to be 1. 00085 00086 On exit, the (triangular) inverse of the original matrix, in 00087 the same storage format. 00088 00089 LDA (input) INTEGER 00090 The leading dimension of the array A. LDA >= max(1,N). 00091 00092 INFO (output) INTEGER 00093 = 0: successful exit 00094 < 0: if INFO = -k, the k-th argument had an illegal value 00095 00096 ===================================================================== 00097 00098 00099 Test the input parameters. 00100 00101 Parameter adjustments */ 00102 /* Table of constant values */ 00103 integer c__1 = 1; 00104 00105 /* System generated locals */ 00106 integer a_dim1, a_offset, i__1, i__2; 00107 /* Local variables */ 00108 integer j; 00109 logical upper; 00110 logical nounit; 00111 Treal ajj; 00112 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 00113 00114 00115 a_dim1 = *lda; 00116 a_offset = 1 + a_dim1 * 1; 00117 a -= a_offset; 00118 00119 /* Function Body */ 00120 *info = 0; 00121 upper = template_blas_lsame(uplo, "U"); 00122 nounit = template_blas_lsame(diag, "N"); 00123 if (! upper && ! template_blas_lsame(uplo, "L")) { 00124 *info = -1; 00125 } else if (! nounit && ! template_blas_lsame(diag, "U")) { 00126 *info = -2; 00127 } else if (*n < 0) { 00128 *info = -3; 00129 } else if (*lda < maxMACRO(1,*n)) { 00130 *info = -5; 00131 } 00132 if (*info != 0) { 00133 i__1 = -(*info); 00134 template_blas_erbla("TRTI2 ", &i__1); 00135 return 0; 00136 } 00137 00138 if (upper) { 00139 00140 /* Compute inverse of upper triangular matrix. */ 00141 00142 i__1 = *n; 00143 for (j = 1; j <= i__1; ++j) { 00144 if (nounit) { 00145 a_ref(j, j) = 1. / a_ref(j, j); 00146 ajj = -a_ref(j, j); 00147 } else { 00148 ajj = -1.; 00149 } 00150 00151 /* Compute elements 1:j-1 of j-th column. */ 00152 00153 i__2 = j - 1; 00154 template_blas_trmv("Upper", "No transpose", diag, &i__2, &a[a_offset], lda, & 00155 a_ref(1, j), &c__1); 00156 i__2 = j - 1; 00157 template_blas_scal(&i__2, &ajj, &a_ref(1, j), &c__1); 00158 /* L10: */ 00159 } 00160 } else { 00161 00162 /* Compute inverse of lower triangular matrix. */ 00163 00164 for (j = *n; j >= 1; --j) { 00165 if (nounit) { 00166 a_ref(j, j) = 1. / a_ref(j, j); 00167 ajj = -a_ref(j, j); 00168 } else { 00169 ajj = -1.; 00170 } 00171 if (j < *n) { 00172 00173 /* Compute elements j+1:n of j-th column. */ 00174 00175 i__1 = *n - j; 00176 template_blas_trmv("Lower", "No transpose", diag, &i__1, &a_ref(j + 1, j 00177 + 1), lda, &a_ref(j + 1, j), &c__1); 00178 i__1 = *n - j; 00179 template_blas_scal(&i__1, &ajj, &a_ref(j + 1, j), &c__1); 00180 } 00181 /* L20: */ 00182 } 00183 } 00184 00185 return 0; 00186 00187 /* End of DTRTI2 */ 00188 00189 } /* dtrti2_ */ 00190 00191 #undef a_ref 00192 00193 00194 #endif