xspacelib.h

Go to the documentation of this file.

#ifndef _DG_XSPACELIB_CUH_
#define _DG_XSPACELIB_CUH_
 
#include <thrust/host_vector.h>
 
#include <cusp/coo_matrix.h>
#include <cassert>
 
#include "dg/backend/typedefs.h"
#include "grid.h"
#include "dlt.h"
#include "operator.h"
#include "operator_tensor.h"
 
 
namespace dg{
 
template< class T>
cusp::csr_matrix< int, T, cusp::host_memory> tensorproduct(
        const cusp::csr_matrix< int, T, cusp::host_memory>& lhs,
        const cusp::csr_matrix< int, T, cusp::host_memory>& rhs)
{
    //dimensions of the matrix
    int num_rows     = lhs.num_rows*rhs.num_rows;
    int num_cols     = lhs.num_cols*rhs.num_cols;
    int num_entries  = lhs.num_entries* rhs.num_entries;
    // allocate output matrix
    cusp::csr_matrix<int, T, cusp::host_memory> A(num_rows, num_cols, num_entries);
    //LHS x RHS
    A.row_offsets[0] = 0;
    int counter = 0;
    for( unsigned i=0; i<lhs.num_rows; i++)
    for( unsigned j=0; j<rhs.num_rows; j++)
    {
        int num_entries_in_row =
            (lhs.row_offsets[i+1] - lhs.row_offsets[i])*
            (rhs.row_offsets[j+1] - rhs.row_offsets[j]);
        A.row_offsets[i*rhs.num_rows+j+1] =
            A.row_offsets[i*rhs.num_rows+j] + num_entries_in_row;
        for( int k=lhs.row_offsets[i]; k<lhs.row_offsets[i+1]; k++)
        for( int l=rhs.row_offsets[j]; l<rhs.row_offsets[j+1]; l++)
        {
            A.column_indices[counter] =
                lhs.column_indices[k]*rhs.num_cols +  rhs.column_indices[l];
            A.values[counter]  = lhs.values[k]*rhs.values[l];
            counter++;
        }
    }
    return A;
}
 
namespace create{
 
template<class real_type>
dg::IHMatrix_t<real_type> backscatter( const RealGrid1d<real_type>& g)
{
    //create equidistant backward transformation
    dg::Operator<real_type> backwardeq( g.dlt().backwardEQ());
    dg::Operator<real_type> forward( g.dlt().forward());
    dg::Operator<real_type> backward1d = backwardeq*forward;
 
    return (dg::IHMatrix_t<real_type>)dg::tensorproduct( g.N(), backward1d);
}
 
template<class real_type>
dg::IHMatrix_t<real_type> backscatter( const aRealTopology2d<real_type>& g)
{
    //create equidistant backward transformation
    auto transformX = backscatter( g.gx());
    auto transformY = backscatter( g.gy());
    return dg::tensorproduct( transformY, transformX);
}
 
template<class real_type>
dg::IHMatrix_t<real_type> backscatter( const aRealTopology3d<real_type>& g)
{
    auto transformX = backscatter( g.gx());
    auto transformY = backscatter( g.gy());
    auto transformZ = backscatter( g.gz());
    return dg::tensorproduct( transformZ, dg::tensorproduct(transformY, transformX));
}
 
template<class real_type>
dg::IHMatrix_t<real_type> inv_backscatter( const RealGrid1d<real_type>& g)
{
    //create equidistant backward transformation
    dg::Operator<real_type> backwardeq( g.dlt().backwardEQ());
    dg::Operator<real_type> backward( g.dlt().backward());
    dg::Operator<real_type> forward1d = backward*dg::invert(backwardeq);
 
    return (dg::IHMatrix_t<real_type>)dg::tensorproduct( g.N(), forward1d);
}
template<class real_type>
dg::IHMatrix_t<real_type> inv_backscatter( const aRealTopology2d<real_type>& g)
{
    //create equidistant backward transformation
    auto transformX = inv_backscatter( g.gx());
    auto transformY = inv_backscatter( g.gy());
    return dg::tensorproduct( transformY, transformX);
}
 
template<class real_type>
dg::IHMatrix_t<real_type> inv_backscatter( const aRealTopology3d<real_type>& g)
{
    auto transformX = inv_backscatter( g.gx());
    auto transformY = inv_backscatter( g.gy());
    auto transformZ = inv_backscatter( g.gz());
    return dg::tensorproduct( transformZ, dg::tensorproduct(transformY, transformX));
}
 
 
} //namespace create
} //namespace dg
#endif // _DG_XSPACELIB_CUH_