Function computes the derivative of Kronecker product A*B
    (or Kronecker sum A+B).

    Input:
    -----------------------

    A: 2D matrix
        Some matrix
    dA: 3D (or 2D matrix)
        Derivarives of A
    B: 2D matrix
        Some matrix
    dB: 3D (or 2D matrix)
        Derivarives of B

    operation: str 'prod' or 'sum'
        Which operation is considered. If the operation is 'sum' it is assumed
        that A and are square matrices.s

    Output:
        dC: 3D matrix
        Derivative of Kronecker product A*B (or Kronecker sum A+B)

def dkron(A,dA,B,dB, operation='prod'):
    """
    Function computes the derivative of Kronecker product A*B
    (or Kronecker sum A+B).

    Input:
    -----------------------

    A: 2D matrix
        Some matrix
    dA: 3D (or 2D matrix)
        Derivarives of A
    B: 2D matrix
        Some matrix
    dB: 3D (or 2D matrix)
        Derivarives of B

    operation: str 'prod' or 'sum'
        Which operation is considered. If the operation is 'sum' it is assumed
        that A and are square matrices.s

    Output:
        dC: 3D matrix
        Derivative of Kronecker product A*B (or Kronecker sum A+B)
    """

    if dA is None:
        dA_param_num = 0
        dA = np.zeros((A.shape[0], A.shape[1],1))
    else:
        dA_param_num = dA.shape[2]

    if dB is None:
        dB_param_num = 0
        dB = np.zeros((B.shape[0], B.shape[1],1))
    else:
        dB_param_num = dB.shape[2]

    # Space allocation for derivative matrix
    dC = np.zeros((A.shape[0]*B.shape[0], A.shape[1]*B.shape[1], dA_param_num +  dB_param_num))

    for k in range(dA_param_num):
        if operation == 'prod':
            dC[:,:,k] = np.kron(dA[:,:,k],B);
        else:
            dC[:,:,k] = np.kron(dA[:,:,k],np.eye( B.shape[0] ))

    for k in range(dB_param_num):
        if operation == 'prod':
            dC[:,:,dA_param_num+k] = np.kron(A,dB[:,:,k])
        else:
            dC[:,:,dA_param_num+k] = np.kron(np.eye( A.shape[0] ),dB[:,:,k])

    return dC