template_blas_syr2.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_SYR2_HEADER
00036 #define TEMPLATE_BLAS_SYR2_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_blas_syr2(const char *uplo, const integer *n, const Treal *alpha, 
00041         const Treal *x, const integer *incx, const Treal *y, const integer *incy, 
00042         Treal *a, const integer *lda)
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     DSYR2  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 n   
00057     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     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     A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).   
00095              Before entry with  UPLO = 'U' or 'u', the leading n by n   
00096              upper triangular part of the array A must contain the upper   
00097              triangular part of the symmetric matrix and the strictly   
00098              lower triangular part of A is not referenced. On exit, the   
00099              upper triangular part of the array A is overwritten by the   
00100              upper triangular part of the updated matrix.   
00101              Before entry with UPLO = 'L' or 'l', the leading n by n   
00102              lower triangular part of the array A must contain the lower   
00103              triangular part of the symmetric matrix and the strictly   
00104              upper triangular part of A is not referenced. On exit, the   
00105              lower triangular part of the array A is overwritten by the   
00106              lower triangular part of the updated matrix.   
00107     LDA    - INTEGER.   
00108              On entry, LDA specifies the first dimension of A as declared   
00109              in the calling (sub) program. LDA must be at least   
00110              max( 1, n ).   
00111              Unchanged on exit.   
00112     Level 2 Blas routine.   
00113     -- Written on 22-October-1986.   
00114        Jack Dongarra, Argonne National Lab.   
00115        Jeremy Du Croz, Nag Central Office.   
00116        Sven Hammarling, Nag Central Office.   
00117        Richard Hanson, Sandia National Labs.   
00118        Test the input parameters.   
00119        Parameter adjustments */
00120     --x;
00121     --y;
00122     a_dim1 = *lda;
00123     a_offset = 1 + a_dim1 * 1;
00124     a -= a_offset;
00125     /* Initialization added by Elias to get rid of compiler warnings. */
00126     jx = jy = kx = ky = 0;
00127     /* Function Body */
00128     info = 0;
00129     if (! template_blas_lsame(uplo, "U") && ! template_blas_lsame(uplo, "L")) {
00130         info = 1;
00131     } else if (*n < 0) {
00132         info = 2;
00133     } else if (*incx == 0) {
00134         info = 5;
00135     } else if (*incy == 0) {
00136         info = 7;
00137     } else if (*lda < maxMACRO(1,*n)) {
00138         info = 9;
00139     }
00140     if (info != 0) {
00141         template_blas_erbla("SYR2  ", &info);
00142         return 0;
00143     }
00144 /*     Quick return if possible. */
00145     if (*n == 0 || *alpha == 0.) {
00146         return 0;
00147     }
00148 /*     Set up the start points in X and Y if the increments are not both   
00149        unity. */
00150     if (*incx != 1 || *incy != 1) {
00151         if (*incx > 0) {
00152             kx = 1;
00153         } else {
00154             kx = 1 - (*n - 1) * *incx;
00155         }
00156         if (*incy > 0) {
00157             ky = 1;
00158         } else {
00159             ky = 1 - (*n - 1) * *incy;
00160         }
00161         jx = kx;
00162         jy = ky;
00163     }
00164 /*     Start the operations. In this version the elements of A are   
00165        accessed sequentially with one pass through the triangular part   
00166        of A. */
00167     if (template_blas_lsame(uplo, "U")) {
00168 /*        Form  A  when A is stored in the upper triangle. */
00169         if (*incx == 1 && *incy == 1) {
00170             i__1 = *n;
00171             for (j = 1; j <= i__1; ++j) {
00172                 if (x[j] != 0. || y[j] != 0.) {
00173                     temp1 = *alpha * y[j];
00174                     temp2 = *alpha * x[j];
00175                     i__2 = j;
00176                     for (i__ = 1; i__ <= i__2; ++i__) {
00177                         a_ref(i__, j) = a_ref(i__, j) + x[i__] * temp1 + y[
00178                                 i__] * temp2;
00179 /* L10: */
00180                     }
00181                 }
00182 /* L20: */
00183             }
00184         } else {
00185             i__1 = *n;
00186             for (j = 1; j <= i__1; ++j) {
00187                 if (x[jx] != 0. || y[jy] != 0.) {
00188                     temp1 = *alpha * y[jy];
00189                     temp2 = *alpha * x[jx];
00190                     ix = kx;
00191                     iy = ky;
00192                     i__2 = j;
00193                     for (i__ = 1; i__ <= i__2; ++i__) {
00194                         a_ref(i__, j) = a_ref(i__, j) + x[ix] * temp1 + y[iy] 
00195                                 * temp2;
00196                         ix += *incx;
00197                         iy += *incy;
00198 /* L30: */
00199                     }
00200                 }
00201                 jx += *incx;
00202                 jy += *incy;
00203 /* L40: */
00204             }
00205         }
00206     } else {
00207 /*        Form  A  when A is stored in the lower triangle. */
00208         if (*incx == 1 && *incy == 1) {
00209             i__1 = *n;
00210             for (j = 1; j <= i__1; ++j) {
00211                 if (x[j] != 0. || y[j] != 0.) {
00212                     temp1 = *alpha * y[j];
00213                     temp2 = *alpha * x[j];
00214                     i__2 = *n;
00215                     for (i__ = j; i__ <= i__2; ++i__) {
00216                         a_ref(i__, j) = a_ref(i__, j) + x[i__] * temp1 + y[
00217                                 i__] * temp2;
00218 /* L50: */
00219                     }
00220                 }
00221 /* L60: */
00222             }
00223         } else {
00224             i__1 = *n;
00225             for (j = 1; j <= i__1; ++j) {
00226                 if (x[jx] != 0. || y[jy] != 0.) {
00227                     temp1 = *alpha * y[jy];
00228                     temp2 = *alpha * x[jx];
00229                     ix = jx;
00230                     iy = jy;
00231                     i__2 = *n;
00232                     for (i__ = j; i__ <= i__2; ++i__) {
00233                         a_ref(i__, j) = a_ref(i__, j) + x[ix] * temp1 + y[iy] 
00234                                 * temp2;
00235                         ix += *incx;
00236                         iy += *incy;
00237 /* L70: */
00238                     }
00239                 }
00240                 jx += *incx;
00241                 jy += *incy;
00242 /* L80: */
00243             }
00244         }
00245     }
00246     return 0;
00247 /*     End of DSYR2 . */
00248 } /* dsyr2_ */
00249 #undef a_ref
00250 
00251 #endif

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