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