[ublas] Kronecker tensor product

From: Nico Galoppo (ngaloppo_at_[hidden])
Date: 2008-03-04 18:24:11

Hi,

I'm currently computing the Kronecker tensor product of two matrices A
and B, both r-by-r to compute the Kronecker tensor product, i.e. in
matrix parlance:

C = A(:) * B(:)';

I'm currently doing this with BLAS and ublas in a way that I think is
pretty dangerous/hacky, but it was the only way I could come up with.
If anyone out there has better options, please let me know.

Note: I'm using Gunther's storage_adaptors.hpp to flatten the matrices
into a vector.

        typedef ublas::matrix<real_type, ublas::column_major> matrix_type;
        size_t r2 = r*r;
        matrix_type A(r,r), B(r,r), C(r2,r2);

        ....

        // tensor product C += A(:) * B(:)'
        // FIXME: this is ugly! (assumes column_major)
        ublas::vector<real_type> A_flat(
ublas::make_vector_from_pointer( r2, &(qiqj.data()[0]) ) );
        ublas::vector<real_type> B_flat(
ublas::make_vector_from_pointer(r2, &(qk.data()[0]) ) );

        atlas::gemm(CblasNoTrans, CblasTrans,
            1.0, A_flat, B_flat,
            1.0, C);

-- 
Nico Galoppo :: http://www.ngaloppo.org