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_SYR2_HEADER 00038 #define TEMPLATE_BLAS_SYR2_HEADER 00039 00040 00041 template<class Treal> 00042 int template_blas_syr2(const char *uplo, const integer *n, const Treal *alpha, 00043 const Treal *x, const integer *incx, const Treal *y, const integer *incy, 00044 Treal *a, const integer *lda) 00045 { 00046 /* System generated locals */ 00047 integer a_dim1, a_offset, i__1, i__2; 00048 /* Local variables */ 00049 integer info; 00050 Treal temp1, temp2; 00051 integer i__, j; 00052 integer ix, iy, jx, jy, kx, ky; 00053 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 00054 /* Purpose 00055 ======= 00056 DSYR2 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 n 00059 by n symmetric matrix. 00060 Parameters 00061 ========== 00062 UPLO - CHARACTER*1. 00063 On entry, UPLO specifies whether the upper or lower 00064 triangular part of the array A is to be referenced as 00065 follows: 00066 UPLO = 'U' or 'u' Only the upper triangular part of A 00067 is to be referenced. 00068 UPLO = 'L' or 'l' Only the lower triangular part of A 00069 is to be referenced. 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 A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). 00097 Before entry with UPLO = 'U' or 'u', the leading n by n 00098 upper triangular part of the array A must contain the upper 00099 triangular part of the symmetric matrix and the strictly 00100 lower triangular part of A is not referenced. On exit, the 00101 upper triangular part of the array A is overwritten by the 00102 upper triangular part of the updated matrix. 00103 Before entry with UPLO = 'L' or 'l', the leading n by n 00104 lower triangular part of the array A must contain the lower 00105 triangular part of the symmetric matrix and the strictly 00106 upper triangular part of A is not referenced. On exit, the 00107 lower triangular part of the array A is overwritten by the 00108 lower triangular part of the updated matrix. 00109 LDA - INTEGER. 00110 On entry, LDA specifies the first dimension of A as declared 00111 in the calling (sub) program. LDA must be at least 00112 max( 1, n ). 00113 Unchanged on exit. 00114 Level 2 Blas routine. 00115 -- Written on 22-October-1986. 00116 Jack Dongarra, Argonne National Lab. 00117 Jeremy Du Croz, Nag Central Office. 00118 Sven Hammarling, Nag Central Office. 00119 Richard Hanson, Sandia National Labs. 00120 Test the input parameters. 00121 Parameter adjustments */ 00122 --x; 00123 --y; 00124 a_dim1 = *lda; 00125 a_offset = 1 + a_dim1 * 1; 00126 a -= a_offset; 00127 /* Initialization added by Elias to get rid of compiler warnings. */ 00128 jx = jy = kx = ky = 0; 00129 /* Function Body */ 00130 info = 0; 00131 if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) { 00132 info = 1; 00133 } else if (*n < 0) { 00134 info = 2; 00135 } else if (*incx == 0) { 00136 info = 5; 00137 } else if (*incy == 0) { 00138 info = 7; 00139 } else if (*lda < maxMACRO(1,*n)) { 00140 info = 9; 00141 } 00142 if (info != 0) { 00143 template_blas_erbla("SYR2 ", &info); 00144 return 0; 00145 } 00146 /* Quick return if possible. */ 00147 if (*n == 0 || *alpha == 0.) { 00148 return 0; 00149 } 00150 /* Set up the start points in X and Y if the increments are not both 00151 unity. */ 00152 if (*incx != 1 || *incy != 1) { 00153 if (*incx > 0) { 00154 kx = 1; 00155 } else { 00156 kx = 1 - (*n - 1) * *incx; 00157 } 00158 if (*incy > 0) { 00159 ky = 1; 00160 } else { 00161 ky = 1 - (*n - 1) * *incy; 00162 } 00163 jx = kx; 00164 jy = ky; 00165 } 00166 /* Start the operations. In this version the elements of A are 00167 accessed sequentially with one pass through the triangular part 00168 of A. */ 00169 if (template_blas_lsame(uplo, "U")) { 00170 /* Form A when A is stored in the upper triangle. */ 00171 if (*incx == 1 && *incy == 1) { 00172 i__1 = *n; 00173 for (j = 1; j <= i__1; ++j) { 00174 if (x[j] != 0. || y[j] != 0.) { 00175 temp1 = *alpha * y[j]; 00176 temp2 = *alpha * x[j]; 00177 i__2 = j; 00178 for (i__ = 1; i__ <= i__2; ++i__) { 00179 a_ref(i__, j) = a_ref(i__, j) + x[i__] * temp1 + y[ 00180 i__] * temp2; 00181 /* L10: */ 00182 } 00183 } 00184 /* L20: */ 00185 } 00186 } else { 00187 i__1 = *n; 00188 for (j = 1; j <= i__1; ++j) { 00189 if (x[jx] != 0. || y[jy] != 0.) { 00190 temp1 = *alpha * y[jy]; 00191 temp2 = *alpha * x[jx]; 00192 ix = kx; 00193 iy = ky; 00194 i__2 = j; 00195 for (i__ = 1; i__ <= i__2; ++i__) { 00196 a_ref(i__, j) = a_ref(i__, j) + x[ix] * temp1 + y[iy] 00197 * temp2; 00198 ix += *incx; 00199 iy += *incy; 00200 /* L30: */ 00201 } 00202 } 00203 jx += *incx; 00204 jy += *incy; 00205 /* L40: */ 00206 } 00207 } 00208 } else { 00209 /* Form A when A is stored in the lower triangle. */ 00210 if (*incx == 1 && *incy == 1) { 00211 i__1 = *n; 00212 for (j = 1; j <= i__1; ++j) { 00213 if (x[j] != 0. || y[j] != 0.) { 00214 temp1 = *alpha * y[j]; 00215 temp2 = *alpha * x[j]; 00216 i__2 = *n; 00217 for (i__ = j; i__ <= i__2; ++i__) { 00218 a_ref(i__, j) = a_ref(i__, j) + x[i__] * temp1 + y[ 00219 i__] * temp2; 00220 /* L50: */ 00221 } 00222 } 00223 /* L60: */ 00224 } 00225 } else { 00226 i__1 = *n; 00227 for (j = 1; j <= i__1; ++j) { 00228 if (x[jx] != 0. || y[jy] != 0.) { 00229 temp1 = *alpha * y[jy]; 00230 temp2 = *alpha * x[jx]; 00231 ix = jx; 00232 iy = jy; 00233 i__2 = *n; 00234 for (i__ = j; i__ <= i__2; ++i__) { 00235 a_ref(i__, j) = a_ref(i__, j) + x[ix] * temp1 + y[iy] 00236 * temp2; 00237 ix += *incx; 00238 iy += *incy; 00239 /* L70: */ 00240 } 00241 } 00242 jx += *incx; 00243 jy += *incy; 00244 /* L80: */ 00245 } 00246 } 00247 } 00248 return 0; 00249 /* End of DSYR2 . */ 00250 } /* dsyr2_ */ 00251 #undef a_ref 00252 00253 #endif