00001 /* Ergo, version 3.7, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2018 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, 00004 * and Anastasia Kruchinina. 00005 * 00006 * This program is free software: you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation, either version 3 of the License, or 00009 * (at your option) any later version. 00010 * 00011 * This program is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00018 * 00019 * Primary academic reference: 00020 * Ergo: An open-source program for linear-scaling electronic structure 00021 * calculations, 00022 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia 00023 * Kruchinina, 00024 * SoftwareX 7, 107 (2018), 00025 * <http://dx.doi.org/10.1016/j.softx.2018.03.005> 00026 * 00027 * For further information about Ergo, see <http://www.ergoscf.org>. 00028 */ 00029 00030 /* This file belongs to the template_lapack part of the Ergo source 00031 * code. The source files in the template_lapack directory are modified 00032 * versions of files originally distributed as CLAPACK, see the 00033 * Copyright/license notice in the file template_lapack/COPYING. 00034 */ 00035 00036 00037 #ifndef TEMPLATE_BLAS_SYMM_HEADER 00038 #define TEMPLATE_BLAS_SYMM_HEADER 00039 00040 00041 template<class Treal> 00042 int template_blas_symm(const char *side, const char *uplo, const integer *m, const integer *n, 00043 const Treal *alpha, const Treal *a, const integer *lda, const Treal *b, 00044 const integer *ldb, const Treal *beta, Treal *c__, const integer *ldc) 00045 { 00046 /* System generated locals */ 00047 integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2, 00048 i__3; 00049 /* Local variables */ 00050 integer info; 00051 Treal temp1, temp2; 00052 integer i__, j, k; 00053 integer nrowa; 00054 logical upper; 00055 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 00056 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1] 00057 #define c___ref(a_1,a_2) c__[(a_2)*c_dim1 + a_1] 00058 /* Purpose 00059 ======= 00060 DSYMM performs one of the matrix-matrix operations 00061 C := alpha*A*B + beta*C, 00062 or 00063 C := alpha*B*A + beta*C, 00064 where alpha and beta are scalars, A is a symmetric matrix and B and 00065 C are m by n matrices. 00066 Parameters 00067 ========== 00068 SIDE - CHARACTER*1. 00069 On entry, SIDE specifies whether the symmetric matrix A 00070 appears on the left or right in the operation as follows: 00071 SIDE = 'L' or 'l' C := alpha*A*B + beta*C, 00072 SIDE = 'R' or 'r' C := alpha*B*A + beta*C, 00073 Unchanged on exit. 00074 UPLO - CHARACTER*1. 00075 On entry, UPLO specifies whether the upper or lower 00076 triangular part of the symmetric matrix A is to be 00077 referenced as follows: 00078 UPLO = 'U' or 'u' Only the upper triangular part of the 00079 symmetric matrix is to be referenced. 00080 UPLO = 'L' or 'l' Only the lower triangular part of the 00081 symmetric matrix is to be referenced. 00082 Unchanged on exit. 00083 M - INTEGER. 00084 On entry, M specifies the number of rows of the matrix C. 00085 M must be at least zero. 00086 Unchanged on exit. 00087 N - INTEGER. 00088 On entry, N specifies the number of columns of the matrix C. 00089 N 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 m when SIDE = 'L' or 'l' and is n otherwise. 00096 Before entry with SIDE = 'L' or 'l', the m by m part of 00097 the array A must contain the symmetric matrix, such that 00098 when UPLO = 'U' or 'u', the leading m by m upper triangular 00099 part of the array A must contain the upper triangular part 00100 of the symmetric matrix and the strictly lower triangular 00101 part of A is not referenced, and when UPLO = 'L' or 'l', 00102 the leading m by m lower triangular part of the array A 00103 must contain the lower triangular part of the symmetric 00104 matrix and the strictly upper triangular part of A is not 00105 referenced. 00106 Before entry with SIDE = 'R' or 'r', the n by n part of 00107 the array A must contain the symmetric matrix, such that 00108 when UPLO = 'U' or 'u', the leading n by n upper triangular 00109 part of the array A must contain the upper triangular part 00110 of the symmetric matrix and the strictly lower triangular 00111 part of A is not referenced, and when UPLO = 'L' or 'l', 00112 the leading n by n lower triangular part of the array A 00113 must contain the lower triangular part of the symmetric 00114 matrix and the strictly upper triangular part of A is not 00115 referenced. 00116 Unchanged on exit. 00117 LDA - INTEGER. 00118 On entry, LDA specifies the first dimension of A as declared 00119 in the calling (sub) program. When SIDE = 'L' or 'l' then 00120 LDA must be at least max( 1, m ), otherwise LDA must be at 00121 least max( 1, n ). 00122 Unchanged on exit. 00123 B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). 00124 Before entry, the leading m by n part of the array B must 00125 contain the matrix B. 00126 Unchanged on exit. 00127 LDB - INTEGER. 00128 On entry, LDB specifies the first dimension of B as declared 00129 in the calling (sub) program. LDB must be at least 00130 max( 1, m ). 00131 Unchanged on exit. 00132 BETA - DOUBLE PRECISION. 00133 On entry, BETA specifies the scalar beta. When BETA is 00134 supplied as zero then C need not be set on input. 00135 Unchanged on exit. 00136 C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). 00137 Before entry, the leading m by n part of the array C must 00138 contain the matrix C, except when beta is zero, in which 00139 case C need not be set on entry. 00140 On exit, the array C is overwritten by the m by n updated 00141 matrix. 00142 LDC - INTEGER. 00143 On entry, LDC specifies the first dimension of C as declared 00144 in the calling (sub) program. LDC must be at least 00145 max( 1, m ). 00146 Unchanged on exit. 00147 Level 3 Blas routine. 00148 -- Written on 8-February-1989. 00149 Jack Dongarra, Argonne National Laboratory. 00150 Iain Duff, AERE Harwell. 00151 Jeremy Du Croz, Numerical Algorithms Group Ltd. 00152 Sven Hammarling, Numerical Algorithms Group Ltd. 00153 Set NROWA as the number of rows of A. 00154 Parameter adjustments */ 00155 a_dim1 = *lda; 00156 a_offset = 1 + a_dim1 * 1; 00157 a -= a_offset; 00158 b_dim1 = *ldb; 00159 b_offset = 1 + b_dim1 * 1; 00160 b -= b_offset; 00161 c_dim1 = *ldc; 00162 c_offset = 1 + c_dim1 * 1; 00163 c__ -= c_offset; 00164 /* Function Body */ 00165 if (template_blas_lsame(side, "L")) { 00166 nrowa = *m; 00167 } else { 00168 nrowa = *n; 00169 } 00170 upper = template_blas_lsame(uplo, "U"); 00171 /* Test the input parameters. */ 00172 info = 0; 00173 if (! template_blas_lsame(side, "L") && ! template_blas_lsame(side, "R")) { 00174 info = 1; 00175 } else if (! upper && ! template_blas_lsame(uplo, "L")) { 00176 info = 2; 00177 } else if (*m < 0) { 00178 info = 3; 00179 } else if (*n < 0) { 00180 info = 4; 00181 } else if (*lda < maxMACRO(1,nrowa)) { 00182 info = 7; 00183 } else if (*ldb < maxMACRO(1,*m)) { 00184 info = 9; 00185 } else if (*ldc < maxMACRO(1,*m)) { 00186 info = 12; 00187 } 00188 if (info != 0) { 00189 template_blas_erbla("SYMM ", &info); 00190 return 0; 00191 } 00192 /* Quick return if possible. */ 00193 if (*m == 0 || *n == 0 || ( *alpha == 0. && *beta == 1. ) ) { 00194 return 0; 00195 } 00196 /* And when alpha.eq.zero. */ 00197 if (*alpha == 0.) { 00198 if (*beta == 0.) { 00199 i__1 = *n; 00200 for (j = 1; j <= i__1; ++j) { 00201 i__2 = *m; 00202 for (i__ = 1; i__ <= i__2; ++i__) { 00203 c___ref(i__, j) = 0.; 00204 /* L10: */ 00205 } 00206 /* L20: */ 00207 } 00208 } else { 00209 i__1 = *n; 00210 for (j = 1; j <= i__1; ++j) { 00211 i__2 = *m; 00212 for (i__ = 1; i__ <= i__2; ++i__) { 00213 c___ref(i__, j) = *beta * c___ref(i__, j); 00214 /* L30: */ 00215 } 00216 /* L40: */ 00217 } 00218 } 00219 return 0; 00220 } 00221 /* Start the operations. */ 00222 if (template_blas_lsame(side, "L")) { 00223 /* Form C := alpha*A*B + beta*C. */ 00224 if (upper) { 00225 i__1 = *n; 00226 for (j = 1; j <= i__1; ++j) { 00227 i__2 = *m; 00228 for (i__ = 1; i__ <= i__2; ++i__) { 00229 temp1 = *alpha * b_ref(i__, j); 00230 temp2 = 0.; 00231 i__3 = i__ - 1; 00232 for (k = 1; k <= i__3; ++k) { 00233 c___ref(k, j) = c___ref(k, j) + temp1 * a_ref(k, i__); 00234 temp2 += b_ref(k, j) * a_ref(k, i__); 00235 /* L50: */ 00236 } 00237 if (*beta == 0.) { 00238 c___ref(i__, j) = temp1 * a_ref(i__, i__) + *alpha * 00239 temp2; 00240 } else { 00241 c___ref(i__, j) = *beta * c___ref(i__, j) + temp1 * 00242 a_ref(i__, i__) + *alpha * temp2; 00243 } 00244 /* L60: */ 00245 } 00246 /* L70: */ 00247 } 00248 } else { 00249 i__1 = *n; 00250 for (j = 1; j <= i__1; ++j) { 00251 for (i__ = *m; i__ >= 1; --i__) { 00252 temp1 = *alpha * b_ref(i__, j); 00253 temp2 = 0.; 00254 i__2 = *m; 00255 for (k = i__ + 1; k <= i__2; ++k) { 00256 c___ref(k, j) = c___ref(k, j) + temp1 * a_ref(k, i__); 00257 temp2 += b_ref(k, j) * a_ref(k, i__); 00258 /* L80: */ 00259 } 00260 if (*beta == 0.) { 00261 c___ref(i__, j) = temp1 * a_ref(i__, i__) + *alpha * 00262 temp2; 00263 } else { 00264 c___ref(i__, j) = *beta * c___ref(i__, j) + temp1 * 00265 a_ref(i__, i__) + *alpha * temp2; 00266 } 00267 /* L90: */ 00268 } 00269 /* L100: */ 00270 } 00271 } 00272 } else { 00273 /* Form C := alpha*B*A + beta*C. */ 00274 i__1 = *n; 00275 for (j = 1; j <= i__1; ++j) { 00276 temp1 = *alpha * a_ref(j, j); 00277 if (*beta == 0.) { 00278 i__2 = *m; 00279 for (i__ = 1; i__ <= i__2; ++i__) { 00280 c___ref(i__, j) = temp1 * b_ref(i__, j); 00281 /* L110: */ 00282 } 00283 } else { 00284 i__2 = *m; 00285 for (i__ = 1; i__ <= i__2; ++i__) { 00286 c___ref(i__, j) = *beta * c___ref(i__, j) + temp1 * b_ref( 00287 i__, j); 00288 /* L120: */ 00289 } 00290 } 00291 i__2 = j - 1; 00292 for (k = 1; k <= i__2; ++k) { 00293 if (upper) { 00294 temp1 = *alpha * a_ref(k, j); 00295 } else { 00296 temp1 = *alpha * a_ref(j, k); 00297 } 00298 i__3 = *m; 00299 for (i__ = 1; i__ <= i__3; ++i__) { 00300 c___ref(i__, j) = c___ref(i__, j) + temp1 * b_ref(i__, k); 00301 /* L130: */ 00302 } 00303 /* L140: */ 00304 } 00305 i__2 = *n; 00306 for (k = j + 1; k <= i__2; ++k) { 00307 if (upper) { 00308 temp1 = *alpha * a_ref(j, k); 00309 } else { 00310 temp1 = *alpha * a_ref(k, j); 00311 } 00312 i__3 = *m; 00313 for (i__ = 1; i__ <= i__3; ++i__) { 00314 c___ref(i__, j) = c___ref(i__, j) + temp1 * b_ref(i__, k); 00315 /* L150: */ 00316 } 00317 /* L160: */ 00318 } 00319 /* L170: */ 00320 } 00321 } 00322 return 0; 00323 /* End of DSYMM . */ 00324 } /* dsymm_ */ 00325 #undef c___ref 00326 #undef b_ref 00327 #undef a_ref 00328 00329 #endif