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_SPR2_HEADER 00038 #define TEMPLATE_BLAS_SPR2_HEADER 00039 00040 #include "template_blas_common.h" 00041 00042 template<class Treal> 00043 int template_blas_spr2(const char *uplo, const integer *n, const Treal *alpha, 00044 const Treal *x, const integer *incx, const Treal *y, const integer *incy, 00045 Treal *ap) 00046 { 00047 /* System generated locals */ 00048 integer i__1, i__2; 00049 /* Local variables */ 00050 integer info; 00051 Treal temp1, temp2; 00052 integer i__, j, k; 00053 integer kk, ix, iy, jx, jy, kx, ky; 00054 /* Purpose 00055 ======= 00056 DSPR2 performs the symmetric rank 2 operation 00057 A := alpha*x*y' + alpha*y*x' + A, 00058 where alpha is a scalar, x and y are n element vectors and A is an 00059 n by n symmetric matrix, supplied in packed form. 00060 Parameters 00061 ========== 00062 UPLO - CHARACTER*1. 00063 On entry, UPLO specifies whether the upper or lower 00064 triangular part of the matrix A is supplied in the packed 00065 array AP as follows: 00066 UPLO = 'U' or 'u' The upper triangular part of A is 00067 supplied in AP. 00068 UPLO = 'L' or 'l' The lower triangular part of A is 00069 supplied in AP. 00070 Unchanged on exit. 00071 N - INTEGER. 00072 On entry, N specifies the order of the matrix A. 00073 N must be at least zero. 00074 Unchanged on exit. 00075 ALPHA - DOUBLE PRECISION. 00076 On entry, ALPHA specifies the scalar alpha. 00077 Unchanged on exit. 00078 X - DOUBLE PRECISION array of dimension at least 00079 ( 1 + ( n - 1 )*abs( INCX ) ). 00080 Before entry, the incremented array X must contain the n 00081 element vector x. 00082 Unchanged on exit. 00083 INCX - INTEGER. 00084 On entry, INCX specifies the increment for the elements of 00085 X. INCX must not be zero. 00086 Unchanged on exit. 00087 Y - DOUBLE PRECISION array of dimension at least 00088 ( 1 + ( n - 1 )*abs( INCY ) ). 00089 Before entry, the incremented array Y must contain the n 00090 element vector y. 00091 Unchanged on exit. 00092 INCY - INTEGER. 00093 On entry, INCY specifies the increment for the elements of 00094 Y. INCY must not be zero. 00095 Unchanged on exit. 00096 AP - DOUBLE PRECISION array of DIMENSION at least 00097 ( ( n*( n + 1 ) )/2 ). 00098 Before entry with UPLO = 'U' or 'u', the array AP must 00099 contain the upper triangular part of the symmetric matrix 00100 packed sequentially, column by column, so that AP( 1 ) 00101 contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) 00102 and a( 2, 2 ) respectively, and so on. On exit, the array 00103 AP is overwritten by the upper triangular part of the 00104 updated matrix. 00105 Before entry with UPLO = 'L' or 'l', the array AP must 00106 contain the lower triangular part of the symmetric matrix 00107 packed sequentially, column by column, so that AP( 1 ) 00108 contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) 00109 and a( 3, 1 ) respectively, and so on. On exit, the array 00110 AP is overwritten by the lower triangular part of the 00111 updated matrix. 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 --ap; 00121 --y; 00122 --x; 00123 /* Initialization added by Elias to get rid of compiler warnings. */ 00124 jx = jy = kx = ky = 0; 00125 /* Function Body */ 00126 info = 0; 00127 if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) { 00128 info = 1; 00129 } else if (*n < 0) { 00130 info = 2; 00131 } else if (*incx == 0) { 00132 info = 5; 00133 } else if (*incy == 0) { 00134 info = 7; 00135 } 00136 if (info != 0) { 00137 template_blas_erbla("SPR2 ", &info); 00138 return 0; 00139 } 00140 /* Quick return if possible. */ 00141 if (*n == 0 || *alpha == 0.) { 00142 return 0; 00143 } 00144 /* Set up the start points in X and Y if the increments are not both 00145 unity. */ 00146 if (*incx != 1 || *incy != 1) { 00147 if (*incx > 0) { 00148 kx = 1; 00149 } else { 00150 kx = 1 - (*n - 1) * *incx; 00151 } 00152 if (*incy > 0) { 00153 ky = 1; 00154 } else { 00155 ky = 1 - (*n - 1) * *incy; 00156 } 00157 jx = kx; 00158 jy = ky; 00159 } 00160 /* Start the operations. In this version the elements of the array AP 00161 are accessed sequentially with one pass through AP. */ 00162 kk = 1; 00163 if (template_blas_lsame(uplo, "U")) { 00164 /* Form A when upper triangle is stored in AP. */ 00165 if (*incx == 1 && *incy == 1) { 00166 i__1 = *n; 00167 for (j = 1; j <= i__1; ++j) { 00168 if (x[j] != 0. || y[j] != 0.) { 00169 temp1 = *alpha * y[j]; 00170 temp2 = *alpha * x[j]; 00171 k = kk; 00172 i__2 = j; 00173 for (i__ = 1; i__ <= i__2; ++i__) { 00174 ap[k] = ap[k] + x[i__] * temp1 + y[i__] * temp2; 00175 ++k; 00176 /* L10: */ 00177 } 00178 } 00179 kk += j; 00180 /* L20: */ 00181 } 00182 } else { 00183 i__1 = *n; 00184 for (j = 1; j <= i__1; ++j) { 00185 if (x[jx] != 0. || y[jy] != 0.) { 00186 temp1 = *alpha * y[jy]; 00187 temp2 = *alpha * x[jx]; 00188 ix = kx; 00189 iy = ky; 00190 i__2 = kk + j - 1; 00191 for (k = kk; k <= i__2; ++k) { 00192 ap[k] = ap[k] + x[ix] * temp1 + y[iy] * temp2; 00193 ix += *incx; 00194 iy += *incy; 00195 /* L30: */ 00196 } 00197 } 00198 jx += *incx; 00199 jy += *incy; 00200 kk += j; 00201 /* L40: */ 00202 } 00203 } 00204 } else { 00205 /* Form A when lower triangle is stored in AP. */ 00206 if (*incx == 1 && *incy == 1) { 00207 i__1 = *n; 00208 for (j = 1; j <= i__1; ++j) { 00209 if (x[j] != 0. || y[j] != 0.) { 00210 temp1 = *alpha * y[j]; 00211 temp2 = *alpha * x[j]; 00212 k = kk; 00213 i__2 = *n; 00214 for (i__ = j; i__ <= i__2; ++i__) { 00215 ap[k] = ap[k] + x[i__] * temp1 + y[i__] * temp2; 00216 ++k; 00217 /* L50: */ 00218 } 00219 } 00220 kk = kk + *n - j + 1; 00221 /* L60: */ 00222 } 00223 } else { 00224 i__1 = *n; 00225 for (j = 1; j <= i__1; ++j) { 00226 if (x[jx] != 0. || y[jy] != 0.) { 00227 temp1 = *alpha * y[jy]; 00228 temp2 = *alpha * x[jx]; 00229 ix = jx; 00230 iy = jy; 00231 i__2 = kk + *n - j; 00232 for (k = kk; k <= i__2; ++k) { 00233 ap[k] = ap[k] + x[ix] * temp1 + y[iy] * temp2; 00234 ix += *incx; 00235 iy += *incy; 00236 /* L70: */ 00237 } 00238 } 00239 jx += *incx; 00240 jy += *incy; 00241 kk = kk + *n - j + 1; 00242 /* L80: */ 00243 } 00244 } 00245 } 00246 return 0; 00247 /* End of DSPR2 . */ 00248 } /* dspr2_ */ 00249 00250 #endif