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_BLAS_TPSV_HEADER 00036 #define TEMPLATE_BLAS_TPSV_HEADER 00037 00038 00039 template<class Treal> 00040 int template_blas_tpsv(const char *uplo, const char *trans, const char *diag, const integer *n, 00041 const Treal *ap, Treal *x, const integer *incx) 00042 { 00043 /* System generated locals */ 00044 integer i__1, i__2; 00045 /* Local variables */ 00046 integer info; 00047 Treal temp; 00048 integer i__, j, k; 00049 integer kk, ix, jx, kx; 00050 logical nounit; 00051 /* Purpose 00052 ======= 00053 DTPSV solves one of the systems of equations 00054 A*x = b, or A'*x = b, 00055 where b and x are n element vectors and A is an n by n unit, or 00056 non-unit, upper or lower triangular matrix, supplied in packed form. 00057 No test for singularity or near-singularity is included in this 00058 routine. Such tests must be performed before calling this routine. 00059 Parameters 00060 ========== 00061 UPLO - CHARACTER*1. 00062 On entry, UPLO specifies whether the matrix is an upper or 00063 lower triangular matrix as follows: 00064 UPLO = 'U' or 'u' A is an upper triangular matrix. 00065 UPLO = 'L' or 'l' A is a lower triangular matrix. 00066 Unchanged on exit. 00067 TRANS - CHARACTER*1. 00068 On entry, TRANS specifies the equations to be solved as 00069 follows: 00070 TRANS = 'N' or 'n' A*x = b. 00071 TRANS = 'T' or 't' A'*x = b. 00072 TRANS = 'C' or 'c' A'*x = b. 00073 Unchanged on exit. 00074 DIAG - CHARACTER*1. 00075 On entry, DIAG specifies whether or not A is unit 00076 triangular as follows: 00077 DIAG = 'U' or 'u' A is assumed to be unit triangular. 00078 DIAG = 'N' or 'n' A is not assumed to be unit 00079 triangular. 00080 Unchanged on exit. 00081 N - INTEGER. 00082 On entry, N specifies the order of the matrix A. 00083 N must be at least zero. 00084 Unchanged on exit. 00085 AP - DOUBLE PRECISION array of DIMENSION at least 00086 ( ( n*( n + 1 ) )/2 ). 00087 Before entry with UPLO = 'U' or 'u', the array AP must 00088 contain the upper triangular matrix packed sequentially, 00089 column by column, so that AP( 1 ) contains a( 1, 1 ), 00090 AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) 00091 respectively, and so on. 00092 Before entry with UPLO = 'L' or 'l', the array AP must 00093 contain the lower triangular matrix packed sequentially, 00094 column by column, so that AP( 1 ) contains a( 1, 1 ), 00095 AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) 00096 respectively, and so on. 00097 Note that when DIAG = 'U' or 'u', the diagonal elements of 00098 A are not referenced, but are assumed to be unity. 00099 Unchanged on exit. 00100 X - DOUBLE PRECISION array of dimension at least 00101 ( 1 + ( n - 1 )*abs( INCX ) ). 00102 Before entry, the incremented array X must contain the n 00103 element right-hand side vector b. On exit, X is overwritten 00104 with the solution vector x. 00105 INCX - INTEGER. 00106 On entry, INCX specifies the increment for the elements of 00107 X. INCX must not be zero. 00108 Unchanged on exit. 00109 Level 2 Blas routine. 00110 -- Written on 22-October-1986. 00111 Jack Dongarra, Argonne National Lab. 00112 Jeremy Du Croz, Nag Central Office. 00113 Sven Hammarling, Nag Central Office. 00114 Richard Hanson, Sandia National Labs. 00115 Test the input parameters. 00116 Parameter adjustments */ 00117 --x; 00118 --ap; 00119 /* Initialization added by Elias to get rid of compiler warnings. */ 00120 kx = 0; 00121 /* Function Body */ 00122 info = 0; 00123 if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) { 00124 info = 1; 00125 } else if (! template_blas_lsame(trans, "N") && ! template_blas_lsame(trans, 00126 "T") && ! template_blas_lsame(trans, "C")) { 00127 info = 2; 00128 } else if (! template_blas_lsame(diag, "U") && ! template_blas_lsame(diag, 00129 "N")) { 00130 info = 3; 00131 } else if (*n < 0) { 00132 info = 4; 00133 } else if (*incx == 0) { 00134 info = 7; 00135 } 00136 if (info != 0) { 00137 template_blas_erbla("DTPSV ", &info); 00138 return 0; 00139 } 00140 /* Quick return if possible. */ 00141 if (*n == 0) { 00142 return 0; 00143 } 00144 nounit = template_blas_lsame(diag, "N"); 00145 /* Set up the start point in X if the increment is not unity. This 00146 will be ( N - 1 )*INCX too small for descending loops. */ 00147 if (*incx <= 0) { 00148 kx = 1 - (*n - 1) * *incx; 00149 } else if (*incx != 1) { 00150 kx = 1; 00151 } 00152 /* Start the operations. In this version the elements of AP are 00153 accessed sequentially with one pass through AP. */ 00154 if (template_blas_lsame(trans, "N")) { 00155 /* Form x := inv( A )*x. */ 00156 if (template_blas_lsame(uplo, "U")) { 00157 kk = *n * (*n + 1) / 2; 00158 if (*incx == 1) { 00159 for (j = *n; j >= 1; --j) { 00160 if (x[j] != 0.) { 00161 if (nounit) { 00162 x[j] /= ap[kk]; 00163 } 00164 temp = x[j]; 00165 k = kk - 1; 00166 for (i__ = j - 1; i__ >= 1; --i__) { 00167 x[i__] -= temp * ap[k]; 00168 --k; 00169 /* L10: */ 00170 } 00171 } 00172 kk -= j; 00173 /* L20: */ 00174 } 00175 } else { 00176 jx = kx + (*n - 1) * *incx; 00177 for (j = *n; j >= 1; --j) { 00178 if (x[jx] != 0.) { 00179 if (nounit) { 00180 x[jx] /= ap[kk]; 00181 } 00182 temp = x[jx]; 00183 ix = jx; 00184 i__1 = kk - j + 1; 00185 for (k = kk - 1; k >= i__1; --k) { 00186 ix -= *incx; 00187 x[ix] -= temp * ap[k]; 00188 /* L30: */ 00189 } 00190 } 00191 jx -= *incx; 00192 kk -= j; 00193 /* L40: */ 00194 } 00195 } 00196 } else { 00197 kk = 1; 00198 if (*incx == 1) { 00199 i__1 = *n; 00200 for (j = 1; j <= i__1; ++j) { 00201 if (x[j] != 0.) { 00202 if (nounit) { 00203 x[j] /= ap[kk]; 00204 } 00205 temp = x[j]; 00206 k = kk + 1; 00207 i__2 = *n; 00208 for (i__ = j + 1; i__ <= i__2; ++i__) { 00209 x[i__] -= temp * ap[k]; 00210 ++k; 00211 /* L50: */ 00212 } 00213 } 00214 kk += *n - j + 1; 00215 /* L60: */ 00216 } 00217 } else { 00218 jx = kx; 00219 i__1 = *n; 00220 for (j = 1; j <= i__1; ++j) { 00221 if (x[jx] != 0.) { 00222 if (nounit) { 00223 x[jx] /= ap[kk]; 00224 } 00225 temp = x[jx]; 00226 ix = jx; 00227 i__2 = kk + *n - j; 00228 for (k = kk + 1; k <= i__2; ++k) { 00229 ix += *incx; 00230 x[ix] -= temp * ap[k]; 00231 /* L70: */ 00232 } 00233 } 00234 jx += *incx; 00235 kk += *n - j + 1; 00236 /* L80: */ 00237 } 00238 } 00239 } 00240 } else { 00241 /* Form x := inv( A' )*x. */ 00242 if (template_blas_lsame(uplo, "U")) { 00243 kk = 1; 00244 if (*incx == 1) { 00245 i__1 = *n; 00246 for (j = 1; j <= i__1; ++j) { 00247 temp = x[j]; 00248 k = kk; 00249 i__2 = j - 1; 00250 for (i__ = 1; i__ <= i__2; ++i__) { 00251 temp -= ap[k] * x[i__]; 00252 ++k; 00253 /* L90: */ 00254 } 00255 if (nounit) { 00256 temp /= ap[kk + j - 1]; 00257 } 00258 x[j] = temp; 00259 kk += j; 00260 /* L100: */ 00261 } 00262 } else { 00263 jx = kx; 00264 i__1 = *n; 00265 for (j = 1; j <= i__1; ++j) { 00266 temp = x[jx]; 00267 ix = kx; 00268 i__2 = kk + j - 2; 00269 for (k = kk; k <= i__2; ++k) { 00270 temp -= ap[k] * x[ix]; 00271 ix += *incx; 00272 /* L110: */ 00273 } 00274 if (nounit) { 00275 temp /= ap[kk + j - 1]; 00276 } 00277 x[jx] = temp; 00278 jx += *incx; 00279 kk += j; 00280 /* L120: */ 00281 } 00282 } 00283 } else { 00284 kk = *n * (*n + 1) / 2; 00285 if (*incx == 1) { 00286 for (j = *n; j >= 1; --j) { 00287 temp = x[j]; 00288 k = kk; 00289 i__1 = j + 1; 00290 for (i__ = *n; i__ >= i__1; --i__) { 00291 temp -= ap[k] * x[i__]; 00292 --k; 00293 /* L130: */ 00294 } 00295 if (nounit) { 00296 temp /= ap[kk - *n + j]; 00297 } 00298 x[j] = temp; 00299 kk -= *n - j + 1; 00300 /* L140: */ 00301 } 00302 } else { 00303 kx += (*n - 1) * *incx; 00304 jx = kx; 00305 for (j = *n; j >= 1; --j) { 00306 temp = x[jx]; 00307 ix = kx; 00308 i__1 = kk - (*n - (j + 1)); 00309 for (k = kk; k >= i__1; --k) { 00310 temp -= ap[k] * x[ix]; 00311 ix -= *incx; 00312 /* L150: */ 00313 } 00314 if (nounit) { 00315 temp /= ap[kk - *n + j]; 00316 } 00317 x[jx] = temp; 00318 jx -= *incx; 00319 kk -= *n - j + 1; 00320 /* L160: */ 00321 } 00322 } 00323 } 00324 } 00325 return 0; 00326 /* End of DTPSV . */ 00327 } /* dtpsv_ */ 00328 00329 #endif