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_TRTRI_HEADER 00036 #define TEMPLATE_LAPACK_TRTRI_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_trtri(const char *uplo, const char *diag, 00041 const integer *n, Treal *a, 00042 const integer *lda, integer *info) 00043 { 00044 /* -- LAPACK routine (version 3.0) -- 00045 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00046 Courant Institute, Argonne National Lab, and Rice University 00047 March 31, 1993 00048 00049 00050 Purpose 00051 ======= 00052 00053 DTRTRI computes the inverse of a real upper or lower triangular 00054 matrix A. 00055 00056 This is the Level 3 BLAS version of the algorithm. 00057 00058 Arguments 00059 ========= 00060 00061 UPLO (input) CHARACTER*1 00062 = 'U': A is upper triangular; 00063 = 'L': A is lower triangular. 00064 00065 DIAG (input) CHARACTER*1 00066 = 'N': A is non-unit triangular; 00067 = 'U': A is unit triangular. 00068 00069 N (input) INTEGER 00070 The order of the matrix A. N >= 0. 00071 00072 A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00073 On entry, the triangular matrix A. If UPLO = 'U', the 00074 leading N-by-N upper triangular part of the array A contains 00075 the upper triangular matrix, and the strictly lower 00076 triangular part of A is not referenced. If UPLO = 'L', the 00077 leading N-by-N lower triangular part of the array A contains 00078 the lower triangular matrix, and the strictly upper 00079 triangular part of A is not referenced. If DIAG = 'U', the 00080 diagonal elements of A are also not referenced and are 00081 assumed to be 1. 00082 On exit, the (triangular) inverse of the original matrix, in 00083 the same storage format. 00084 00085 LDA (input) INTEGER 00086 The leading dimension of the array A. LDA >= max(1,N). 00087 00088 INFO (output) INTEGER 00089 = 0: successful exit 00090 < 0: if INFO = -i, the i-th argument had an illegal value 00091 > 0: if INFO = i, A(i,i) is exactly zero. The triangular 00092 matrix is singular and its inverse can not be computed. 00093 00094 ===================================================================== 00095 00096 00097 Test the input parameters. 00098 00099 Parameter adjustments */ 00100 /* Table of constant values */ 00101 integer c__1 = 1; 00102 integer c_n1 = -1; 00103 integer c__2 = 2; 00104 Treal c_b18 = 1.; 00105 Treal c_b22 = -1.; 00106 00107 /* System generated locals */ 00108 address a__1[2]; 00109 integer a_dim1, a_offset, i__1, i__2[2], i__3, i__4, i__5; 00110 char ch__1[2]; 00111 /* Local variables */ 00112 integer j; 00113 logical upper; 00114 integer jb, nb, nn; 00115 logical nounit; 00116 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 00117 00118 00119 a_dim1 = *lda; 00120 a_offset = 1 + a_dim1 * 1; 00121 a -= a_offset; 00122 00123 /* Function Body */ 00124 *info = 0; 00125 upper = template_blas_lsame(uplo, "U"); 00126 nounit = template_blas_lsame(diag, "N"); 00127 if (! upper && ! template_blas_lsame(uplo, "L")) { 00128 *info = -1; 00129 } else if (! nounit && ! template_blas_lsame(diag, "U")) { 00130 *info = -2; 00131 } else if (*n < 0) { 00132 *info = -3; 00133 } else if (*lda < maxMACRO(1,*n)) { 00134 *info = -5; 00135 } 00136 if (*info != 0) { 00137 i__1 = -(*info); 00138 template_blas_erbla("TRTRI ", &i__1); 00139 return 0; 00140 } 00141 00142 /* Quick return if possible */ 00143 00144 if (*n == 0) { 00145 return 0; 00146 } 00147 00148 /* Check for singularity if non-unit. */ 00149 00150 if (nounit) { 00151 i__1 = *n; 00152 for (*info = 1; *info <= i__1; ++(*info)) { 00153 if (a_ref(*info, *info) == 0.) { 00154 return 0; 00155 } 00156 /* L10: */ 00157 } 00158 *info = 0; 00159 } 00160 00161 /* Determine the block size for this environment. 00162 00163 Writing concatenation */ 00164 i__2[0] = 1, a__1[0] = (char*)uplo; 00165 i__2[1] = 1, a__1[1] = (char*)diag; 00166 template_blas_s_cat(ch__1, a__1, i__2, &c__2, (ftnlen)2); 00167 nb = template_lapack_ilaenv(&c__1, "DTRTRI", ch__1, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, ( 00168 ftnlen)2); 00169 if (nb <= 1 || nb >= *n) { 00170 00171 /* Use unblocked code */ 00172 00173 template_lapack_trti2(uplo, diag, n, &a[a_offset], lda, info); 00174 } else { 00175 00176 /* Use blocked code */ 00177 00178 if (upper) { 00179 00180 /* Compute inverse of upper triangular matrix */ 00181 00182 i__1 = *n; 00183 i__3 = nb; 00184 for (j = 1; i__3 < 0 ? j >= i__1 : j <= i__1; j += i__3) { 00185 /* Computing MIN */ 00186 i__4 = nb, i__5 = *n - j + 1; 00187 jb = minMACRO(i__4,i__5); 00188 00189 /* Compute rows 1:j-1 of current block column */ 00190 00191 i__4 = j - 1; 00192 template_blas_trmm("Left", "Upper", "No transpose", diag, &i__4, &jb, & 00193 c_b18, &a[a_offset], lda, &a_ref(1, j), lda); 00194 i__4 = j - 1; 00195 template_blas_trsm("Right", "Upper", "No transpose", diag, &i__4, &jb, & 00196 c_b22, &a_ref(j, j), lda, &a_ref(1, j), lda); 00197 00198 /* Compute inverse of current diagonal block */ 00199 00200 template_lapack_trti2("Upper", diag, &jb, &a_ref(j, j), lda, info); 00201 /* L20: */ 00202 } 00203 } else { 00204 00205 /* Compute inverse of lower triangular matrix */ 00206 00207 nn = (*n - 1) / nb * nb + 1; 00208 i__3 = -nb; 00209 for (j = nn; i__3 < 0 ? j >= 1 : j <= 1; j += i__3) { 00210 /* Computing MIN */ 00211 i__1 = nb, i__4 = *n - j + 1; 00212 jb = minMACRO(i__1,i__4); 00213 if (j + jb <= *n) { 00214 00215 /* Compute rows j+jb:n of current block column */ 00216 00217 i__1 = *n - j - jb + 1; 00218 template_blas_trmm("Left", "Lower", "No transpose", diag, &i__1, &jb, 00219 &c_b18, &a_ref(j + jb, j + jb), lda, &a_ref(j + 00220 jb, j), lda); 00221 i__1 = *n - j - jb + 1; 00222 template_blas_trsm("Right", "Lower", "No transpose", diag, &i__1, &jb, 00223 &c_b22, &a_ref(j, j), lda, &a_ref(j + jb, j), 00224 lda); 00225 } 00226 00227 /* Compute inverse of current diagonal block */ 00228 00229 template_lapack_trti2("Lower", diag, &jb, &a_ref(j, j), lda, info); 00230 /* L30: */ 00231 } 00232 } 00233 } 00234 00235 return 0; 00236 00237 /* End of DTRTRI */ 00238 00239 } /* dtrtri_ */ 00240 00241 #undef a_ref 00242 00243 00244 #endif