1 // Copyright (c) 2020 John Abbott, Anna M. Bigatti 2 3 // This file is part of the source of CoCoALib, the CoCoA Library. 4 5 // CoCoALib is free software: you can redistribute it and/or modify 6 // it under the terms of the GNU General Public License as published by 7 // the Free Software Foundation, either version 3 of the License, or 8 // (at your option) any later version. 9 10 // CoCoALib is distributed in the hope that it will be useful, 11 // but WITHOUT ANY WARRANTY; without even the implied warranty of 12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 // GNU General Public License for more details. 14 15 // You should have received a copy of the GNU General Public License 16 // along with CoCoALib. If not, see <http://www.gnu.org/licenses/>. 17 18 #include "CoCoA/MatrixOps.H" 19 20 #include "CoCoA/apply.H" 21 #include "CoCoA/error.H" 22 #include "CoCoA/DenseMatrix.H" 23 #include "CoCoA/MachineInt.H" 24 #include "CoCoA/matrix.H" 25 #include "CoCoA/ring.H" 26 #include "CoCoA/ring-AutomaticConversion.H" 27 #include "CoCoA/RingHom.H" 28 29 #include <limits> 30 using std::numeric_limits; 31 32 namespace CoCoA 33 { 34 35 // Simple "dense" impl KroneckerProd(ConstMatrixView M1,ConstMatrixView M2)36 matrix KroneckerProd(ConstMatrixView M1, ConstMatrixView M2) 37 { 38 const char* const FnName = "KroneckerProd"; 39 CoCoA_STATIC_ERROR_MESG(ErrMixed, ERR::MixedRings,FnName); 40 const ring& R1 = RingOf(M1); 41 const ring& R2 = RingOf(M2); 42 if (R1 != R2) 43 { 44 const RingHom promote = AutomaticConversionHom(R1,R2,ErrMixed); // throws ErrMixed if auto-conv not possible 45 if (codomain(promote) == R1) 46 return KroneckerProd(M1, apply(promote,M2)); 47 return KroneckerProd(apply(promote,M1), M2); 48 } 49 // Here M1 and M2 are over the same ring. 50 const long r1 = NumRows(M1); // int or long??? 51 const long c1 = NumCols(M1); // 52 const long r2 = NumRows(M2); // 53 const long c2 = NumCols(M2); // 54 if (c1 > 1 && c2 > 1 && c1 > numeric_limits<long>::max() / c2) //avoid overflow 55 CoCoA_THROW_ERROR(ERR::BadColIndex, FnName); // ERR::TooBig??? 56 if (r1 > 1 && r2 > 1 && r1 > numeric_limits<long>::max() / r2) //avoid overflow 57 CoCoA_THROW_ERROR(ERR::BadRowIndex, FnName); // ERR::TooBig??? 58 59 const long r = r1*r2; 60 const long c = c1*c2; 61 62 matrix M = NewDenseMat(R1, r,c); 63 64 for (long i1=0; i1 < r1; ++i1) 65 for (long i2=0; i2 < r2; ++i2) 66 for (long j1=0; j1 < c1; ++j1) 67 for (long j2=0; j2 < c2; ++j2) 68 SetEntry(M,i1*r2+i2, j1*c2+j2, M1(i1,j1)*M2(i2,j2)); 69 70 return M; 71 } 72 73 // ORIG DEFN; better local var names? 74 // //--------- tensor matrix 75 76 // matrix TensorMat(ConstMatrixView A, ConstMatrixView B) 77 // { 78 // if (RingOf(A) != RingOf(B)) 79 // CoCoA_THROW_ERROR(ERR::MixedRings, "TensorMat(A,B)"); 80 // if (NumCols(A) > numeric_limits<long>::max() /NumCols(B)) //avoid overflow 81 // CoCoA_THROW_ERROR(ERR::BadColIndex, "TensorMat(A,B)"); 82 // if (NumRows(A) > numeric_limits<long>::max() /NumRows(B)) //avoid overflow 83 // CoCoA_THROW_ERROR(ERR::BadRowIndex, "TensorMat(A,B)"); 84 // matrix ans = NewDenseMat(RingOf(A), NumRows(A)*NumRows(B), NumCols(A)*NumCols(B)); 85 // for (long iA=0; iA < NumRows(A); ++iA) 86 // for (long jA=0; jA < NumCols(A); ++jA) 87 // for (long iB=0; iB < NumRows(B); ++iB) 88 // for (long jB=0; jB < NumCols(B); ++jB) 89 // SetEntry(ans, iA*NumRows(B)+iB, jA*NumCols(B)+jB, A(iA, jA)*B(iB, jB)); 90 // return ans; 91 // } 92 93 94 95 } // end of namespace CoCoA 96 97 98 // RCS header/log in the next few lines 99 // $Header: /Volumes/Home_1/cocoa/cvs-repository/CoCoALib-0.99/src/AlgebraicCore/MatrixOps-KroneckerProd.C,v 1.3 2020/06/22 15:43:37 abbott Exp $ 100 // $Log: MatrixOps-KroneckerProd.C,v $ 101 // Revision 1.3 2020/06/22 15:43:37 abbott 102 // Summary: Use new CoCoA_STATIC_ERROR_MESG macro 103 // 104 // Revision 1.2 2020/06/17 15:49:23 abbott 105 // Summary: Changed CoCoA_ERROR into CoCoA_THROW_ERROR 106 // 107 // Revision 1.1 2020/05/26 12:06:18 abbott 108 // Summary: Renamed TensorMat to KroneckerProd; doc & tests updated 109 // 110 // 111 // 112