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_SYMV_HEADER 00036 #define TEMPLATE_BLAS_SYMV_HEADER 00037 00038 00039 template<class Treal> 00040 int template_blas_symv(const char *uplo, const integer *n, const Treal *alpha, 00041 const Treal *a, const integer *lda, const Treal *x, const integer *incx, const Treal 00042 *beta, Treal *y, const integer *incy) 00043 { 00044 /* System generated locals */ 00045 integer a_dim1, a_offset, i__1, i__2; 00046 /* Local variables */ 00047 integer info; 00048 Treal temp1, temp2; 00049 integer i__, j; 00050 integer ix, iy, jx, jy, kx, ky; 00051 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 00052 /* Purpose 00053 ======= 00054 DSYMV performs the matrix-vector operation 00055 y := alpha*A*x + beta*y, 00056 where alpha and beta are scalars, x and y are n element vectors and 00057 A is an n by n symmetric matrix. 00058 Parameters 00059 ========== 00060 UPLO - CHARACTER*1. 00061 On entry, UPLO specifies whether the upper or lower 00062 triangular part of the array A is to be referenced as 00063 follows: 00064 UPLO = 'U' or 'u' Only the upper triangular part of A 00065 is to be referenced. 00066 UPLO = 'L' or 'l' Only the lower triangular part of A 00067 is to be referenced. 00068 Unchanged on exit. 00069 N - INTEGER. 00070 On entry, N specifies the order of the matrix A. 00071 N must be at least zero. 00072 Unchanged on exit. 00073 ALPHA - DOUBLE PRECISION. 00074 On entry, ALPHA specifies the scalar alpha. 00075 Unchanged on exit. 00076 A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). 00077 Before entry with UPLO = 'U' or 'u', the leading n by n 00078 upper triangular part of the array A must contain the upper 00079 triangular part of the symmetric matrix and the strictly 00080 lower triangular part of A is not referenced. 00081 Before entry with UPLO = 'L' or 'l', the leading n by n 00082 lower triangular part of the array A must contain the lower 00083 triangular part of the symmetric matrix and the strictly 00084 upper triangular part of A is not referenced. 00085 Unchanged on exit. 00086 LDA - INTEGER. 00087 On entry, LDA specifies the first dimension of A as declared 00088 in the calling (sub) program. LDA must be at least 00089 max( 1, n ). 00090 Unchanged on exit. 00091 X - DOUBLE PRECISION array of dimension at least 00092 ( 1 + ( n - 1 )*abs( INCX ) ). 00093 Before entry, the incremented array X must contain the n 00094 element vector x. 00095 Unchanged on exit. 00096 INCX - INTEGER. 00097 On entry, INCX specifies the increment for the elements of 00098 X. INCX must not be zero. 00099 Unchanged on exit. 00100 BETA - DOUBLE PRECISION. 00101 On entry, BETA specifies the scalar beta. When BETA is 00102 supplied as zero then Y need not be set on input. 00103 Unchanged on exit. 00104 Y - DOUBLE PRECISION array of dimension at least 00105 ( 1 + ( n - 1 )*abs( INCY ) ). 00106 Before entry, the incremented array Y must contain the n 00107 element vector y. On exit, Y is overwritten by the updated 00108 vector y. 00109 INCY - INTEGER. 00110 On entry, INCY specifies the increment for the elements of 00111 Y. INCY must not be zero. 00112 Unchanged on exit. 00113 Level 2 Blas routine. 00114 -- Written on 22-October-1986. 00115 Jack Dongarra, Argonne National Lab. 00116 Jeremy Du Croz, Nag Central Office. 00117 Sven Hammarling, Nag Central Office. 00118 Richard Hanson, Sandia National Labs. 00119 Test the input parameters. 00120 Parameter adjustments */ 00121 a_dim1 = *lda; 00122 a_offset = 1 + a_dim1 * 1; 00123 a -= a_offset; 00124 --x; 00125 --y; 00126 /* Function Body */ 00127 info = 0; 00128 if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) { 00129 info = 1; 00130 } else if (*n < 0) { 00131 info = 2; 00132 } else if (*lda < maxMACRO(1,*n)) { 00133 info = 5; 00134 } else if (*incx == 0) { 00135 info = 7; 00136 } else if (*incy == 0) { 00137 info = 10; 00138 } 00139 if (info != 0) { 00140 template_blas_erbla("SYMV ", &info); 00141 return 0; 00142 } 00143 /* Quick return if possible. */ 00144 if (*n == 0 || (*alpha == 0. && *beta == 1.) ) { 00145 return 0; 00146 } 00147 /* Set up the start points in X and Y. */ 00148 if (*incx > 0) { 00149 kx = 1; 00150 } else { 00151 kx = 1 - (*n - 1) * *incx; 00152 } 00153 if (*incy > 0) { 00154 ky = 1; 00155 } else { 00156 ky = 1 - (*n - 1) * *incy; 00157 } 00158 /* Start the operations. In this version the elements of A are 00159 accessed sequentially with one pass through the triangular part 00160 of A. 00161 First form y := beta*y. */ 00162 if (*beta != 1.) { 00163 if (*incy == 1) { 00164 if (*beta == 0.) { 00165 i__1 = *n; 00166 for (i__ = 1; i__ <= i__1; ++i__) { 00167 y[i__] = 0.; 00168 /* L10: */ 00169 } 00170 } else { 00171 i__1 = *n; 00172 for (i__ = 1; i__ <= i__1; ++i__) { 00173 y[i__] = *beta * y[i__]; 00174 /* L20: */ 00175 } 00176 } 00177 } else { 00178 iy = ky; 00179 if (*beta == 0.) { 00180 i__1 = *n; 00181 for (i__ = 1; i__ <= i__1; ++i__) { 00182 y[iy] = 0.; 00183 iy += *incy; 00184 /* L30: */ 00185 } 00186 } else { 00187 i__1 = *n; 00188 for (i__ = 1; i__ <= i__1; ++i__) { 00189 y[iy] = *beta * y[iy]; 00190 iy += *incy; 00191 /* L40: */ 00192 } 00193 } 00194 } 00195 } 00196 if (*alpha == 0.) { 00197 return 0; 00198 } 00199 if (template_blas_lsame(uplo, "U")) { 00200 /* Form y when A is stored in upper triangle. */ 00201 if (*incx == 1 && *incy == 1) { 00202 i__1 = *n; 00203 for (j = 1; j <= i__1; ++j) { 00204 temp1 = *alpha * x[j]; 00205 temp2 = 0.; 00206 i__2 = j - 1; 00207 for (i__ = 1; i__ <= i__2; ++i__) { 00208 y[i__] += temp1 * a_ref(i__, j); 00209 temp2 += a_ref(i__, j) * x[i__]; 00210 /* L50: */ 00211 } 00212 y[j] = y[j] + temp1 * a_ref(j, j) + *alpha * temp2; 00213 /* L60: */ 00214 } 00215 } else { 00216 jx = kx; 00217 jy = ky; 00218 i__1 = *n; 00219 for (j = 1; j <= i__1; ++j) { 00220 temp1 = *alpha * x[jx]; 00221 temp2 = 0.; 00222 ix = kx; 00223 iy = ky; 00224 i__2 = j - 1; 00225 for (i__ = 1; i__ <= i__2; ++i__) { 00226 y[iy] += temp1 * a_ref(i__, j); 00227 temp2 += a_ref(i__, j) * x[ix]; 00228 ix += *incx; 00229 iy += *incy; 00230 /* L70: */ 00231 } 00232 y[jy] = y[jy] + temp1 * a_ref(j, j) + *alpha * temp2; 00233 jx += *incx; 00234 jy += *incy; 00235 /* L80: */ 00236 } 00237 } 00238 } else { 00239 /* Form y when A is stored in lower triangle. */ 00240 if (*incx == 1 && *incy == 1) { 00241 i__1 = *n; 00242 for (j = 1; j <= i__1; ++j) { 00243 temp1 = *alpha * x[j]; 00244 temp2 = 0.; 00245 y[j] += temp1 * a_ref(j, j); 00246 i__2 = *n; 00247 for (i__ = j + 1; i__ <= i__2; ++i__) { 00248 y[i__] += temp1 * a_ref(i__, j); 00249 temp2 += a_ref(i__, j) * x[i__]; 00250 /* L90: */ 00251 } 00252 y[j] += *alpha * temp2; 00253 /* L100: */ 00254 } 00255 } else { 00256 jx = kx; 00257 jy = ky; 00258 i__1 = *n; 00259 for (j = 1; j <= i__1; ++j) { 00260 temp1 = *alpha * x[jx]; 00261 temp2 = 0.; 00262 y[jy] += temp1 * a_ref(j, j); 00263 ix = jx; 00264 iy = jy; 00265 i__2 = *n; 00266 for (i__ = j + 1; i__ <= i__2; ++i__) { 00267 ix += *incx; 00268 iy += *incy; 00269 y[iy] += temp1 * a_ref(i__, j); 00270 temp2 += a_ref(i__, j) * x[ix]; 00271 /* L110: */ 00272 } 00273 y[jy] += *alpha * temp2; 00274 jx += *incx; 00275 jy += *incy; 00276 /* L120: */ 00277 } 00278 } 00279 } 00280 return 0; 00281 /* End of DSYMV . */ 00282 } /* dsymv_ */ 00283 #undef a_ref 00284 00285 #endif