template_blas_syrk.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_SYRK_HEADER
00036 #define TEMPLATE_BLAS_SYRK_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_blas_syrk(const char *uplo, const char *trans, const integer *n, const integer *k, 
00041         const Treal *alpha, const Treal *a, const integer *lda, const Treal *beta, 
00042         Treal *c__, const integer *ldc)
00043 {
00044     /* System generated locals */
00045     integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2, i__3;
00046     /* Local variables */
00047      integer info;
00048      Treal temp;
00049      integer i__, j, l;
00050      integer nrowa;
00051      logical upper;
00052 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
00053 #define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1]
00054 /*  Purpose   
00055     =======   
00056     DSYRK  performs one of the symmetric rank k operations   
00057        C := alpha*A*A' + beta*C,   
00058     or   
00059        C := alpha*A'*A + beta*C,   
00060     where  alpha and beta  are scalars, C is an  n by n  symmetric matrix   
00061     and  A  is an  n by k  matrix in the first case and a  k by n  matrix   
00062     in the second case.   
00063     Parameters   
00064     ==========   
00065     UPLO   - CHARACTER*1.   
00066              On  entry,   UPLO  specifies  whether  the  upper  or  lower   
00067              triangular  part  of the  array  C  is to be  referenced  as   
00068              follows:   
00069                 UPLO = 'U' or 'u'   Only the  upper triangular part of  C   
00070                                     is to be referenced.   
00071                 UPLO = 'L' or 'l'   Only the  lower triangular part of  C   
00072                                     is to be referenced.   
00073              Unchanged on exit.   
00074     TRANS  - CHARACTER*1.   
00075              On entry,  TRANS  specifies the operation to be performed as   
00076              follows:   
00077                 TRANS = 'N' or 'n'   C := alpha*A*A' + beta*C.   
00078                 TRANS = 'T' or 't'   C := alpha*A'*A + beta*C.   
00079                 TRANS = 'C' or 'c'   C := alpha*A'*A + beta*C.   
00080              Unchanged on exit.   
00081     N      - INTEGER.   
00082              On entry,  N specifies the order of the matrix C.  N must be   
00083              at least zero.   
00084              Unchanged on exit.   
00085     K      - INTEGER.   
00086              On entry with  TRANS = 'N' or 'n',  K  specifies  the number   
00087              of  columns   of  the   matrix   A,   and  on   entry   with   
00088              TRANS = 'T' or 't' or 'C' or 'c',  K  specifies  the  number   
00089              of rows of the matrix  A.  K must be at least zero.   
00090              Unchanged on exit.   
00091     ALPHA  - DOUBLE PRECISION.   
00092              On entry, ALPHA specifies the scalar alpha.   
00093              Unchanged on exit.   
00094     A      - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is   
00095              k  when  TRANS = 'N' or 'n',  and is  n  otherwise.   
00096              Before entry with  TRANS = 'N' or 'n',  the  leading  n by k   
00097              part of the array  A  must contain the matrix  A,  otherwise   
00098              the leading  k by n  part of the array  A  must contain  the   
00099              matrix A.   
00100              Unchanged on exit.   
00101     LDA    - INTEGER.   
00102              On entry, LDA specifies the first dimension of A as declared   
00103              in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'   
00104              then  LDA must be at least  max( 1, n ), otherwise  LDA must   
00105              be at least  max( 1, k ).   
00106              Unchanged on exit.   
00107     BETA   - DOUBLE PRECISION.   
00108              On entry, BETA specifies the scalar beta.   
00109              Unchanged on exit.   
00110     C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ).   
00111              Before entry  with  UPLO = 'U' or 'u',  the leading  n by n   
00112              upper triangular part of the array C must contain the upper   
00113              triangular part  of the  symmetric matrix  and the strictly   
00114              lower triangular part of C is not referenced.  On exit, the   
00115              upper triangular part of the array  C is overwritten by the   
00116              upper triangular part of the updated matrix.   
00117              Before entry  with  UPLO = 'L' or 'l',  the leading  n by n   
00118              lower triangular part of the array C must contain the lower   
00119              triangular part  of the  symmetric matrix  and the strictly   
00120              upper triangular part of C is not referenced.  On exit, the   
00121              lower triangular part of the array  C is overwritten by the   
00122              lower triangular part of the updated matrix.   
00123     LDC    - INTEGER.   
00124              On entry, LDC specifies the first dimension of C as declared   
00125              in  the  calling  (sub)  program.   LDC  must  be  at  least   
00126              max( 1, n ).   
00127              Unchanged on exit.   
00128     Level 3 Blas routine.   
00129     -- Written on 8-February-1989.   
00130        Jack Dongarra, Argonne National Laboratory.   
00131        Iain Duff, AERE Harwell.   
00132        Jeremy Du Croz, Numerical Algorithms Group Ltd.   
00133        Sven Hammarling, Numerical Algorithms Group Ltd.   
00134        Test the input parameters.   
00135        Parameter adjustments */
00136     a_dim1 = *lda;
00137     a_offset = 1 + a_dim1 * 1;
00138     a -= a_offset;
00139     c_dim1 = *ldc;
00140     c_offset = 1 + c_dim1 * 1;
00141     c__ -= c_offset;
00142     /* Function Body */
00143     if (template_blas_lsame(trans, "N")) {
00144         nrowa = *n;
00145     } else {
00146         nrowa = *k;
00147     }
00148     upper = template_blas_lsame(uplo, "U");
00149     info = 0;
00150     if (! upper && ! template_blas_lsame(uplo, "L")) {
00151         info = 1;
00152     } else if (! template_blas_lsame(trans, "N") && ! template_blas_lsame(trans, 
00153             "T") && ! template_blas_lsame(trans, "C")) {
00154         info = 2;
00155     } else if (*n < 0) {
00156         info = 3;
00157     } else if (*k < 0) {
00158         info = 4;
00159     } else if (*lda < maxMACRO(1,nrowa)) {
00160         info = 7;
00161     } else if (*ldc < maxMACRO(1,*n)) {
00162         info = 10;
00163     }
00164     if (info != 0) {
00165         template_blas_erbla("DSYRK ", &info);
00166         return 0;
00167     }
00168 /*     Quick return if possible. */
00169     if (*n == 0 || ( (*alpha == 0. || *k == 0) && *beta == 1. ) ) {
00170         return 0;
00171     }
00172 /*     And when  alpha.eq.zero. */
00173     if (*alpha == 0.) {
00174         if (upper) {
00175             if (*beta == 0.) {
00176                 i__1 = *n;
00177                 for (j = 1; j <= i__1; ++j) {
00178                     i__2 = j;
00179                     for (i__ = 1; i__ <= i__2; ++i__) {
00180                         c___ref(i__, j) = 0.;
00181 /* L10: */
00182                     }
00183 /* L20: */
00184                 }
00185             } else {
00186                 i__1 = *n;
00187                 for (j = 1; j <= i__1; ++j) {
00188                     i__2 = j;
00189                     for (i__ = 1; i__ <= i__2; ++i__) {
00190                         c___ref(i__, j) = *beta * c___ref(i__, j);
00191 /* L30: */
00192                     }
00193 /* L40: */
00194                 }
00195             }
00196         } else {
00197             if (*beta == 0.) {
00198                 i__1 = *n;
00199                 for (j = 1; j <= i__1; ++j) {
00200                     i__2 = *n;
00201                     for (i__ = j; i__ <= i__2; ++i__) {
00202                         c___ref(i__, j) = 0.;
00203 /* L50: */
00204                     }
00205 /* L60: */
00206                 }
00207             } else {
00208                 i__1 = *n;
00209                 for (j = 1; j <= i__1; ++j) {
00210                     i__2 = *n;
00211                     for (i__ = j; i__ <= i__2; ++i__) {
00212                         c___ref(i__, j) = *beta * c___ref(i__, j);
00213 /* L70: */
00214                     }
00215 /* L80: */
00216                 }
00217             }
00218         }
00219         return 0;
00220     }
00221 /*     Start the operations. */
00222     if (template_blas_lsame(trans, "N")) {
00223 /*        Form  C := alpha*A*A' + beta*C. */
00224         if (upper) {
00225             i__1 = *n;
00226             for (j = 1; j <= i__1; ++j) {
00227                 if (*beta == 0.) {
00228                     i__2 = j;
00229                     for (i__ = 1; i__ <= i__2; ++i__) {
00230                         c___ref(i__, j) = 0.;
00231 /* L90: */
00232                     }
00233                 } else if (*beta != 1.) {
00234                     i__2 = j;
00235                     for (i__ = 1; i__ <= i__2; ++i__) {
00236                         c___ref(i__, j) = *beta * c___ref(i__, j);
00237 /* L100: */
00238                     }
00239                 }
00240                 i__2 = *k;
00241                 for (l = 1; l <= i__2; ++l) {
00242                     if (a_ref(j, l) != 0.) {
00243                         temp = *alpha * a_ref(j, l);
00244                         i__3 = j;
00245                         for (i__ = 1; i__ <= i__3; ++i__) {
00246                             c___ref(i__, j) = c___ref(i__, j) + temp * a_ref(
00247                                     i__, l);
00248 /* L110: */
00249                         }
00250                     }
00251 /* L120: */
00252                 }
00253 /* L130: */
00254             }
00255         } else {
00256             i__1 = *n;
00257             for (j = 1; j <= i__1; ++j) {
00258                 if (*beta == 0.) {
00259                     i__2 = *n;
00260                     for (i__ = j; i__ <= i__2; ++i__) {
00261                         c___ref(i__, j) = 0.;
00262 /* L140: */
00263                     }
00264                 } else if (*beta != 1.) {
00265                     i__2 = *n;
00266                     for (i__ = j; i__ <= i__2; ++i__) {
00267                         c___ref(i__, j) = *beta * c___ref(i__, j);
00268 /* L150: */
00269                     }
00270                 }
00271                 i__2 = *k;
00272                 for (l = 1; l <= i__2; ++l) {
00273                     if (a_ref(j, l) != 0.) {
00274                         temp = *alpha * a_ref(j, l);
00275                         i__3 = *n;
00276                         for (i__ = j; i__ <= i__3; ++i__) {
00277                             c___ref(i__, j) = c___ref(i__, j) + temp * a_ref(
00278                                     i__, l);
00279 /* L160: */
00280                         }
00281                     }
00282 /* L170: */
00283                 }
00284 /* L180: */
00285             }
00286         }
00287     } else {
00288 /*        Form  C := alpha*A'*A + beta*C. */
00289         if (upper) {
00290             i__1 = *n;
00291             for (j = 1; j <= i__1; ++j) {
00292                 i__2 = j;
00293                 for (i__ = 1; i__ <= i__2; ++i__) {
00294                     temp = 0.;
00295                     i__3 = *k;
00296                     for (l = 1; l <= i__3; ++l) {
00297                         temp += a_ref(l, i__) * a_ref(l, j);
00298 /* L190: */
00299                     }
00300                     if (*beta == 0.) {
00301                         c___ref(i__, j) = *alpha * temp;
00302                     } else {
00303                         c___ref(i__, j) = *alpha * temp + *beta * c___ref(i__,
00304                                  j);
00305                     }
00306 /* L200: */
00307                 }
00308 /* L210: */
00309             }
00310         } else {
00311             i__1 = *n;
00312             for (j = 1; j <= i__1; ++j) {
00313                 i__2 = *n;
00314                 for (i__ = j; i__ <= i__2; ++i__) {
00315                     temp = 0.;
00316                     i__3 = *k;
00317                     for (l = 1; l <= i__3; ++l) {
00318                         temp += a_ref(l, i__) * a_ref(l, j);
00319 /* L220: */
00320                     }
00321                     if (*beta == 0.) {
00322                         c___ref(i__, j) = *alpha * temp;
00323                     } else {
00324                         c___ref(i__, j) = *alpha * temp + *beta * c___ref(i__,
00325                                  j);
00326                     }
00327 /* L230: */
00328                 }
00329 /* L240: */
00330             }
00331         }
00332     }
00333     return 0;
00334 /*     End of DSYRK . */
00335 } /* dsyrk_ */
00336 #undef c___ref
00337 #undef a_ref
00338 
00339 #endif

Generated on Wed Nov 21 09:32:11 2012 for ergo by  doxygen 1.4.7