template_blas_symv.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_SYMV_HEADER
00036 #define TEMPLATE_BLAS_SYMV_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_blas_symv(const char *uplo, const integer *n, const Treal *alpha, 
00041         const Treal *a, const integer *lda, const Treal *x, const integer *incx, const Treal 
00042         *beta, Treal *y, const integer *incy)
00043 {
00044     /* System generated locals */
00045     integer a_dim1, a_offset, i__1, i__2;
00046     /* Local variables */
00047      integer info;
00048      Treal temp1, temp2;
00049      integer i__, j;
00050      integer ix, iy, jx, jy, kx, ky;
00051 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00052 /*  Purpose   
00053     =======   
00054     DSYMV  performs the matrix-vector  operation   
00055        y := alpha*A*x + beta*y,   
00056     where alpha and beta are scalars, x and y are n element vectors and   
00057     A is an n by n symmetric matrix.   
00058     Parameters   
00059     ==========   
00060     UPLO   - CHARACTER*1.   
00061              On entry, UPLO specifies whether the upper or lower   
00062              triangular part of the array A is to be referenced as   
00063              follows:   
00064                 UPLO = 'U' or 'u'   Only the upper triangular part of A   
00065                                     is to be referenced.   
00066                 UPLO = 'L' or 'l'   Only the lower triangular part of A   
00067                                     is to be referenced.   
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     A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).   
00077              Before entry with  UPLO = 'U' or 'u', the leading n by n   
00078              upper triangular part of the array A must contain the upper   
00079              triangular part of the symmetric matrix and the strictly   
00080              lower triangular part of A is not referenced.   
00081              Before entry with UPLO = 'L' or 'l', the leading n by n   
00082              lower triangular part of the array A must contain the lower   
00083              triangular part of the symmetric matrix and the strictly   
00084              upper triangular part of A is not referenced.   
00085              Unchanged on exit.   
00086     LDA    - INTEGER.   
00087              On entry, LDA specifies the first dimension of A as declared   
00088              in the calling (sub) program. LDA must be at least   
00089              max( 1, n ).   
00090              Unchanged on exit.   
00091     X      - DOUBLE PRECISION array of dimension at least   
00092              ( 1 + ( n - 1 )*abs( INCX ) ).   
00093              Before entry, the incremented array X must contain the n   
00094              element vector x.   
00095              Unchanged on exit.   
00096     INCX   - INTEGER.   
00097              On entry, INCX specifies the increment for the elements of   
00098              X. INCX must not be zero.   
00099              Unchanged on exit.   
00100     BETA   - DOUBLE PRECISION.   
00101              On entry, BETA specifies the scalar beta. When BETA is   
00102              supplied as zero then Y need not be set on input.   
00103              Unchanged on exit.   
00104     Y      - DOUBLE PRECISION array of dimension at least   
00105              ( 1 + ( n - 1 )*abs( INCY ) ).   
00106              Before entry, the incremented array Y must contain the n   
00107              element vector y. On exit, Y is overwritten by the updated   
00108              vector y.   
00109     INCY   - INTEGER.   
00110              On entry, INCY specifies the increment for the elements of   
00111              Y. INCY must not be zero.   
00112              Unchanged on exit.   
00113     Level 2 Blas routine.   
00114     -- Written on 22-October-1986.   
00115        Jack Dongarra, Argonne National Lab.   
00116        Jeremy Du Croz, Nag Central Office.   
00117        Sven Hammarling, Nag Central Office.   
00118        Richard Hanson, Sandia National Labs.   
00119        Test the input parameters.   
00120        Parameter adjustments */
00121     a_dim1 = *lda;
00122     a_offset = 1 + a_dim1 * 1;
00123     a -= a_offset;
00124     --x;
00125     --y;
00126     /* Function Body */
00127     info = 0;
00128     if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) {
00129         info = 1;
00130     } else if (*n < 0) {
00131         info = 2;
00132     } else if (*lda < maxMACRO(1,*n)) {
00133         info = 5;
00134     } else if (*incx == 0) {
00135         info = 7;
00136     } else if (*incy == 0) {
00137         info = 10;
00138     }
00139     if (info != 0) {
00140         template_blas_erbla("SYMV  ", &info);
00141         return 0;
00142     }
00143 /*     Quick return if possible. */
00144     if (*n == 0 || (*alpha == 0. && *beta == 1.) ) {
00145         return 0;
00146     }
00147 /*     Set up the start points in  X  and  Y. */
00148     if (*incx > 0) {
00149         kx = 1;
00150     } else {
00151         kx = 1 - (*n - 1) * *incx;
00152     }
00153     if (*incy > 0) {
00154         ky = 1;
00155     } else {
00156         ky = 1 - (*n - 1) * *incy;
00157     }
00158 /*     Start the operations. In this version the elements of A are   
00159        accessed sequentially with one pass through the triangular part   
00160        of A.   
00161        First form  y := beta*y. */
00162     if (*beta != 1.) {
00163         if (*incy == 1) {
00164             if (*beta == 0.) {
00165                 i__1 = *n;
00166                 for (i__ = 1; i__ <= i__1; ++i__) {
00167                     y[i__] = 0.;
00168 /* L10: */
00169                 }
00170             } else {
00171                 i__1 = *n;
00172                 for (i__ = 1; i__ <= i__1; ++i__) {
00173                     y[i__] = *beta * y[i__];
00174 /* L20: */
00175                 }
00176             }
00177         } else {
00178             iy = ky;
00179             if (*beta == 0.) {
00180                 i__1 = *n;
00181                 for (i__ = 1; i__ <= i__1; ++i__) {
00182                     y[iy] = 0.;
00183                     iy += *incy;
00184 /* L30: */
00185                 }
00186             } else {
00187                 i__1 = *n;
00188                 for (i__ = 1; i__ <= i__1; ++i__) {
00189                     y[iy] = *beta * y[iy];
00190                     iy += *incy;
00191 /* L40: */
00192                 }
00193             }
00194         }
00195     }
00196     if (*alpha == 0.) {
00197         return 0;
00198     }
00199     if (template_blas_lsame(uplo, "U")) {
00200 /*        Form  y  when A is stored in upper triangle. */
00201         if (*incx == 1 && *incy == 1) {
00202             i__1 = *n;
00203             for (j = 1; j <= i__1; ++j) {
00204                 temp1 = *alpha * x[j];
00205                 temp2 = 0.;
00206                 i__2 = j - 1;
00207                 for (i__ = 1; i__ <= i__2; ++i__) {
00208                     y[i__] += temp1 * a_ref(i__, j);
00209                     temp2 += a_ref(i__, j) * x[i__];
00210 /* L50: */
00211                 }
00212                 y[j] = y[j] + temp1 * a_ref(j, j) + *alpha * temp2;
00213 /* L60: */
00214             }
00215         } else {
00216             jx = kx;
00217             jy = ky;
00218             i__1 = *n;
00219             for (j = 1; j <= i__1; ++j) {
00220                 temp1 = *alpha * x[jx];
00221                 temp2 = 0.;
00222                 ix = kx;
00223                 iy = ky;
00224                 i__2 = j - 1;
00225                 for (i__ = 1; i__ <= i__2; ++i__) {
00226                     y[iy] += temp1 * a_ref(i__, j);
00227                     temp2 += a_ref(i__, j) * x[ix];
00228                     ix += *incx;
00229                     iy += *incy;
00230 /* L70: */
00231                 }
00232                 y[jy] = y[jy] + temp1 * a_ref(j, j) + *alpha * temp2;
00233                 jx += *incx;
00234                 jy += *incy;
00235 /* L80: */
00236             }
00237         }
00238     } else {
00239 /*        Form  y  when A is stored in lower triangle. */
00240         if (*incx == 1 && *incy == 1) {
00241             i__1 = *n;
00242             for (j = 1; j <= i__1; ++j) {
00243                 temp1 = *alpha * x[j];
00244                 temp2 = 0.;
00245                 y[j] += temp1 * a_ref(j, j);
00246                 i__2 = *n;
00247                 for (i__ = j + 1; i__ <= i__2; ++i__) {
00248                     y[i__] += temp1 * a_ref(i__, j);
00249                     temp2 += a_ref(i__, j) * x[i__];
00250 /* L90: */
00251                 }
00252                 y[j] += *alpha * temp2;
00253 /* L100: */
00254             }
00255         } else {
00256             jx = kx;
00257             jy = ky;
00258             i__1 = *n;
00259             for (j = 1; j <= i__1; ++j) {
00260                 temp1 = *alpha * x[jx];
00261                 temp2 = 0.;
00262                 y[jy] += temp1 * a_ref(j, j);
00263                 ix = jx;
00264                 iy = jy;
00265                 i__2 = *n;
00266                 for (i__ = j + 1; i__ <= i__2; ++i__) {
00267                     ix += *incx;
00268                     iy += *incy;
00269                     y[iy] += temp1 * a_ref(i__, j);
00270                     temp2 += a_ref(i__, j) * x[ix];
00271 /* L110: */
00272                 }
00273                 y[jy] += *alpha * temp2;
00274                 jx += *incx;
00275                 jy += *incy;
00276 /* L120: */
00277             }
00278         }
00279     }
00280     return 0;
00281 /*     End of DSYMV . */
00282 } /* dsymv_ */
00283 #undef a_ref
00284 
00285 #endif

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