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