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