operator_tensor.h

Go to the documentation of this file.

#pragma once
 
#include <algorithm>
 
#include "operator.h"
 
namespace dg
{
 
 
 
template< class T>
SquareMatrix<T> tensorproduct( const SquareMatrix<T>& op1, const SquareMatrix<T>& op2)
{
    assert( op1.size() == op2.size());
    unsigned n = op1.size();
    SquareMatrix<T> prod( n*n);
    for( unsigned i=0; i<n; i++)
        for( unsigned j=0; j<n; j++)
            for( unsigned k=0; k<n; k++)
                for( unsigned l=0; l<n; l++)
                    prod(i*n+k, j*n+l) = op1(i,j)*op2(k,l);
    return prod;
}
 
 
template< class T>
dg::SparseMatrix<int,T, thrust::host_vector> tensorproduct( unsigned N, const SquareMatrix<T>& op)
{
    assert( N>0);
    unsigned n = op.size();
    // allocate output matrix
    thrust::host_vector<int> A_row_offsets(n*N+1), A_column_indices( N*n*n);
    thrust::host_vector<T> A_values( N*n*n);
    A_row_offsets[0] = 0;
    for( unsigned k=0; k<N; k++)
        for( unsigned i=0; i<n; i++)
        {
            A_row_offsets[k*n+i+1]      = A_row_offsets[k*n+i] + n;
            for( unsigned j=0; j<n; j++)
                //if( op(i,j) != 0)
                {
                    A_column_indices[(k*n+i)*n+j]   = k*n+j;
                    A_values[(k*n+i)*n+j]           = op(i,j);
                }
        }
    return {n*N, n*N, A_row_offsets, A_column_indices, A_values};
}
 
 
template< class T>
dg::SparseMatrix<int, T, thrust::host_vector> sandwich( const SquareMatrix<T>& left,  const dg::SparseMatrix<int, T, thrust::host_vector>& m, const SquareMatrix<T>& right)
{
    assert( left.size() == right.size());
    unsigned n = left.size();
    unsigned N = m.num_rows()/n;
    return tensorproduct( N, left)*(m*tensorproduct( N, right));
}
 
 
 
 
 
}//namespace dg