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_SYR2_HEADER 00036 #define TEMPLATE_BLAS_SYR2_HEADER 00037 00038 00039 template<class Treal> 00040 int template_blas_syr2(const char *uplo, const integer *n, const Treal *alpha, 00041 const Treal *x, const integer *incx, const Treal *y, const integer *incy, 00042 Treal *a, const integer *lda) 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 DSYR2 performs the symmetric rank 2 operation 00055 A := alpha*x*y' + alpha*y*x' + A, 00056 where alpha is a scalar, x and y are n element vectors and A is an n 00057 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 X - DOUBLE PRECISION array of dimension at least 00077 ( 1 + ( n - 1 )*abs( INCX ) ). 00078 Before entry, the incremented array X must contain the n 00079 element vector x. 00080 Unchanged on exit. 00081 INCX - INTEGER. 00082 On entry, INCX specifies the increment for the elements of 00083 X. INCX must not be zero. 00084 Unchanged on exit. 00085 Y - DOUBLE PRECISION array of dimension at least 00086 ( 1 + ( n - 1 )*abs( INCY ) ). 00087 Before entry, the incremented array Y must contain the n 00088 element vector y. 00089 Unchanged on exit. 00090 INCY - INTEGER. 00091 On entry, INCY specifies the increment for the elements of 00092 Y. INCY must not be zero. 00093 Unchanged on exit. 00094 A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). 00095 Before entry with UPLO = 'U' or 'u', the leading n by n 00096 upper triangular part of the array A must contain the upper 00097 triangular part of the symmetric matrix and the strictly 00098 lower triangular part of A is not referenced. On exit, the 00099 upper triangular part of the array A is overwritten by the 00100 upper triangular part of the updated matrix. 00101 Before entry with UPLO = 'L' or 'l', the leading n by n 00102 lower triangular part of the array A must contain the lower 00103 triangular part of the symmetric matrix and the strictly 00104 upper triangular part of A is not referenced. On exit, the 00105 lower triangular part of the array A is overwritten by the 00106 lower triangular part of the updated matrix. 00107 LDA - INTEGER. 00108 On entry, LDA specifies the first dimension of A as declared 00109 in the calling (sub) program. LDA must be at least 00110 max( 1, n ). 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 --x; 00121 --y; 00122 a_dim1 = *lda; 00123 a_offset = 1 + a_dim1 * 1; 00124 a -= a_offset; 00125 /* Initialization added by Elias to get rid of compiler warnings. */ 00126 jx = jy = kx = ky = 0; 00127 /* Function Body */ 00128 info = 0; 00129 if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) { 00130 info = 1; 00131 } else if (*n < 0) { 00132 info = 2; 00133 } else if (*incx == 0) { 00134 info = 5; 00135 } else if (*incy == 0) { 00136 info = 7; 00137 } else if (*lda < maxMACRO(1,*n)) { 00138 info = 9; 00139 } 00140 if (info != 0) { 00141 template_blas_erbla("SYR2 ", &info); 00142 return 0; 00143 } 00144 /* Quick return if possible. */ 00145 if (*n == 0 || *alpha == 0.) { 00146 return 0; 00147 } 00148 /* Set up the start points in X and Y if the increments are not both 00149 unity. */ 00150 if (*incx != 1 || *incy != 1) { 00151 if (*incx > 0) { 00152 kx = 1; 00153 } else { 00154 kx = 1 - (*n - 1) * *incx; 00155 } 00156 if (*incy > 0) { 00157 ky = 1; 00158 } else { 00159 ky = 1 - (*n - 1) * *incy; 00160 } 00161 jx = kx; 00162 jy = ky; 00163 } 00164 /* Start the operations. In this version the elements of A are 00165 accessed sequentially with one pass through the triangular part 00166 of A. */ 00167 if (template_blas_lsame(uplo, "U")) { 00168 /* Form A when A is stored in the upper triangle. */ 00169 if (*incx == 1 && *incy == 1) { 00170 i__1 = *n; 00171 for (j = 1; j <= i__1; ++j) { 00172 if (x[j] != 0. || y[j] != 0.) { 00173 temp1 = *alpha * y[j]; 00174 temp2 = *alpha * x[j]; 00175 i__2 = j; 00176 for (i__ = 1; i__ <= i__2; ++i__) { 00177 a_ref(i__, j) = a_ref(i__, j) + x[i__] * temp1 + y[ 00178 i__] * temp2; 00179 /* L10: */ 00180 } 00181 } 00182 /* L20: */ 00183 } 00184 } else { 00185 i__1 = *n; 00186 for (j = 1; j <= i__1; ++j) { 00187 if (x[jx] != 0. || y[jy] != 0.) { 00188 temp1 = *alpha * y[jy]; 00189 temp2 = *alpha * x[jx]; 00190 ix = kx; 00191 iy = ky; 00192 i__2 = j; 00193 for (i__ = 1; i__ <= i__2; ++i__) { 00194 a_ref(i__, j) = a_ref(i__, j) + x[ix] * temp1 + y[iy] 00195 * temp2; 00196 ix += *incx; 00197 iy += *incy; 00198 /* L30: */ 00199 } 00200 } 00201 jx += *incx; 00202 jy += *incy; 00203 /* L40: */ 00204 } 00205 } 00206 } else { 00207 /* Form A when A is stored in the lower triangle. */ 00208 if (*incx == 1 && *incy == 1) { 00209 i__1 = *n; 00210 for (j = 1; j <= i__1; ++j) { 00211 if (x[j] != 0. || y[j] != 0.) { 00212 temp1 = *alpha * y[j]; 00213 temp2 = *alpha * x[j]; 00214 i__2 = *n; 00215 for (i__ = j; i__ <= i__2; ++i__) { 00216 a_ref(i__, j) = a_ref(i__, j) + x[i__] * temp1 + y[ 00217 i__] * temp2; 00218 /* L50: */ 00219 } 00220 } 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 = *n; 00232 for (i__ = j; i__ <= i__2; ++i__) { 00233 a_ref(i__, j) = a_ref(i__, j) + x[ix] * temp1 + y[iy] 00234 * temp2; 00235 ix += *incx; 00236 iy += *incy; 00237 /* L70: */ 00238 } 00239 } 00240 jx += *incx; 00241 jy += *incy; 00242 /* L80: */ 00243 } 00244 } 00245 } 00246 return 0; 00247 /* End of DSYR2 . */ 00248 } /* dsyr2_ */ 00249 #undef a_ref 00250 00251 #endif