template_blas_spr.h

Go to the documentation of this file.
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

Generated on Mon Sep 17 14:30:40 2012 for ergo by  doxygen 1.4.7