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