ergo
|
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