1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
30 /** @file matrix_algebra.cc
31
32 @brief A few matrix algebra routines for dense matrices.
33
34 @author: Elias Rudberg <em>responsible</em>
35 */
36
37 #include <stdlib.h>
38
39 #include "matrix_algebra.h"
40 #include "memorymanag.h"
41 #include "output.h"
42
43 #include "../matrix/mat_gblas.h"
44
45
46 #define USE_BLAS_MM
47
48
49 #if 0
50 #ifdef USE_BLAS_MM
51 #ifdef __cplusplus
52 extern "C"
53 #endif
54 void dgemm_(const char *ta,const char *tb,
55 const int *n, const int *k, const int *l,
56 const double *alpha,const double *A,const int *lda,
57 const double *B, const int *ldb,
58 const double *beta,const double *C, const int *ldc);
59 #endif
60 #endif
61
62
63 /*
64 Standard matrix multiplication.
65 */
66 void
multiply_matrices_general(int An1,int An2,int Bn1,int Bn2,const ergo_real * A,const ergo_real * B,ergo_real * AB)67 multiply_matrices_general(int An1, int An2, int Bn1, int Bn2, const ergo_real* A, const ergo_real* B, ergo_real* AB)
68 {
69 int i, j;
70 if(An2 != Bn1)
71 {
72 do_output(LOG_CAT_ERROR, LOG_AREA_LOWLEVEL,
73 "error in multiply_matrices_general: (An2 != Bn1)");
74 exit(0);
75 }
76 #ifdef USE_BLAS_MM
77 if(An1 == 0 || An2 == 0 || Bn1 == 0 || Bn2 == 0)
78 return;
79 /* call gemm */
80 ergo_real alpha = 1;
81 ergo_real beta = 0;
82 ergo_real* ABtemp = (ergo_real*)ergo_malloc(An1*Bn2*sizeof(ergo_real));
83 memset(ABtemp, 0, An1*Bn2*sizeof(ergo_real));
84 mat::gemm("T", "T", &An1, &Bn2, &An2, &alpha,
85 A, &An2,
86 B, &Bn2,
87 &beta,
88 ABtemp, &An1);
89 /* do transpose of result */
90 for(i = 0; i < An1; i++)
91 for(j = 0; j < Bn2; j++)
92 {
93 AB[i*Bn2+j] = ABtemp[j*An1+i];
94 }
95 ergo_free(ABtemp);
96 #else
97 for(i = 0; i < An1; i++)
98 for(j = 0; j < Bn2; j++)
99 {
100 ergo_real sum = 0;
101 for(int k = 0; k < An2; k++)
102 sum += A[i*An2+k] * B[k*Bn2+j];
103 AB[i*Bn2+j] = sum;
104 }
105 #endif
106 }
107
108
109 /*
110 Matrix multiplication when the first matrix is transposed.
111 */
112 void
multiply_matrices_general_T_1(int An1,int An2,int Bn1,int Bn2,const ergo_real * A,const ergo_real * B,ergo_real * AB)113 multiply_matrices_general_T_1(int An1, int An2, int Bn1, int Bn2, const ergo_real* A, const ergo_real* B, ergo_real* AB)
114 {
115 int i, j;
116 if(An1 != Bn1)
117 {
118 do_output(LOG_CAT_ERROR, LOG_AREA_LOWLEVEL, "error in multiply_matrices_general_T_1: (An1 != Bn1)");
119 exit(0);
120 }
121 ergo_real* At = (ergo_real*)ergo_malloc(An1*An2*sizeof(ergo_real));
122 for(i = 0; i < An1; i++)
123 for(j = 0; j < An2; j++)
124 At[j*An1+i] = A[i*An2+j];
125 multiply_matrices_general(An2, An1, Bn1, Bn2, At, B, AB);
126 ergo_free(At);
127 }
128