1 #pragma once
2 // blas_l3.hpp: BLAS Level 3 functions
3 //
4 // Copyright (C) 2017-2021 Stillwater Supercomputing, Inc.
5 //
6 // This file is part of the universal numbers project, which is released under an MIT Open Source license.
7
8 namespace sw::universal::blas {
9
10 // sum entire matrix (dim == 0), all rows (dim == 1), or all columns (dim == 2)
11 template<typename Matrix>
sum(Matrix & A,int dim=0)12 vector<typename Matrix::value_type> sum(Matrix& A, int dim = 0) {
13 using value_type = typename Matrix::value_type;
14 using size_type = typename Matrix::size_type;
15
16 size_type rows = num_rows(A);
17 size_type cols = num_cols(A);
18
19 switch (dim) {
20 case 0:
21 {
22 value_type sum{ 0 };
23 for (size_type i = 0; i < rows; ++i) {
24 for (size_type j = 0; j < cols; ++j) {
25 sum += A(i, j);
26 }
27 }
28 return vector<value_type>{sum};
29 }
30
31 case 1:
32 {
33 vector<value_type> rowSums(rows);
34 for (size_type i = 0; i < rows; ++i) {
35 for (size_type j = 0; j < cols; ++j) {
36 rowSums[i] += A(i, j);
37 }
38 }
39 return rowSums;
40 }
41
42 case 2:
43 {
44 vector<value_type> colSums(cols);
45 for (size_type i = 0; i < rows; ++i) {
46 for (size_type j = 0; j < cols; ++j) {
47 colSums[j] += A(i, j);
48 }
49 }
50 return colSums;
51 }
52
53 default:
54 break;
55 }
56
57 return vector<value_type>{0};
58 }
59
60 enum NormalizationMethod {
61 Norm2,
62 Center,
63 Zscore,
64 Scale,
65 Range
66 };
67
68 // normalize entire matrix (dim == 0), all rows (dim == 1), or all columns (dim == 2)
69 template<typename Matrix>
normalize(Matrix & A,int dim=0)70 void normalize(Matrix& A, int dim = 0) {
71 using value_type = typename Matrix::value_type;
72 using size_type = typename Matrix::size_type;
73
74 size_type rows = num_rows(A);
75 size_type cols = num_cols(A);
76
77 switch (dim) {
78 case 0:
79 {
80 value_type sos{ 0 };
81 for (size_type i = 0; i < rows; ++i) {
82 for (size_type j = 0; j < cols; ++j) {
83 sos += A(i, j) * A(i, j);
84 }
85 }
86 for (size_type i = 0; i < rows; ++i) {
87 for (size_type j = 0; j < cols; ++j) {
88 A(i, j) /= sqrt(sos);
89 }
90 }
91 }
92 break;
93 case 1:
94 {
95 vector<value_type> rowSos(rows);
96 for (size_type i = 0; i < rows; ++i) {
97 for (size_type j = 0; j < cols; ++j) {
98 rowSos[i] += A(i, j) * A(i, j);
99 }
100 }
101 for (size_type i = 0; i < rows; ++i) {
102 for (size_type j = 0; j < cols; ++j) {
103 A(i, j) /= sqrt(rowSos[i]);
104 }
105 }
106 }
107 break;
108 case 2:
109 {
110 vector<value_type> colSos(cols);
111 for (size_type i = 0; i < rows; ++i) {
112 for (size_type j = 0; j < cols; ++j) {
113 colSos[j] += A(i, j) * A(i, j);
114 }
115 }
116 for (size_type i = 0; i < rows; ++i) {
117 for (size_type j = 0; j < cols; ++j) {
118 A(i, j) /= sqrt(colSos[j]);
119 }
120 }
121 }
122 break;
123 default:
124 break;
125 }
126 }
127
128 // 2-norm of entire matrix (dim == 0), each rows (dim == 1), or each columns (dim == 2)
129 template<typename Matrix>
matrixNorm(Matrix & A,int dim=0)130 vector<typename Matrix::value_type> matrixNorm(Matrix& A, int dim = 0) {
131 using value_type = typename Matrix::value_type;
132 using size_type = typename Matrix::size_type;
133
134 size_type rows = num_rows(A);
135 size_type cols = num_cols(A);
136
137 switch (dim) {
138 case 0:
139 {
140 value_type sos{ 0 };
141 for (size_type i = 0; i < rows; ++i) {
142 for (size_type j = 0; j < cols; ++j) {
143 sos += A(i, j) * A(i,j);
144 }
145 }
146 return vector<value_type>{sqrt(sos)};
147 }
148
149 case 1:
150 {
151 vector<value_type> rowSoS(rows);
152 for (size_type i = 0; i < rows; ++i) {
153 for (size_type j = 0; j < cols; ++j) {
154 rowSoS[i] += A(i, j) * A(i, j);
155 }
156 rowSoS[i] = sqrt(rowSoS[i]);
157 }
158 return rowSoS;
159 }
160
161 case 2:
162 {
163 vector<value_type> colSoS(cols);
164 for (size_type i = 0; i < rows; ++i) {
165 for (size_type j = 0; j < cols; ++j) {
166 colSoS[j] += A(i, j) * A(i, j);
167 }
168 }
169 for (size_type j = 0; j < cols; ++j) {
170 colSoS[j] = sqrt(colSoS[j]);
171 }
172 return colSoS;
173 }
174
175 default:
176 break;
177 }
178 return vector<value_type>{-1};
179 }
180
181 }
182