template_blas_spr2.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_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

Generated on 21 Nov 2012 for ergo by  doxygen 1.6.1