Loading...
Searching...
No Matches
multiplication.tpp
1/*---------------------------------------------------------------------------*\
2 *
3 * bitpit
4 *
5 * Copyright (C) 2015-2021 OPTIMAD engineering Srl
6 *
7 * -------------------------------------------------------------------------
8 * License
9 * This file is part of bitpit.
10 *
11 * bitpit is free software: you can redistribute it and/or modify it
12 * under the terms of the GNU Lesser General Public License v3 (LGPL)
13 * as published by the Free Software Foundation.
14 *
15 * bitpit is distributed in the hope that it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18 * License for more details.
19 *
20 * You should have received a copy of the GNU Lesser General Public License
21 * along with bitpit. If not, see <http://www.gnu.org/licenses/>.
22 *
23\*---------------------------------------------------------------------------*/
24
25namespace bitpit{
26namespace linearalgebra{
31// -------------------------------------------------------------------------- //
39template <class T>
40void matmul(
41 T A,
42 const std::vector< std::vector< T > > &B,
43 std::vector< std::vector< T > > &C
44) {
45
46// ========================================================================== //
47// VARIABLE DECLARATION //
48// ========================================================================== //
49
50// Local variables
51std::size_t m, n;
52
53// ========================================================================== //
54// CHECK INPUT //
55// ========================================================================== //
56
57// Matrix
58m = B.size();
59if (m == 0) {
60 return;
61}
62n = B[0].size();
63if (n == 0) {
64 return;
65}
66
67// ========================================================================== //
68// PERFORM PRODUCT //
69// ========================================================================== //
70
71// Resize output variable
72C.resize(m);
73
74// Perform product
75for (std::size_t i = 0; i < m; i++) {
76 C[i].resize(n, T{});
77 for (std::size_t j = 0; j < n; j++) {
78 C[i][j] = A * B[i][j];
79 } //next j
80} //next i
81
82return; };
83
84// -------------------------------------------------------------------------- //
93template <class T, std::size_t m, std::size_t n>
94void matmul(
95 T A,
96 const std::array< std::array< T, n >, m > &B,
97 std::array< std::array< T, n >, m > &C
98) {
99
100// ========================================================================== //
101// VARIABLE DECLARATION //
102// ========================================================================== //
103
104// Local variables
105// none
106
107// ========================================================================== //
108// CHECK INPUT //
109// ========================================================================== //
110
111// Matrix
112if (m == 0) {
113 return;
114}
115if (n == 0) {
116 return;
117}
118
119// ========================================================================== //
120// PERFORM PRODUCT //
121// ========================================================================== //
122for (std::size_t i = 0; i < m; i++) {
123 for (std::size_t j = 0; j < n; j++) {
124 C[i][j] = A * B[i][j];
125 } //next j
126} //next i
127
128return; };
129
130// -------------------------------------------------------------------------- //
138template <class T>
140 const std::vector< std::vector< T > > &B,
141 T A,
142 std::vector< std::vector< T > > &C
143) {
144
145// ========================================================================== //
146// VARIABLE DECLARATION //
147// ========================================================================== //
148
149// Local variables
150std::size_t m, n;
151
152// Counters
153// none
154
155// ========================================================================== //
156// CHECK INPUT //
157// ========================================================================== //
158
159// Matrix
160m = B.size();
161if (m == 0) {
162 return;
163}
164n = B[0].size();
165if (n == 0) {
166 return;
167}
168
169// ========================================================================== //
170// PERFORM PRODUCT //
171// ========================================================================== //
172matmul(A, B, C);
173
174return; };
175
176// -------------------------------------------------------------------------- //
185template <class T, std::size_t m, std::size_t n>
187 const std::array< std::array< T, n >, m > &B,
188 T A,
189 std::array< std::array< T, n >, m > &C
190) {
191
192// ========================================================================== //
193// VARIABLE DECLARATION //
194// ========================================================================== //
195
196// Local variables
197// none
198
199// Counters
200// none
201
202// ========================================================================== //
203// CHECK INPUT //
204// ========================================================================== //
205
206// Matrix
207if (m == 0) {
208 return;
209}
210if (n == 0) {
211 return;
212}
213
214// ========================================================================== //
215// PERFORM PRODUCT //
216// ========================================================================== //
217matmul(A, B, C);
218
219return; };
220
221// -------------------------------------------------------------------------- //
229template <class T>
231 const std::vector< T > &A,
232 const std::vector< std::vector < T > > &B,
233 std::vector< T > &C
234) {
235
236// ========================================================================== //
237// VARIABLES DECLARATION //
238// ========================================================================== //
239
240// Local variables
241std::size_t l, m, n;
242
243// ========================================================================== //
244// CHECK INPUT //
245// ========================================================================== //
246
247// Input vector
248l = A.size();
249if (l == 0) {
250 return;
251}
252
253// Input matrix
254m = B.size();
255if (m == 0) {
256 return;
257}
258n = B[0].size();
259if (n == 0) {
260 return;
261}
262
263// Check dimensions coherency
264if (l != m) {
265 return;
266}
267
268// ========================================================================== //
269// COMPUTE THE MATRIX PRODUCT //
270// ========================================================================== //
271
272// Resize vector
273C.resize(n, T{});
274
275// Compute matrix product
276for (std::size_t i = 0; i < n; i++) {
277 C[i] = T{};
278 for (std::size_t j = 0; j < m; j++) {
279 C[i] += A[j]*B[j][i];
280 } //next j
281} //next i
282
283return; }
284
285// -------------------------------------------------------------------------- //
294template <class T, std::size_t m, std::size_t n>
296 const std::array< T, m > &A,
297 const std::array< std::array < T, n >, m > &B,
298 std::array< T, n > &C
299) {
300
301// ========================================================================== //
302// VARIABLES DECLARATION //
303// ========================================================================== //
304
305// Local variables
306// none
307
308// ========================================================================== //
309// CHECK INPUT //
310// ========================================================================== //
311
312// Input matrix
313if (m == 0) {
314 return;
315}
316if (n == 0) {
317 return;
318}
319
320// ========================================================================== //
321// COMPUTE THE MATRIX PRODUCT //
322// ========================================================================== //
323for (std::size_t i = 0; i < n; i++) {
324 C[i] = T{};
325 for (std::size_t j = 0; j < m; j++) {
326 C[i] += A[j]*B[j][i];
327 } //next j
328} //next i
329
330return; }
331
332// -------------------------------------------------------------------------- //
340template <class T>
342 const std::vector< std::vector < T > > &A,
343 const std::vector< T > &B,
344 std::vector< T > &C
345) {
346
347// ========================================================================== //
348// VARIABLES DECLARATION //
349// ========================================================================== //
350
351// Local variables
352std::size_t l, m, n;
353
354// ========================================================================== //
355// CHECK INPUT //
356// ========================================================================== //
357
358// Input vector
359l = B.size();
360if (l == 0) {
361 return;
362}
363
364// Input matrix
365m = A.size();
366if (m == 0) {
367 return;
368}
369n = A[0].size();
370if (n == 0) {
371 return;
372}
373
374// Check dimensions coherency
375if (l != n) {
376 return;
377}
378
379// ========================================================================== //
380// COMPUTE THE MATRIX PRODUCT //
381// ========================================================================== //
382
383// Resize vector
384C.resize(m, T{});
385
386// Compute matrix product
387for (std::size_t i = 0; i < m; i++) {
388 C[i] = T{};
389 for (std::size_t j = 0; j < n; j++) {
390 C[i] += B[j]*A[i][j];
391 } //next j
392} //next i
393
394return; }
395
396// -------------------------------------------------------------------------- //
405template <class T, std::size_t m, std::size_t n>
407 const std::array< std::array < T, n >, m > &A,
408 const std::array< T, n > &B,
409 std::array< T, m > &C
410) {
411
412// ========================================================================== //
413// VARIABLES DECLARATION //
414// ========================================================================== //
415
416// Local variables
417// none
418
419// ========================================================================== //
420// CHECK INPUT //
421// ========================================================================== //
422
423// Input matrix
424if (m == 0) {
425 return;
426}
427if (n == 0) {
428 return;
429}
430
431// ========================================================================== //
432// COMPUTE THE MATRIX PRODUCT //
433// ========================================================================== //
434for (std::size_t i = 0; i < m; i++) {
435 C[i] = T{};
436 for (std::size_t j = 0; j < n; j++) {
437 C[i] += B[j]*A[i][j];
438 } //next j
439} //next i
440
441return; }
442
443// -------------------------------------------------------------------------- //
451template <class T>
453 const std::vector< std::vector< T > > &A,
454 const std::vector< std::vector< T > > &B,
455 std::vector< std::vector< T > > &C
456) {
457
458// ========================================================================== //
459// VARIABLE DECLARATION //
460// ========================================================================== //
461
462// Local variables
463std::size_t m1, n1, n2, m2;
464
465// ========================================================================== //
466// CHECK INPUT //
467// ========================================================================== //
468
469// 1st matrix
470m1 = A.size();
471if (m1 == 0) {
472 return;
473}
474n1 = A[0].size();
475if (n1 == 0) {
476 return;
477}
478
479// 2nd matrix
480m2 = B.size();
481if (m2 == 0) {
482 return;
483}
484n2 = B[0].size();
485if (n2 == 0) {
486 return;
487}
488
489// Check dimensions coherency
490if (n1 != m2) {
491 return;
492}
493
494// ========================================================================== //
495// PERFORM PRODUCT //
496// ========================================================================== //
497
498// Resiz output variable
499C.resize(m1);
500
501for (std::size_t i = 0; i < m1; i++) {
502 C[i].resize(n2, T{});
503 for (std::size_t j = 0; j < n2; j++) {
504 C[i][j] = T{};
505 for (std::size_t k = 0; k < n1; k++) {
506 C[i][j] += A[i][k] * B[k][j];
507 } //next k
508 } //next j
509} //next i
510
511return; };
512
513// -------------------------------------------------------------------------- //
522template <class T, std::size_t m, std::size_t n, std::size_t l>
524 const std::array< std::array< T, n >, m > &A,
525 const std::array< std::array< T, l >, n > &B,
526 std::array< std::array< T, l >, m > &C
527) {
528
529// ========================================================================== //
530// VARIABLE DECLARATION //
531// ========================================================================== //
532
533// Local variables
534// none
535
536// ========================================================================== //
537// CHECK INPUT //
538// ========================================================================== //
539
540// 1st matrix
541if (m == 0) {
542 return;
543}
544if (n == 0) {
545 return;
546}
547
548// 2nd matrix
549if (l == 0) {
550 return;
551}
552
553// ========================================================================== //
554// PERFORM PRODUCT //
555// ========================================================================== //
556for (std::size_t i = 0; i < m; i++) {
557 for (std::size_t j = 0; j < l; j++) {
558 C[i][j] = T{};
559 for (std::size_t k = 0; k < n; k++) {
560 C[i][j] += A[i][k] * B[k][j];
561 } //next k
562 } //next j
563} //next i
564
565return; };
566
567// ----------------------------------------------------------------------------------- //
576template <class T>
577std::vector< std::vector<T> > matmul(
578 const std::vector< std::vector<T> > &M,
579 const std::vector<std::vector<T> > &N
580) {
581
582 std::vector< std::vector<T> > Tr = transpose(N);
583
584 std::size_t d1= M.size();
585 std::size_t d2= Tr.size();
586
587 std::vector< std::vector<T> > Q(d1, std::vector<T> (d2, T()) );
588
589 for( std::size_t i=0; i<d1; i++){
590 for( std::size_t j=0; j<d2; j++){
591 Q[i][j]= dotProduct( M[i], Tr[j] );
592 };
593 };
594
595 return (Q);
596};
597
598// ----------------------------------------------------------------------------------- //
606template <class T, std::size_t d1, std::size_t d2, std::size_t d3>
607std::array< std::array<T, d2> , d1> matmul(
608 const std::array< std::array<T, d3>, d1> &M,
609 const std::array<std::array<T, d2>, d3> &N
610){
611 std::array< std::array<T, d2> , d1> Q;
612 std::array< std::array<T, d3> , d2> Tr;
613
614 Tr = transpose( N ) ;
615
616 for( std::size_t i=0; i<d1; i++){
617 for( std::size_t j=0; j<d2; j++){
618 Q[i][j]= dotProduct( M[i], Tr[j] );
619 };
620 };
621
622 return (Q);
623};
624
625// ----------------------------------------------------------------------------------- //
634template <class T>
635std::vector< std::vector<T> > matmulDiag(
636 const std::vector<T> &M,
637 const std::vector<std::vector<T> > &N
638) {
639
640 std::size_t d1= M.size();
641 std::vector< std::vector<T> > Q( N );
642
643 for( std::size_t i=0; i<d1; i++){
644 Q[i] *= M[i] ;
645 };
646
647 return (Q);
648};
649
650// ----------------------------------------------------------------------------------- //
659template <class T>
660std::vector< std::vector<T> > matmulDiag(
661 const std::vector< std::vector<T> > &M,
662 const std::vector<T> &N
663) {
664
665 std::vector< std::vector<T> > Q( M );
666
667 std::size_t d1= M.size() ;
668
669 for( std::size_t i=0; i<d1; i++ ){
670 Q[i] = M[i] * N ;
671 };
672
673 return (Q);
674
675};
676
677// ----------------------------------------------------------------------------------- //
686template <class T, std::size_t d1, std::size_t d2>
687std::array< std::array<T, d2> , d1> matmulDiag(
688 const std::array< T, d1> &M,
689 const std::array<std::array<T, d2>, d1> &N
690){
691
692 std::size_t i;
693 std::array< std::array<T, d2> , d1> Q(N);
694
695 for( i=0; i<d1; i++){
696 Q[i] *= M[i] ;
697 };
698
699 return (Q);
700};
701
702// ----------------------------------------------------------------------------------- //
711template <class T, std::size_t d1, std::size_t d2>
712std::array< std::array<T, d2> , d1> matmulDiag(
713 const std::array<std::array<T, d2>, d1> &M,
714 const std::array< T, d2> &N
715) {
716
717 std::size_t i;
718 std::array< std::array<T, d2> , d1> Q;
719
720 for( i=0; i<d1; i++){
721 Q[i] = M[i] *N ;
722 };
723
724 return (Q);
725};
726
727// Matrix Vector Multiplication ====================================================== //
728
729// ----------------------------------------------------------------------------------- //
738template <class T>
739std::vector<T> matmul(
740 const std::vector< std::vector<T>> &M,
741 const std::vector<T> &x
742) {
743
744 std::size_t d1 = M.size();
745
746 std::vector<T> z(d1, T{});
747
748 for( std::size_t i=0; i<d1; i++){
749 z[i]= dotProduct( M[i], x );
750 }
751
752 return (z);
753};
754
755// ----------------------------------------------------------------------------------- //
764template <class T, std::size_t d1, std::size_t d2>
765std::array<T, d1> matmul(
766 const std::array< std::array<T, d2>, d1> &M,
767 const std::array<T, d2> &x
768) {
769
770 std::array<T, d1> z;
771
772 for( std::size_t i=0; i<d1; i++){
773 z[i]= dotProduct( M[i], x);
774 }
775
776 return (z);
777};
778
779// ----------------------------------------------------------------------------------- //
788template <class T>
789std::vector<std::vector<T>> tensorProduct(
790 const std::vector<T> &x,
791 const std::vector<T> &y
792) {
793
794 std::size_t n = x.size();
795 std::size_t m = y.size();
796 std::vector<T> row(m, T{});
797 std::vector<std::vector<T>> z(n,row) ;
798
799 for( std::size_t i=0; i<n; i++){
800 for( std::size_t j=0; j<m; j++){
801 z[i][j] = x[i] *y[j] ;
802 };
803 };
804
805return (z);}
806
807// ----------------------------------------------------------------------------------- //
816template <class T, std::size_t n, std::size_t m>
817std::array<std::array<T,m>,n> tensorProduct(
818 const std::array<T,n> &x,
819 const std::array<T,m> &y
820) {
821 std::array<std::array<T,m>,n> z ;
822
823 for( std::size_t i=0; i<n; i++){
824 for( std::size_t j=0; j<m; j++){
825 z[i][j] = x[i] *y[j] ;
826 };
827 };
828
829 return (z);
830}
831
835}
836}
T dotProduct(const std::array< T, d > &x, const std::array< T, d > &y)
void transpose(std::vector< std::vector< T > > &, std::vector< std::vector< T > > &)
std::vector< std::vector< T > > tensorProduct(std::vector< T > const &, std::vector< T > const &)
void matmul(T, const std::vector< std::vector< T > > &, std::vector< std::vector< T > > &)
std::vector< std::vector< T > > matmulDiag(const std::vector< T > &, const std::vector< std::vector< T > > &)
--- layout: doxygen_footer ---