MiniTensor Version of the Day
Loading...
Searching...
No Matches
MiniTensor_Matrix.h
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// MiniTensor Package
4//
5// Copyright 2016 NTESS and the MiniTensor contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#if !defined(MiniTensor_Matrix_h)
11#define MiniTensor_Matrix_h
12
13#include <algorithm>
14#include <cassert>
15#include <iomanip>
16#include <iostream>
17#include <vector>
18
19#include <Kokkos_Core.hpp>
20#include "MiniTensor_Vector.h"
21
22namespace minitensor {
23
26
30template<typename T, Index M, Index N>
32
36template<typename T, Index M = DYNAMIC, Index N = DYNAMIC>
37class Matrix: public TensorBase<T, matrix_store<T, M, N>>
38{
39public:
40
44 static constexpr
45 Index
46 ORDER = 1;
47
51 static constexpr
52 bool
53 IS_DYNAMIC = (M == DYNAMIC || N == DYNAMIC);
54
59
64 static constexpr
65 Index
67 {
68 return ORDER;
69 }
70
75 explicit
76 Matrix();
77
84 explicit
85 Matrix(Index const rows, Index const cols);
86
92 explicit
93 Matrix(Filler const value);
94
102 explicit
103 Matrix(Index const rows, Index cols, Filler const value);
104
111 template<class ArrayT>
113 Matrix(
114 Source const source,
115 ArrayT & data,
116 Index index1);
117
125 template<class ArrayT>
127 Matrix(
128 Source const source,
129 ArrayT & data,
130 Index index1,
131 Index index2);
132
141 template<class ArrayT>
143 Matrix(
144 Source const source,
145 ArrayT & data,
146 Index index1,
147 Index index2,
148 Index index3);
149
159 template<class ArrayT>
161 Matrix(
162 Source const source,
163 ArrayT & data,
164 Index index1,
165 Index index2,
166 Index index3,
167 Index index4);
168
179 template<class ArrayT>
181 Matrix(
182 Source const source,
183 ArrayT & data,
184 Index index1,
185 Index index2,
186 Index index3,
187 Index index4,
188 Index index5);
189
201 template<class ArrayT>
203 Matrix(
204 Source const source,
205 ArrayT & data,
206 Index index1,
207 Index index2,
208 Index index3,
209 Index index4,
210 Index index5,
211 Index index6);
212
221 template<class ArrayT>
223 Matrix(
224 Source const source,
225 Index const rows,
226 Index const cols,
227 ArrayT & data,
228 Index index1);
229
239 template<class ArrayT>
241 Matrix(
242 Source const source,
243 Index const rows,
244 Index const cols,
245 ArrayT & data,
246 Index index1,
247 Index index2);
248
259 template<class ArrayT>
261 Matrix(
262 Source const source,
263 Index const rows,
264 Index const cols,
265 ArrayT & data,
266 Index index1,
267 Index index2,
268 Index index3);
269
281 template<class ArrayT>
283 Matrix(
284 Source const source,
285 Index const rows,
286 Index const cols,
287 ArrayT & data,
288 Index index1,
289 Index index2,
290 Index index3,
291 Index index4);
292
305 template<class ArrayT>
307 Matrix(
308 Source const source,
309 Index const rows,
310 Index const cols,
311 ArrayT & data,
312 Index index1,
313 Index index2,
314 Index index3,
315 Index index4,
316 Index index5);
317
331 template<class ArrayT>
333 Matrix(
334 Source const source,
335 Index const rows,
336 Index const cols,
337 ArrayT & data,
338 Index index1,
339 Index index2,
340 Index index3,
341 Index index4,
342 Index index5,
343 Index index6);
344
350 explicit
351 Matrix(T const * data_ptr);
352
360 explicit
361 Matrix(Index const rows, Index cols, T const * data_ptr);
362
367 Matrix(Matrix<T, M, N> const & A);
368
372 virtual
374 ~Matrix();
375
382 T const &
383 operator()(Index const i, Index const j) const;
384
391 T &
392 operator()(Index const i, Index const j);
393
397 std::pair<Index, Index>
398 get_dimensions() const;
399
404 Index
405 get_num_rows() const;
406
411 Index
412 get_num_cols() const;
413
419 void
420 set_dimensions(Index const rows, Index const cols);
421
422private:
423
427 Index
429
433 Index
435
436};
437
442template<typename S, typename T, Index M, Index N>
444Matrix<typename Promote<S, T>::type, M, N>
445operator+(Matrix<S, M, N> const & A, Matrix<T, M, N> const & B);
446
451template<typename S, typename T, Index M, Index N>
453Matrix<typename Promote<S, T>::type, M, N>
454operator-(Matrix<S, M, N> const & A, Matrix<T, M, N> const & B);
455
460template<typename T, Index M, Index N>
462Matrix<T, M, N>
463operator-(Matrix<T, M, N> const & A);
464
470template<typename T, Index M, Index N>
472bool
473operator==(Matrix<T, M, N> const & A, Matrix<T, M, N> const & B);
474
480template<typename T, Index M, Index N>
482bool
483operator!=(Matrix<T, M, N> const & A, Matrix<T, M, N> const & B);
484
491template<typename S, typename T, Index M, Index N>
493Vector<typename Promote<S, T>::type, M>
494operator*(Matrix<T, M, N> const & A, Vector<S, N> const & u);
495
502template<typename S, typename T, Index M, Index N>
504Vector<typename Promote<S, T>::type, N>
505operator*(Vector<S, M> const & u, Matrix<T, M, N> const & A);
506
513template<typename S, typename T, Index M, Index N>
515Matrix<typename Promote<S, T>::type, M, N>
516operator*(Matrix<T, M, N> const & A, Tensor<S, N> const & B);
517
524template<typename S, typename T, Index M, Index N>
526Matrix<typename Promote<S, T>::type, M, N>
527operator*(Tensor<S, M> const & A, Matrix<T, M, N> const & B);
528
533template<typename S, typename T, Index M, Index P, Index N>
535Matrix<typename Promote<S, T>::type, M, N>
536operator*(Matrix<S, M, P> const & A, Matrix<T, P, N> const & B);
537
544template<typename S, typename T, Index M, Index N>
546typename lazy_disable_if<order_1234<S>, apply_matrix<Promote<S, T>, M, N>>::type
547operator*(S const & s, Matrix<T, M, N> const & A);
548
555template<typename S, typename T, Index M, Index N>
557typename lazy_disable_if<order_1234<S>, apply_matrix<Promote<S, T>, M, N>>::type
558operator*(Matrix<T, M, N> const & A, S const & s);
559
566template<typename S, typename T, Index M, Index N>
568Matrix<typename Promote<S, T>::type, M, N>
569operator/(Matrix<T, M, N> const & A, S const & s);
570
577template<typename S, typename T, Index M, Index N>
579Matrix<typename Promote<S, T>::type, M, N>
580operator/(S const & s, Matrix<T, M, N> const & A);
581
588template<typename T, Index M, Index N>
589std::istream &
590operator>>(std::istream & is, Matrix<T, M, N> & A);
591
598template<typename T, Index M, Index N>
599std::ostream &
600operator<<(std::ostream & os, Matrix<T, M, N> const & A);
601
608template<typename T, Index M, Index N>
610Vector<T, N>
611row(Matrix<T, M, N> const & A, Index const i);
612
619template<typename T, Index M, Index N>
621Vector<T, M>
622col(Matrix<T, M, N> const & A, Index const j);
623
630template<typename S, typename T, Index M, Index N>
632Vector<typename Promote<S, T>::type, M>
633dot(Matrix<T, M, N> const & A, Vector<S, N> const & u);
634
641template<typename S, typename T, Index M, Index N>
643Vector<typename Promote<S, T>::type, N>
644dot(Vector<S, M> const & u, Matrix<T, M, N> const & A);
645
652template<typename S, typename T, Index M, Index N>
654Matrix<typename Promote<S, T>::type, M, N>
655dot(Matrix<T, M, N> const & A, Tensor<S, N> const & B);
656
663template<typename S, typename T, Index M, Index N>
665Matrix<typename Promote<S, T>::type, M, N>
666dot(Tensor<S, M> const & A, Matrix<T, M, N> const & B);
667
674template<typename S, typename T, Index M, Index P, Index N>
676Matrix<typename Promote<S, T>::type, M, N>
677dot(Matrix<S, M, P> const & A, Matrix<T, P, N> const & B);
678
685template<typename S, typename T, Index M, Index P, Index N>
687Matrix<typename Promote<S, T>::type, M, N>
688t_dot(Matrix<S, P, M> const & A, Matrix<T, P, N> const & B);
689
696template<typename S, typename T, Index M, Index P, Index N>
698Matrix<typename Promote<S, T>::type, M, N>
699dot_t(Matrix<S, M, P> const & A, Matrix<T, N, P> const & B);
700
707template<typename S, typename T, Index M, Index P, Index N>
709Matrix<typename Promote<S, T>::type, M, N>
710t_dot_t(Matrix<S, P, M> const & A, Matrix<T, N, P> const & B);
711
718template<typename S, typename T, Index M, Index N>
720typename Promote<S, T>::type
721dotdot(Matrix<S, M, N> const & A, Matrix<T, M, N> const & B);
722
729template<typename S, typename T, Index M, Index N>
731Matrix<typename Promote<S, T>::type, M, N>
732matrix(Vector<S, M> const & u, Vector<T, N> const & v);
733
737template<typename T, Index M, Index N>
739Matrix<T, M, N>
740transpose(Matrix<T, N, M> const & A);
741
745template<typename T, Index M, Index N>
747Matrix<T, M, N>
749
750} // namespace minitensor
751
752namespace minitensor {
753
754//
755// Constructor that initializes to NaNs
756//
757template<typename T, Index M, Index N>
761{
762 set_dimensions(M, N);
763 return;
764}
765
766template<typename T, Index M, Index N>
768Matrix<T, M, N>::Matrix(Index const rows, Index const cols) :
769 TensorBase<T, Store>::TensorBase(rows * cols, ORDER),
770 rows_(rows), cols_(cols)
771{
772 return;
773}
774
775//
776// Create matrix from a specified value
777//
778template<typename T, Index M, Index N>
781 TensorBase<T, Store>::TensorBase(M * N, ORDER, value)
782{
783 return;
784}
785
786template<typename T, Index M, Index N>
789 Index const rows,
790 Index const cols ,
791 Filler const value) :
792 TensorBase<T, Store>::TensorBase(rows * cols, ORDER, value),
793 rows_(rows), cols_(cols)
794{
795 return;
796}
797
798//
799// Create matrix from array
800//
801//
802template<typename T, Index M, Index N>
803template<class ArrayT>
806 Source const,
807 ArrayT & data,
808 Index index1) :
809 TensorBase<T, Store>::TensorBase(M * N, ORDER, data, index1)
810{
811 return;
812}
813
814template<typename T, Index M, Index N>
815template<class ArrayT>
818 Source const,
819 ArrayT & data,
820 Index index1,
821 Index index2) :
822 TensorBase<T, Store>::TensorBase(M * N, ORDER, data, index1, index2)
823{
824 return;
825}
826
827template<typename T, Index M, Index N>
828template<class ArrayT>
831 Source const,
832 ArrayT & data,
833 Index index1,
834 Index index2,
835 Index index3) :
836 TensorBase<T, Store>::TensorBase(M * N, ORDER, data, index1, index2, index3)
837{
838 return;
839}
840
841template<typename T, Index M, Index N>
842template<class ArrayT>
845 Source const,
846 ArrayT & data,
847 Index index1,
848 Index index2,
849 Index index3,
850 Index index4) :
852 M * N,
853 ORDER,
854 data,
855 index1,
856 index2,
857 index3,
858 index4)
859{
860 return;
861}
862
863template<typename T, Index M, Index N>
864template<class ArrayT>
867 Source const,
868 ArrayT & data,
869 Index index1,
870 Index index2,
871 Index index3,
872 Index index4,
873 Index index5) :
875 M * N,
876 ORDER,
877 data,
878 index1,
879 index2,
880 index3,
881 index4,
882 index5)
883{
884 return;
885}
886
887template<typename T, Index M, Index N>
888template<class ArrayT>
891 Source const,
892 ArrayT & data,
893 Index index1,
894 Index index2,
895 Index index3,
896 Index index4,
897 Index index5,
898 Index index6) :
900 M * N,
901 ORDER,
902 data,
903 index1,
904 index2,
905 index3,
906 index4,
907 index5,
908 index6)
909{
910 return;
911}
912
913template<typename T, Index M, Index N>
914template<class ArrayT>
917 Source const,
918 Index const rows,
919 Index const cols,
920 ArrayT & data,
921 Index index1) :
922 TensorBase<T, Store>::TensorBase(rows * cols, ORDER, data, index1),
923 rows_(rows), cols_(cols)
924{
925 return;
926}
927
928template<typename T, Index M, Index N>
929template<class ArrayT>
932 Source const,
933 Index const rows,
934 Index const cols,
935 ArrayT & data,
936 Index index1,
937 Index index2) :
938 TensorBase<T, Store>::TensorBase(rows * cols, ORDER, data, index1, index2),
939 rows_(rows), cols_(cols)
940{
941 return;
942}
943
944template<typename T, Index M, Index N>
945template<class ArrayT>
948 Source const,
949 Index const rows,
950 Index const cols,
951 ArrayT & data,
952 Index index1,
953 Index index2,
954 Index index3) :
956 rows * cols,
957 ORDER,
958 data,
959 index1,
960 index2,
961 index3),
962 rows_(rows), cols_(cols)
963{
964 return;
965}
966
967template<typename T, Index M, Index N>
968template<class ArrayT>
971 Source const,
972 Index const rows,
973 Index const cols,
974 ArrayT & data,
975 Index index1,
976 Index index2,
977 Index index3,
978 Index index4) :
980 rows * cols,
981 ORDER,
982 data,
983 index1,
984 index2,
985 index3,
986 index4),
987 rows_(rows), cols_(cols)
988{
989 return;
990}
991
992template<typename T, Index M, Index N>
993template<class ArrayT>
996 Source const,
997 Index const rows,
998 Index const cols,
999 ArrayT & data,
1000 Index index1,
1001 Index index2,
1002 Index index3,
1003 Index index4,
1004 Index index5) :
1006 rows * cols,
1007 ORDER,
1008 data,
1009 index1,
1010 index2,
1011 index3,
1012 index4,
1013 index5),
1014 rows_(rows), cols_(cols)
1015{
1016 return;
1017}
1018
1019template<typename T, Index M, Index N>
1020template<class ArrayT>
1023 Source const,
1024 Index const rows,
1025 Index const cols,
1026 ArrayT & data,
1027 Index index1,
1028 Index index2,
1029 Index index3,
1030 Index index4,
1031 Index index5,
1032 Index index6) :
1034 rows * cols,
1035 ORDER,
1036 data,
1037 index1,
1038 index2,
1039 index3,
1040 index4,
1041 index5,
1042 index6),
1043 rows_(rows), cols_(cols)
1044{
1045 return;
1046}
1047
1048template<typename T, Index M, Index N>
1050Matrix<T, M, N>::Matrix(T const * data_ptr) :
1051 TensorBase<T, Store>::TensorBase(M * N, ORDER, data_ptr)
1052{
1053 return;
1054}
1055
1056template<typename T, Index M, Index N>
1059 Index const rows,
1060 Index const cols,
1061 T const * data_ptr) :
1062 TensorBase<T, Store>::TensorBase(rows * cols, ORDER, data_ptr),
1063 rows_(rows), cols_(cols)
1064{
1065 return;
1066}
1067
1068//
1069// Copy constructor
1070//
1071template<typename T, Index M, Index N>
1074 TensorBase<T, Store>::TensorBase(A),
1075 rows_(A.get_num_rows()), cols_(A.get_num_cols())
1076{
1077 return;
1078}
1079
1080//
1081// Simple destructor
1082//
1083template<typename T, Index M, Index N>
1086{
1087 return;
1088}
1089
1090//
1091// Get dimensions
1092//
1093template<typename T, Index M, Index N>
1094inline
1095std::pair<Index, Index>
1097{
1098 return std::make_pair(rows_, cols_);
1099}
1100
1101//
1102// Get number rows
1103//
1104template<typename T, Index M, Index N>
1106Index
1108{
1109 return rows_;
1110}
1111
1112//
1113// Get number cols
1114//
1115template<typename T, Index M, Index N>
1117Index
1119{
1120 return cols_;
1121}
1122
1123//
1124// Set dimensions
1125//
1126template<typename T, Index M, Index N>
1128void
1130{
1131 TensorBase<T, Store>::set_dimension(rows * cols, ORDER);
1132 rows_ = rows;
1133 cols_ = cols;
1134
1135 return;
1136}
1137
1138//
1139// Indexing for constant matrix
1140//
1141template<typename T, Index M, Index N>
1143T const &
1145{
1146 Matrix<T, M, N> const &
1147 self = (*this);
1148
1149 Index const
1150 num_cols = self.get_num_cols();
1151
1152 return self[i * num_cols + j];
1153}
1154
1155//
1156// Matrix indexing
1157//
1158template<typename T, Index M, Index N>
1160T &
1162{
1164 self = (*this);
1165
1166 Index const
1167 num_cols = self.get_num_cols();
1168
1169 return self[i * num_cols + j];
1170}
1171
1172//
1173// Extract a row as a vector
1174//
1175template<typename T, Index M, Index N>
1178row(Matrix<T, M, N> const & A, Index const i)
1179{
1180 Index const
1181 num_cols = A.get_num_cols();
1182
1184 v(num_cols);
1185
1186 switch (num_cols) {
1187 default:
1188 for (Index j = 0; j < num_cols; ++j) {
1189 v(j) = A(i, j);
1190 }
1191 break;
1192
1193 case 2:
1194 v(0) = A(i, 0);
1195 v(1) = A(i, 1);
1196 break;
1197
1198 case 3:
1199 v(0) = A(i, 0);
1200 v(1) = A(i, 1);
1201 v(2) = A(i, 2);
1202 break;
1203 }
1204
1205 return v;
1206}
1207
1208//
1209// Extract a column as a vector
1210//
1211template<typename T, Index M, Index N>
1213Vector<T, M>
1214col(Matrix<T, M, N> const & A, Index const j)
1215{
1216 Index const
1217 num_rows = A.get_num_rows();
1218
1220 v(num_rows);
1221
1222 switch (num_rows) {
1223 default:
1224 for (Index i = 0; i < num_rows; ++i) {
1225 v(i) = A(i, j);
1226 }
1227 break;
1228
1229 case 2:
1230 v(0) = A(0, j);
1231 v(1) = A(1, j);
1232 break;
1233
1234 case 3:
1235 v(0) = A(0, j);
1236 v(1) = A(1, j);
1237 v(2) = A(2, j);
1238 break;
1239 }
1240
1241 return v;
1242}
1243
1244//
1245// Matrix addition
1246//
1247template<typename S, typename T, Index M, Index N>
1249Matrix<typename Promote<S, T>::type, M, N>
1251{
1253 C(A.get_num_rows(), A.get_num_cols());
1254
1255 add(A, B, C);
1256
1257 return C;
1258}
1259
1260//
1261// Matrix subtraction
1262//
1263template<typename S, typename T, Index M, Index N>
1265Matrix<typename Promote<S, T>::type, M, N>
1267{
1269 C(A.get_num_rows(), A.get_num_cols());
1270
1271 subtract(A, B, C);
1272
1273 return C;
1274}
1275
1276//
1277// Matrix minus
1278//
1279template<typename T, Index M, Index N>
1281Matrix<T, M, N>
1283{
1285 B(A.get_num_rows(), A.get_num_cols());
1286
1287 minus(A, B);
1288
1289 return B;
1290}
1291
1292//
1293// Matrix equality
1294//
1295template<typename T, Index M, Index N>
1297bool
1299{
1300 return equal(A, B);
1301}
1302
1303//
1304// Matrix inequality
1305//
1306template<typename T, Index M, Index N>
1308bool
1310{
1311 return not_equal(A, B);
1312}
1313
1314//
1315// Scalar matrix product
1316//
1317template<typename S, typename T, Index M, Index N>
1319typename
1320lazy_disable_if<order_1234<S>, apply_matrix<Promote<S, T>, M, N>>::type
1321operator*(S const & s, Matrix<T, M, N> const & A)
1322{
1324 B(A.get_num_rows(), A.get_num_cols());
1325
1326 scale(A, s, B);
1327
1328 return B;
1329}
1330
1331//
1332// Matrix scalar product
1333//
1334template<typename S, typename T, Index M, Index N>
1336typename
1337lazy_disable_if<order_1234<S>, apply_matrix<Promote<S, T>, M, N>>::type
1338operator*(Matrix<T, M, N> const & A, S const & s)
1339{
1341 B(A.get_num_rows(), A.get_num_cols());
1342
1343 scale(A, s, B);
1344
1345 return B;
1346}
1347
1348//
1349// Matrix scalar division
1350//
1351template<typename S, typename T, Index M, Index N>
1353Matrix<typename Promote<S, T>::type, M, N>
1354operator/(Matrix<T, M, N> const & A, S const & s)
1355{
1357 B(A.get_num_rows(), A.get_num_cols());
1358
1359 divide(A, s, B);
1360
1361 return B;
1362}
1363
1364//
1365// Scalar matrix division
1366//
1367template<typename S, typename T, Index M, Index N>
1369Matrix<typename Promote<S, T>::type, M, N>
1370operator/(S const & s, Matrix<T, M, N> const & A)
1371{
1373 B(A.get_num_rows(), A.get_num_cols());
1374
1375 split(A, s, B);
1376
1377 return B;
1378}
1379
1380//
1381// Matrix vector product v = A u
1382//
1383template<typename S, typename T, Index M, Index N>
1385Vector<typename Promote<S, T>::type, M>
1387{
1388 return dot(A, u);
1389}
1390
1391//
1392// Vector matrix product v = u A
1393//
1394template<typename S, typename T, Index M, Index N>
1396Vector<typename Promote<S, T>::type, N>
1398{
1399 return dot(u, A);
1400}
1401
1402//
1403// Matrix tensor product C = A B
1404//
1405template<typename S, typename T, Index M, Index N>
1407Matrix<typename Promote<S, T>::type, M, N>
1409{
1410 return dot(A, B);
1411}
1412
1413//
1414// Tensor matrix product C = A B
1415//
1416template<typename S, typename T, Index M, Index N>
1418Matrix<typename Promote<S, T>::type, M, N>
1420{
1421 return dot(A, B);
1422}
1423
1424//
1425// Matrix dot product C = A B
1426//
1427template<typename S, typename T, Index M, Index P, Index N>
1429Matrix<typename Promote<S, T>::type, M, N>
1431{
1432 return dot(A, B);
1433}
1434
1435//
1436// Matrix vector product v = A u
1437//
1438template<typename S, typename T, Index M, Index N>
1440Vector<typename Promote<S, T>::type, M>
1441dot(Matrix<T, M, N> const & A, Vector<S, N> const & u)
1442{
1443 Index const
1444 num_rows = A.get_num_rows();
1445
1446 Index const
1447 num_cols = A.get_num_cols();
1448
1449 assert(u.get_dimension() == num_cols);
1450
1452 v(num_rows);
1453
1454 for (Index i = 0; i < num_rows; ++i) {
1455
1456 typename Promote<S, T>::type
1457 s = 0.0;
1458
1459 for (Index p = 0; p < num_cols; ++p) {
1460 s += A(i, p) * u(p);
1461 }
1462 v(i) = s;
1463 }
1464
1465 return v;
1466}
1467
1468//
1469// Matrix vector product v = A^T u
1470//
1471template<typename S, typename T, Index M, Index N>
1473Vector<typename Promote<S, T>::type, N>
1474dot(Vector<S, M> const & u, Matrix<T, M, N> const & A)
1475{
1476 Index const
1477 num_rows = A.get_num_rows();
1478
1479 Index const
1480 num_cols = A.get_num_cols();
1481
1482 assert(u.get_dimension() == num_rows);
1483
1485 v(num_cols);
1486
1487 for (Index i = 0; i < num_cols; ++i) {
1488
1489 typename Promote<S, T>::type
1490 s = 0.0;
1491
1492 for (Index p = 0; p < num_rows; ++p) {
1493 s += A(p, i) * u(p);
1494 }
1495 v(i) = s;
1496 }
1497
1498 return v;
1499}
1500
1501//
1502// Matrix tensor product C = A B
1503//
1504template<typename S, typename T, Index M, Index N>
1506Matrix<typename Promote<S, T>::type, M, N>
1507dot(Matrix<T, M, N> const & A, Tensor<S, N> const & B)
1508{
1509 Index const
1510 num_rows = A.get_num_rows();
1511
1512 Index const
1513 num_cols = A.get_num_cols();
1514
1515 assert(B.get_dimension() == num_cols);
1516
1518 C(num_rows, num_cols);
1519
1520 for (Index i = 0; i < num_rows; ++i) {
1521 for (Index j = 0; j < num_cols; ++j) {
1522
1523 typename Promote<S, T>::type
1524 s = 0.0;
1525
1526 for (Index p = 0; p < num_cols; ++p) {
1527 s += A(i, p) * B(p, j);
1528 }
1529 C(i, j) = s;
1530 }
1531 }
1532
1533 return C;
1534}
1535
1536//
1537// Tensor matrix product C = A B
1538//
1539template<typename S, typename T, Index M, Index N>
1541Matrix<typename Promote<S, T>::type, M, N>
1542dot(Tensor<S, M> const & A, Matrix<T, M, N> const & B)
1543{
1544 Index const
1545 num_rows = B.get_num_rows();
1546
1547 Index const
1548 num_cols = B.get_num_cols();
1549
1550 assert(A.get_dimension() == num_rows);
1551
1553 C(num_rows, num_cols);
1554
1555 for (Index i = 0; i < num_rows; ++i) {
1556 for (Index j = 0; j < num_cols; ++j) {
1557
1558 typename Promote<S, T>::type
1559 s = 0.0;
1560
1561 for (Index p = 0; p < num_rows; ++p) {
1562 s += A(i, p) * B(p, j);
1563 }
1564 C(i, j) = s;
1565 }
1566 }
1567
1568 return C;
1569}
1570
1571//
1572// Matrix product C = A B
1573//
1574template<typename S, typename T, Index M, Index P, Index N>
1576Matrix<typename Promote<S, T>::type, M, N>
1578{
1579 Index const
1580 num_rows = A.get_num_rows();
1581
1582 Index const
1583 num_cols = B.get_num_cols();
1584
1585 Index const
1586 dimension = A.get_num_cols();
1587
1588 assert(B.get_num_rows() == dimension);
1589
1591 C(num_rows, num_cols);
1592
1593 for (Index i = 0; i < num_rows; ++i) {
1594 for (Index j = 0; j < num_cols; ++j) {
1595
1596 typename Promote<S, T>::type
1597 s = 0.0;
1598
1599 for (Index p = 0; p < dimension; ++p) {
1600 s += A(i, p) * B(p, j);
1601 }
1602 C(i, j) = s;
1603 }
1604 }
1605
1606 return C;
1607}
1608
1609//
1610// Matrix product C = A^T B
1611//
1612template<typename S, typename T, Index M, Index P, Index N>
1614Matrix<typename Promote<S, T>::type, M, N>
1616{
1617 Index const
1618 num_rows = A.get_num_cols();
1619
1620 Index const
1621 num_cols = B.get_num_cols();
1622
1623 Index const
1624 dimension = A.get_num_rows();
1625
1626 assert(B.get_num_rows() == dimension);
1627
1629 C(num_rows, num_cols);
1630
1631 for (Index i = 0; i < num_rows; ++i) {
1632 for (Index j = 0; j < num_cols; ++j) {
1633
1634 typename Promote<S, T>::type
1635 s = 0.0;
1636
1637 for (Index p = 0; p < dimension; ++p) {
1638 s += A(p, i) * B(p, j);
1639 }
1640 C(i, j) = s;
1641 }
1642 }
1643
1644 return C;
1645}
1646
1647//
1648// Matrix product C = A B^T
1649//
1650template<typename S, typename T, Index M, Index P, Index N>
1652Matrix<typename Promote<S, T>::type, M, N>
1654{
1655 Index const
1656 num_rows = A.get_num_rows();
1657
1658 Index const
1659 num_cols = B.get_num_rows();
1660
1661 Index const
1662 dimension = A.get_num_cols();
1663
1664 assert(B.get_num_cols() == dimension);
1665
1667 C(num_rows, num_cols);
1668
1669 for (Index i = 0; i < num_rows; ++i) {
1670 for (Index j = 0; j < num_cols; ++j) {
1671
1672 typename Promote<S, T>::type
1673 s = 0.0;
1674
1675 for (Index p = 0; p < dimension; ++p) {
1676 s += A(i, p) * B(j, p);
1677 }
1678 C(i, j) = s;
1679 }
1680 }
1681
1682 return C;
1683}
1684
1685//
1686// Matrix product C = A^T B^T
1687// \param A tensor
1688// \param B tensor
1689// \return a tensor \f$ A^T \cdot B^T \f$
1690//
1691template<typename S, typename T, Index M, Index P, Index N>
1693Matrix<typename Promote<S, T>::type, M, N>
1695{
1696 Index const
1697 num_rows = A.get_num_cols();
1698
1699 Index const
1700 num_cols = B.get_num_rows();
1701
1702 Index const
1703 dimension = A.get_num_rows();
1704
1705 assert(B.get_num_cols() == dimension);
1706
1708 C(num_rows, num_cols);
1709
1710 for (Index i = 0; i < num_rows; ++i) {
1711 for (Index j = 0; j < num_cols; ++j) {
1712
1713 typename Promote<S, T>::type
1714 s = 0.0;
1715
1716 for (Index p = 0; p < dimension; ++p) {
1717 s += A(p, i) * B(j, p);
1718 }
1719 C(i, j) = s;
1720 }
1721 }
1722
1723 return C;
1724}
1725
1726//
1727// Matrix double dot product (contraction)
1728//
1729template<typename S, typename T, Index M, Index N>
1731typename Promote<S, T>::type
1733{
1734 Index const
1735 num_rows = A.get_num_rows();
1736
1737 Index const
1738 num_cols = B.get_num_cols();
1739
1740 assert(B.get_num_rows() == num_rows);
1741 assert(B.get_num_cols() == num_cols);
1742
1743 typename Promote<S, T>::type
1744 s = 0.0;
1745
1746 for (Index p = 0; p < num_rows; ++p) {
1747 for (Index q = 0; q < num_cols; ++q) {
1748 s += A(p, q) * B(p, q);
1749 }
1750 }
1751
1752 return s;
1753}
1754
1755//
1756// exterior matrix product
1757//
1758template<typename S, typename T, Index M, Index N>
1760Matrix<typename Promote<S, T>::type, M, N>
1761matrix(Vector<S, M> const & u, Vector<T, N> const & v)
1762{
1763 Index const
1764 num_rows = u.get_dimension();
1765
1766 Index const
1767 num_cols = v.get_dimension();
1768
1770 A(num_rows, num_cols);
1771
1772 for (Index i = 0; i < num_rows; ++i) {
1773
1774 typename Promote<S, T>::type const
1775 s = u(i);
1776
1777 for (Index j = 0; j < num_cols; ++j) {
1778 A(i, j) = s * v(j);
1779 }
1780 }
1781
1782 return A;
1783}
1784
1785//
1786// matrix transpose
1787//
1788template<typename T, Index M, Index N>
1790Matrix<T, M, N>
1792{
1793 Index const
1794 num_rows = A.get_num_cols();
1795
1796 Index const
1797 num_cols = A.get_num_rows();
1798
1800 B(num_rows, num_cols);
1801
1802 for (Index i = 0; i < num_rows; ++i) {
1803 for (Index j = i + 1; j < num_cols; ++j) {
1804 B(i, j) = A(j, i);
1805 }
1806 }
1807
1808 return B;
1809}
1810
1811} // namespace minitensor
1812namespace minitensor {
1813
1814//
1815// Matrix input
1816//
1817template<typename T, Index M, Index N>
1818std::istream &
1819operator>>(std::istream & is, Matrix<T, M, N> & A)
1820{
1821 Index const
1822 num_rows = A.get_num_rows();
1823
1824 Index const
1825 num_cols = A.get_num_cols();
1826
1827 for (Index i = 0; i < num_rows; ++i) {
1828 for (Index j = 0; j < num_cols; ++j) {
1829 is >> A(i,j);
1830 }
1831 }
1832
1833 return is;
1834}
1835
1836//
1837// Matrix output
1838//
1839template<typename T, Index M, Index N>
1840std::ostream &
1841operator<<(std::ostream & os, Matrix<T, M, N> const & A)
1842{
1843 Index const
1844 num_rows = A.get_num_rows();
1845
1846 Index const
1847 num_cols = A.get_num_cols();
1848
1849 Index const
1850 dimension = num_rows * num_cols;
1851
1852 if (dimension == 0) {
1853 return os;
1854 }
1855
1856 os << std::scientific << std::setprecision(17);
1857
1858 for (Index i = 0; i < num_rows; ++i) {
1859
1860 os << std::setw(24) << A(i,0);
1861
1862 for (Index j = 1; j < num_cols; ++j) {
1863 os << "," << std::setw(24) << A(i,j);
1864 }
1865
1866 os << std::endl;
1867 }
1868
1869 return os;
1870}
1871
1873} // namespace minitensor
1874
1875#endif //MiniTensor_Matrix_h
#define KOKKOS_INLINE_FUNCTION
static KOKKOS_INLINE_FUNCTION constexpr Index get_order()
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > operator+(Matrix< S, M, N > const &A, Matrix< T, M, N > const &B)
KOKKOS_INLINE_FUNCTION void minus(TensorBase< T, ST > const &A, TensorBase< T, ST > &B)
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > operator-(Matrix< S, M, N > const &A, Matrix< T, M, N > const &B)
KOKKOS_INLINE_FUNCTION T const & operator()(Index const i, Index const j) const
std::ostream & operator<<(std::ostream &os, Matrix< T, M, N > const &A)
KOKKOS_INLINE_FUNCTION void add(TensorBase< R, SR > const &A, TensorBase< S, SS > const &B, TensorBase< T, ST > &C)
KOKKOS_INLINE_FUNCTION Matrix()
static constexpr bool IS_DYNAMIC
std::istream & operator>>(std::istream &is, Matrix< T, M, N > &A)
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > operator/(Matrix< T, M, N > const &A, S const &s)
KOKKOS_INLINE_FUNCTION Vector< typename Promote< S, T >::type, M > operator*(Matrix< T, M, N > const &A, Vector< S, N > const &u)
KOKKOS_INLINE_FUNCTION Promote< S, T >::type dotdot(Matrix< S, M, N > const &A, Matrix< T, M, N > const &B)
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > t_dot_t(Matrix< S, P, M > const &A, Matrix< T, N, P > const &B)
KOKKOS_INLINE_FUNCTION Index get_dimension() const
KOKKOS_INLINE_FUNCTION bool equal(TensorBase< T, ST > const &A, TensorBase< T, ST > const &B)
KOKKOS_INLINE_FUNCTION Vector< T, M > col(Matrix< T, M, N > const &A, Index const j)
KOKKOS_INLINE_FUNCTION Index get_dimension() const
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > dot_t(Matrix< S, M, P > const &A, Matrix< T, N, P > const &B)
static constexpr Index ORDER
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > matrix(Vector< S, M > const &u, Vector< T, N > const &v)
KOKKOS_INLINE_FUNCTION void scale(TensorBase< R, SR > const &A, S const &s, TensorBase< T, ST > &B)
KOKKOS_INLINE_FUNCTION bool operator!=(Matrix< T, M, N > const &A, Matrix< T, M, N > const &B)
KOKKOS_INLINE_FUNCTION Matrix< T, M, N > adjoint(Matrix< T, N, M > const &A)
KOKKOS_INLINE_FUNCTION void set_dimensions(Index const rows, Index const cols)
KOKKOS_INLINE_FUNCTION bool operator==(Matrix< T, M, N > const &A, Matrix< T, M, N > const &B)
KOKKOS_INLINE_FUNCTION Vector< typename Promote< S, T >::type, M > dot(Matrix< T, M, N > const &A, Vector< S, N > const &u)
KOKKOS_INLINE_FUNCTION Index get_num_cols() const
KOKKOS_INLINE_FUNCTION void split(TensorBase< R, SR > const &A, S const &s, TensorBase< T, ST > &B)
KOKKOS_INLINE_FUNCTION Index get_num_rows() const
std::pair< Index, Index > get_dimensions() const
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > t_dot(Matrix< S, P, M > const &A, Matrix< T, P, N > const &B)
KOKKOS_INLINE_FUNCTION bool not_equal(TensorBase< T, ST > const &A, TensorBase< T, ST > const &B)
KOKKOS_INLINE_FUNCTION void set_dimension(Index const dimension, Index const order)
KOKKOS_INLINE_FUNCTION Matrix< T, M, N > transpose(Matrix< T, N, M > const &A)
virtual KOKKOS_INLINE_FUNCTION ~Matrix()
KOKKOS_INLINE_FUNCTION void divide(TensorBase< R, SR > const &A, S const &s, TensorBase< T, ST > &B)
KOKKOS_INLINE_FUNCTION Vector< T, N > row(Matrix< T, M, N > const &A, Index const i)
KOKKOS_INLINE_FUNCTION void subtract(TensorBase< R, SR > const &A, TensorBase< S, SS > const &B, TensorBase< T, ST > &C)
uint32_t Index
Indexing type.
constexpr Index DYNAMIC
Indicator for dynamic storage.