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_SPR_HEADER 00036 #define TEMPLATE_BLAS_SPR_HEADER 00037 00038 #include "template_blas_common.h" 00039 00040 template<class Treal> 00041 int template_blas_spr(const char *uplo, integer *n, Treal *alpha, 00042 Treal *x, integer *incx, Treal *ap) 00043 { 00044 /* System generated locals */ 00045 integer i__1, i__2; 00046 /* Local variables */ 00047 integer info; 00048 Treal temp; 00049 integer i__, j, k; 00050 integer kk, ix, jx, kx; 00051 /* Purpose 00052 ======= 00053 DSPR performs the symmetric rank 1 operation 00054 A := alpha*x*x' + A, 00055 where alpha is a real scalar, x is an n element vector and A is an 00056 n by n symmetric matrix, supplied in packed form. 00057 Parameters 00058 ========== 00059 UPLO - CHARACTER*1. 00060 On entry, UPLO specifies whether the upper or lower 00061 triangular part of the matrix A is supplied in the packed 00062 array AP as follows: 00063 UPLO = 'U' or 'u' The upper triangular part of A is 00064 supplied in AP. 00065 UPLO = 'L' or 'l' The lower triangular part of A is 00066 supplied in AP. 00067 Unchanged on exit. 00068 N - INTEGER. 00069 On entry, N specifies the order of the matrix A. 00070 N must be at least zero. 00071 Unchanged on exit. 00072 ALPHA - DOUBLE PRECISION. 00073 On entry, ALPHA specifies the scalar alpha. 00074 Unchanged on exit. 00075 X - DOUBLE PRECISION array of dimension at least 00076 ( 1 + ( n - 1 )*abs( INCX ) ). 00077 Before entry, the incremented array X must contain the n 00078 element vector x. 00079 Unchanged on exit. 00080 INCX - INTEGER. 00081 On entry, INCX specifies the increment for the elements of 00082 X. INCX must not be zero. 00083 Unchanged on exit. 00084 AP - DOUBLE PRECISION array of DIMENSION at least 00085 ( ( n*( n + 1 ) )/2 ). 00086 Before entry with UPLO = 'U' or 'u', the array AP must 00087 contain the upper triangular part of the symmetric matrix 00088 packed sequentially, column by column, so that AP( 1 ) 00089 contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) 00090 and a( 2, 2 ) respectively, and so on. On exit, the array 00091 AP is overwritten by the upper triangular part of the 00092 updated matrix. 00093 Before entry with UPLO = 'L' or 'l', the array AP must 00094 contain the lower triangular part of the symmetric matrix 00095 packed sequentially, column by column, so that AP( 1 ) 00096 contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) 00097 and a( 3, 1 ) respectively, and so on. On exit, the array 00098 AP is overwritten by the lower triangular part of the 00099 updated matrix. 00100 Level 2 Blas routine. 00101 -- Written on 22-October-1986. 00102 Jack Dongarra, Argonne National Lab. 00103 Jeremy Du Croz, Nag Central Office. 00104 Sven Hammarling, Nag Central Office. 00105 Richard Hanson, Sandia National Labs. 00106 Test the input parameters. 00107 Parameter adjustments */ 00108 --ap; 00109 --x; 00110 /* Initialization added by Elias to get rid of compiler warnings. */ 00111 kx = 0; 00112 /* Function Body */ 00113 info = 0; 00114 if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) { 00115 info = 1; 00116 } else if (*n < 0) { 00117 info = 2; 00118 } else if (*incx == 0) { 00119 info = 5; 00120 } 00121 if (info != 0) { 00122 template_blas_erbla("SPR ", &info); 00123 return 0; 00124 } 00125 /* Quick return if possible. */ 00126 if (*n == 0 || *alpha == 0.) { 00127 return 0; 00128 } 00129 /* Set the start point in X if the increment is not unity. */ 00130 if (*incx <= 0) { 00131 kx = 1 - (*n - 1) * *incx; 00132 } else if (*incx != 1) { 00133 kx = 1; 00134 } 00135 /* Start the operations. In this version the elements of the array AP 00136 are accessed sequentially with one pass through AP. */ 00137 kk = 1; 00138 if (template_blas_lsame(uplo, "U")) { 00139 /* Form A when upper triangle is stored in AP. */ 00140 if (*incx == 1) { 00141 i__1 = *n; 00142 for (j = 1; j <= i__1; ++j) { 00143 if (x[j] != 0.) { 00144 temp = *alpha * x[j]; 00145 k = kk; 00146 i__2 = j; 00147 for (i__ = 1; i__ <= i__2; ++i__) { 00148 ap[k] += x[i__] * temp; 00149 ++k; 00150 /* L10: */ 00151 } 00152 } 00153 kk += j; 00154 /* L20: */ 00155 } 00156 } else { 00157 jx = kx; 00158 i__1 = *n; 00159 for (j = 1; j <= i__1; ++j) { 00160 if (x[jx] != 0.) { 00161 temp = *alpha * x[jx]; 00162 ix = kx; 00163 i__2 = kk + j - 1; 00164 for (k = kk; k <= i__2; ++k) { 00165 ap[k] += x[ix] * temp; 00166 ix += *incx; 00167 /* L30: */ 00168 } 00169 } 00170 jx += *incx; 00171 kk += j; 00172 /* L40: */ 00173 } 00174 } 00175 } else { 00176 /* Form A when lower triangle is stored in AP. */ 00177 if (*incx == 1) { 00178 i__1 = *n; 00179 for (j = 1; j <= i__1; ++j) { 00180 if (x[j] != 0.) { 00181 temp = *alpha * x[j]; 00182 k = kk; 00183 i__2 = *n; 00184 for (i__ = j; i__ <= i__2; ++i__) { 00185 ap[k] += x[i__] * temp; 00186 ++k; 00187 /* L50: */ 00188 } 00189 } 00190 kk = kk + *n - j + 1; 00191 /* L60: */ 00192 } 00193 } else { 00194 jx = kx; 00195 i__1 = *n; 00196 for (j = 1; j <= i__1; ++j) { 00197 if (x[jx] != 0.) { 00198 temp = *alpha * x[jx]; 00199 ix = jx; 00200 i__2 = kk + *n - j; 00201 for (k = kk; k <= i__2; ++k) { 00202 ap[k] += x[ix] * temp; 00203 ix += *incx; 00204 /* L70: */ 00205 } 00206 } 00207 jx += *incx; 00208 kk = kk + *n - j + 1; 00209 /* L80: */ 00210 } 00211 } 00212 } 00213 return 0; 00214 /* End of DSPR . */ 00215 } /* dspr_ */ 00216 00217 #endif