Eigen writing a function that takes all complex matrices as inputs?

I'm trying to write a function for the Kronecker product which takes all the various types of complex matrices (to start), e.g. I can put in a MatrixXcf and a Matrix2cf and it will output the Kronecker product of the two as a MatrixXcf, but I'm running up against a problem. As it stands I have

#pragma once
#include <Eigen/Dense>
#include <iostream>
#include <cmath>

Eigen::MatrixXcf otimes(Eigen::MatrixXcf& A, Eigen::MatrixXcf& B) {
    //Get Number of Columns
    int C_col = A.cols() * B.cols();
    int C_row = A.rows() * B.rows();
    //Define new matrix
    Eigen::MatrixXcf C(C_row, C_col);
    //Populate C
    Eigen::Index i_A = 0;
    Eigen::Index i_B = 0;
    Eigen::Index j_A = 0;
    Eigen::Index j_B = 0;
    for (int i = 0; i < C_row; i++) {
        for (int j = 0; j < C_col; j++) {
            i_A = Eigen::Index(floor(i / B.cols()));
            j_A = Eigen::Index(floor(j / B.cols()));
            i_B = Eigen::Index(i % A.cols());
            j_B = Eigen::Index(j % A.cols());
            C(i,j) = A(i_A,j_A)* B(i_B,j_B);
        }
    }
    return C;
}

Which I've tested with two Eigen::MatrixXcfs and gives the expected result.

However, whenever I try to put in some other types of matrices, e.g. one 'Matrix2cf' I get compile errors. For example,

#pragma once
#include <Eigen/Dense>
#include <complex>

using namespace std::complex_literals;

namespace Pauli {
    Eigen::Matrix2cf Px{
        {0, 1},
        {1, 0}
    };

    Eigen::Matrix2cf Py{
        {0, -1.0if},
        {1.0if, 0}
    };

    Eigen::Matrix2cf Pz{
        {1, 0},
        {0, -1}
    };
}

otimes(Pauli::Pz, Pauli::Pz)

This throws the error

Severity    Code    Description Project File    Line    Suppression State   Details
Error (active)  E0434   a reference of type "Eigen::MatrixXcf &" (aka "Eigen::Matrix<std::complex<float>, -1, -1, 0, -1, -1> &") (not const-qualified) cannot be initialized with a value of type "Eigen::Matrix2cf" (aka "Eigen::Matrix<std::complex<float>, 2, 2, 0, 2, 2>")  DMRG.exe - x64 Debug

for the line where I try to otimes two matrices of fixed dimension 2. I'd like to avoid just switching them to MatrixXd , since I know their size is essentially fixed to be a 2 by 2 array, and I'd like to avoid additional memory from using dynamic matrices.

Additionally, trying to pass through an identity matrix

int n = 2;
std::cout << otimes(Pauli::Pz, Eigen::Matrix2cf::Identity(n,n)) << std::endl;

gives this error

Severity    Code    Description Project File    Line    Suppression State   Details
Error (active)  E0304   no instance of overloaded function "Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::Identity [with _Scalar=std::complex<float>, _Rows=2, _Cols=2, _Options=0, _MaxRows=2, _MaxCols=2]" matches the argument list    DMRG.exe - x64 Debug

I'm also not sure what to make of it.

Would I have to overload the function for each possible combination of two complex matrix types? Also how can I implement taking the Identity matrix as an input to a function, since it seems to behave quite differently but I'd like to keep my code as general as possible. Thanks!

Edit: I've changed the otimes function to

    template <typename Derived>
Eigen::MatrixXcf otimes(Eigen::EigenBase<Derived>& A, Eigen::EigenBase<Derived>& B) {
    //Get Number of Columns
    int C_col = A.cols() * B.cols();
    int C_row = A.rows() * B.rows();
    //Define new matrix
    Eigen::MatrixXcf C(C_row, C_col);
    //Populate C
    Eigen::Index i_A = 0;
    Eigen::Index i_B = 0;
    Eigen::Index j_A = 0;
    Eigen::Index j_B = 0;
    for (int i = 0; i < C_row; i++) {
        for (int j = 0; j < C_col; j++) {
            i_A = Eigen::Index(floor(i / B.cols()));
            j_A = Eigen::Index(floor(j / B.cols()));
            i_B = Eigen::Index(i % A.cols());
            j_B = Eigen::Index(j % A.cols());
            C(i,j) = A(i_A,j_A) * B(i_B,j_B);
        }
    }
    return static_cast<Eigen::MatrixXcf> (C);
}

This throws the error Severity Code Description Project File Line Suppression State Details Error C2064 term does not evaluate to a function taking 2 arguments This error is from the line C(i,j) = A(i_A, j_A) * B(i_B, j_B). I think I need a way to access the elements of the matrix. I've also tried using Eigen::Ref<Eigen::MatrixXcf> MatA = A; to let it know it's a matrix, but that gives the error Severity Code Description Project File Line Suppression State Details Error C2440 'initializing': cannot convert from 'Eigen::EigenBase<Derived>' to 'Eigen::Ref<Eigen::Matrix<std::complex<float>,-1,-1,0,-1,-1>,0,Eigen::OuterStride<-1>>'.
I've also tried having the function return just an Eigen::EigenBase<Derived>, but I'm not able to std::cout the function to check the calculation is correct.

Solution

I'll limit myself to a basic fix, nothing fancy. Most importantly, I don't try to make an efficient Kronecker product implementation.

Relevant fixes:

There should be separate template parameters Derived1 and Derived2, so that both matrices can be separate types, e.g. one MatrixXcf and one Matrix2cf
The argument type should be MatrixBase<Derived> to ensure that matrix semantics apply. Technically, DenseBase<> is enough since that specifies the matrix(row, col) index operator that the function uses but for clarity and future extension that may use the matrix products for its implementation, MatrixBase is a better fit
Further, they should be passed as const reference so that it can be used with unnamed and constant input parameters, e.g. otimes(a * 2.f, b + c)
All integer types should be Eigen::Index, not a mix between int and Index. Personally I declare using long_t = Eigen::Index; because I don't like writing Eigen::Index all the time, but that's just me
That whole floor(x / y) stuff is superfluous. Integer division in C/C++ always rounds down. You just convert to floating point without reason, then back to Eigen::Index
The order of the loop counters should be reversed so that the row is the inner loop. Eigen uses column-major storage order. Therefore adjacent row indices are adjacent in memory
I have no idea what the idea was behind that static_cast in the return statement. Don't do that

template <typename Derived1, typename Derived2>
Eigen::MatrixXcf otimes(const Eigen::MatrixBase<Derived1>& A,
                        const Eigen::MatrixBase<Derived2>& B) {
    //Get Number of Columns
    Eigen::Index C_col = A.cols() * B.cols();
    Eigen::Index C_row = A.rows() * B.rows();
    //Define new matrix
    Eigen::MatrixXcf C(C_row, C_col);
    //Populate C
    for (Eigen::Index  j = 0; j < C_col; j++) {
        for (Eigen::Index i = 0; i < C_row; i++) {
            Eigen::Index i_A = i / B.cols();
            Eigen::Index j_A = j / B.cols();
            Eigen::Index i_B = i % A.cols();
            Eigen::Index j_B = j % A.cols();
            C(i,j) = A(i_A,j_A) * B(i_B,j_B);
        }
    }
    return C;
}

Just so you know, there is an unsupported Kronecker product module as part of Eigen. I have no idea how usable it is.

Fixing the actual product

There is a bug remaining in indexing. Just going off the Wikipedia definition, the loop body should be

            Eigen::Index i_A = i / B.rows();
            Eigen::Index j_A = j / B.cols();
            Eigen::Index i_B = i % B.rows();
            Eigen::Index j_B = j % B.cols();
            C(i,j) = A(i_A,j_A) * B(i_B,j_B);

You can make your code shorter and more efficient (eliminating the divisions, allowing vectorization) by writing something like this:

template <typename Derived1, typename Derived2>
Eigen::MatrixXcf otimes(const Eigen::MatrixBase<Derived1>& A,
                        const Eigen::MatrixBase<Derived2>& B) {
    Eigen::Index A_rows = A.rows(), B_rows = B.rows();
    Eigen::Index A_cols = A.cols(), B_cols = B.cols();
    Eigen::MatrixXcf C(A_rows * B_rows, A_cols * B_cols);
    for(Eigen::Index A_col = 0; A_col < A_cols; ++A_col)
        for(Eigen::Index A_row = 0; A_row < A_rows; ++A_row)
            C.block(A_row * B_rows, A_col * B_cols, B_rows, B_cols) =
                  A(A_row, A_col) * B;
    return C;
}