All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
common.hpp
Go to the documentation of this file.
#ifndef VIENNACL_LINALG_DETAIL_ILU_COMMON_HPP_
#define VIENNACL_LINALG_DETAIL_ILU_COMMON_HPP_
/* =========================================================================
Copyright (c) 2010-2016, Institute for Microelectronics,
Institute for Analysis and Scientific Computing,
TU Wien.
Portions of this software are copyright by UChicago Argonne, LLC.
-----------------
ViennaCL - The Vienna Computing Library
-----------------
Project Head: Karl Rupp rupp@iue.tuwien.ac.at
(A list of authors and contributors can be found in the manual)
License: MIT (X11), see file LICENSE in the base directory
============================================================================= */
#include <vector>
#include <cmath>
#include <iostream>
#include <map>
#include <list>
#include "viennacl/forwards.h"
#include "viennacl/tools/tools.hpp"
#include "viennacl/backend/memory.hpp"
#include "viennacl/linalg/host_based/common.hpp"
#include "viennacl/linalg/misc_operations.hpp"
namespace viennacl
{
namespace linalg
{
namespace detail
{
//
// Level Scheduling Setup for ILU:
//
template<typename NumericT, unsigned int AlignmentV>
void level_scheduling_setup_impl(viennacl::compressed_matrix<NumericT, AlignmentV> const & LU,
viennacl::vector<NumericT> const & diagonal_LU,
std::list<viennacl::backend::mem_handle> & row_index_arrays,
std::list<viennacl::backend::mem_handle> & row_buffers,
std::list<viennacl::backend::mem_handle> & col_buffers,
std::list<viennacl::backend::mem_handle> & element_buffers,
std::list<vcl_size_t> & row_elimination_num_list,
bool setup_U)
{
NumericT const * diagonal_buf = viennacl::linalg::host_based::detail::extract_raw_pointer<NumericT>(diagonal_LU.handle());
NumericT const * elements = viennacl::linalg::host_based::detail::extract_raw_pointer<NumericT>(LU.handle());
unsigned int const * row_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(LU.handle1());
unsigned int const * col_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(LU.handle2());
//
// Step 1: Determine row elimination order for each row and build up meta information about the number of entries taking part in each elimination step:
//
std::vector<vcl_size_t> row_elimination(LU.size1());
std::map<vcl_size_t, std::map<vcl_size_t, vcl_size_t> > row_entries_per_elimination_step;
vcl_size_t max_elimination_runs = 0;
for (vcl_size_t row2 = 0; row2 < LU.size1(); ++row2)
{
vcl_size_t row = setup_U ? (LU.size1() - row2) - 1 : row2;
vcl_size_t row_begin = row_buffer[row];
vcl_size_t row_end = row_buffer[row+1];
vcl_size_t elimination_index = 0; //Note: first run corresponds to elimination_index = 1 (otherwise, type issues with int <-> unsigned int would arise
for (vcl_size_t i = row_begin; i < row_end; ++i)
{
unsigned int col = col_buffer[i];
if ( (!setup_U && col < row) || (setup_U && col > row) )
{
elimination_index = std::max<vcl_size_t>(elimination_index, row_elimination[col]);
row_entries_per_elimination_step[row_elimination[col]][row] += 1;
}
}
row_elimination[row] = elimination_index + 1;
max_elimination_runs = std::max<vcl_size_t>(max_elimination_runs, elimination_index + 1);
}
//std::cout << "Number of elimination runs: " << max_elimination_runs << std::endl;
//
// Step 2: Build row-major elimination matrix for each elimination step
//
//std::cout << "Elimination order: " << std::endl;
//for (vcl_size_t i=0; i<row_elimination.size(); ++i)
// std::cout << row_elimination[i] << ", ";
//std::cout << std::endl;
//vcl_size_t summed_rows = 0;
for (vcl_size_t elimination_run = 1; elimination_run <= max_elimination_runs; ++elimination_run)
{
std::map<vcl_size_t, vcl_size_t> const & current_elimination_info = row_entries_per_elimination_step[elimination_run];
// count cols and entries handled in this elimination step
vcl_size_t num_tainted_cols = current_elimination_info.size();
vcl_size_t num_entries = 0;
for (std::map<vcl_size_t, vcl_size_t>::const_iterator it = current_elimination_info.begin();
it != current_elimination_info.end();
++it)
num_entries += it->second;
//std::cout << "num_entries: " << num_entries << std::endl;
//std::cout << "num_tainted_cols: " << num_tainted_cols << std::endl;
if (num_tainted_cols > 0)
{
row_index_arrays.push_back(viennacl::backend::mem_handle());
viennacl::backend::switch_memory_context<unsigned int>(row_index_arrays.back(), viennacl::traits::context(LU));
viennacl::backend::typesafe_host_array<unsigned int> elim_row_index_array(row_index_arrays.back(), num_tainted_cols);
row_buffers.push_back(viennacl::backend::mem_handle());
viennacl::backend::switch_memory_context<unsigned int>(row_buffers.back(), viennacl::traits::context(LU));
viennacl::backend::typesafe_host_array<unsigned int> elim_row_buffer(row_buffers.back(), num_tainted_cols + 1);
col_buffers.push_back(viennacl::backend::mem_handle());
viennacl::backend::switch_memory_context<unsigned int>(col_buffers.back(), viennacl::traits::context(LU));
viennacl::backend::typesafe_host_array<unsigned int> elim_col_buffer(col_buffers.back(), num_entries);
element_buffers.push_back(viennacl::backend::mem_handle());
viennacl::backend::switch_memory_context<NumericT>(element_buffers.back(), viennacl::traits::context(LU));
std::vector<NumericT> elim_elements_buffer(num_entries);
row_elimination_num_list.push_back(num_tainted_cols);
vcl_size_t k=0;
vcl_size_t nnz_index = 0;
elim_row_buffer.set(0, 0);
for (std::map<vcl_size_t, vcl_size_t>::const_iterator it = current_elimination_info.begin();
it != current_elimination_info.end();
++it)
{
//vcl_size_t col = setup_U ? (elimination_matrix.size() - it->first) - 1 : col2;
vcl_size_t row = it->first;
elim_row_index_array.set(k, row);
vcl_size_t row_begin = row_buffer[row];
vcl_size_t row_end = row_buffer[row+1];
for (vcl_size_t i = row_begin; i < row_end; ++i)
{
unsigned int col = col_buffer[i];
if ( (!setup_U && col < row) || (setup_U && col > row) ) //entry of L/U
{
if (row_elimination[col] == elimination_run) // this entry is substituted in this run
{
elim_col_buffer.set(nnz_index, col);
elim_elements_buffer[nnz_index] = setup_U ? elements[i] / diagonal_buf[it->first] : elements[i];
++nnz_index;
}
}
}
elim_row_buffer.set(++k, nnz_index);
}
//
// Wrap in memory_handles:
//
viennacl::backend::memory_create(row_index_arrays.back(), elim_row_index_array.raw_size(), viennacl::traits::context(row_index_arrays.back()), elim_row_index_array.get());
viennacl::backend::memory_create(row_buffers.back(), elim_row_buffer.raw_size(), viennacl::traits::context(row_buffers.back()), elim_row_buffer.get());
viennacl::backend::memory_create(col_buffers.back(), elim_col_buffer.raw_size(), viennacl::traits::context(col_buffers.back()), elim_col_buffer.get());
viennacl::backend::memory_create(element_buffers.back(), sizeof(NumericT) * elim_elements_buffer.size(), viennacl::traits::context(element_buffers.back()), &(elim_elements_buffer[0]));
}
// Print some info:
//std::cout << "Eliminated columns in run " << elimination_run << ": " << num_tainted_cols << " (tainted columns: " << num_tainted_cols << ")" << std::endl;
//summed_rows += eliminated_rows_in_run;
//if (eliminated_rows_in_run == 0)
// break;
}
//std::cout << "Eliminated rows: " << summed_rows << " out of " << row_elimination.size() << std::endl;
}
template<typename NumericT, unsigned int AlignmentV>
void level_scheduling_setup_L(viennacl::compressed_matrix<NumericT, AlignmentV> const & LU,
viennacl::vector<NumericT> const & diagonal_LU,
std::list<viennacl::backend::mem_handle> & row_index_arrays,
std::list<viennacl::backend::mem_handle> & row_buffers,
std::list<viennacl::backend::mem_handle> & col_buffers,
std::list<viennacl::backend::mem_handle> & element_buffers,
std::list<vcl_size_t> & row_elimination_num_list)
{
level_scheduling_setup_impl(LU, diagonal_LU, row_index_arrays, row_buffers, col_buffers, element_buffers, row_elimination_num_list, false);
}
//
// Multifrontal setup of U:
//
template<typename NumericT, unsigned int AlignmentV>
void level_scheduling_setup_U(viennacl::compressed_matrix<NumericT, AlignmentV> const & LU,
viennacl::vector<NumericT> const & diagonal_LU,
std::list<viennacl::backend::mem_handle> & row_index_arrays,
std::list<viennacl::backend::mem_handle> & row_buffers,
std::list<viennacl::backend::mem_handle> & col_buffers,
std::list<viennacl::backend::mem_handle> & element_buffers,
std::list<vcl_size_t> & row_elimination_num_list)
{
level_scheduling_setup_impl(LU, diagonal_LU, row_index_arrays, row_buffers, col_buffers, element_buffers, row_elimination_num_list, true);
}
//
// Multifrontal substitution (both L and U). Will partly be moved to single_threaded/opencl/cuda implementations
//
template<typename NumericT>
void level_scheduling_substitute(viennacl::vector<NumericT> & vec,
std::list<viennacl::backend::mem_handle> const & row_index_arrays,
std::list<viennacl::backend::mem_handle> const & row_buffers,
std::list<viennacl::backend::mem_handle> const & col_buffers,
std::list<viennacl::backend::mem_handle> const & element_buffers,
std::list<vcl_size_t> const & row_elimination_num_list)
{
typedef typename std::list< viennacl::backend::mem_handle >::const_iterator ListIterator;
ListIterator row_index_array_it = row_index_arrays.begin();
ListIterator row_buffers_it = row_buffers.begin();
ListIterator col_buffers_it = col_buffers.begin();
ListIterator element_buffers_it = element_buffers.begin();
typename std::list< vcl_size_t>::const_iterator row_elimination_num_it = row_elimination_num_list.begin();
for (vcl_size_t i=0; i<row_index_arrays.size(); ++i)
{
viennacl::linalg::detail::level_scheduling_substitute(vec, *row_index_array_it, *row_buffers_it, *col_buffers_it, *element_buffers_it, *row_elimination_num_it);
++row_index_array_it;
++row_buffers_it;
++col_buffers_it;
++element_buffers_it;
++row_elimination_num_it;
}
}
} // namespace detail
} // namespace linalg
} // namespace viennacl
#endif
viennacl::backend::typesafe_host_array
Helper class implementing an array on the host. Default case: No conversion necessary.
Definition: util.hpp:92
viennacl::compressed_matrix::size1
const vcl_size_t & size1() const
Returns the number of rows.
Definition: compressed_matrix.hpp:927
tools.hpp
Various little tools used here and there in ViennaCL.
forwards.h
This file provides the forward declarations for the main types used within ViennaCL.
viennacl::linalg::detail::level_scheduling_setup_U
void level_scheduling_setup_U(viennacl::compressed_matrix< NumericT, AlignmentV > const &LU, viennacl::vector< NumericT > const &diagonal_LU, std::list< viennacl::backend::mem_handle > &row_index_arrays, std::list< viennacl::backend::mem_handle > &row_buffers, std::list< viennacl::backend::mem_handle > &col_buffers, std::list< viennacl::backend::mem_handle > &element_buffers, std::list< vcl_size_t > &row_elimination_num_list)
Definition: common.hpp:208
viennacl::linalg::detail::level_scheduling_substitute
void level_scheduling_substitute(viennacl::vector< NumericT > &vec, std::list< viennacl::backend::mem_handle > const &row_index_arrays, std::list< viennacl::backend::mem_handle > const &row_buffers, std::list< viennacl::backend::mem_handle > const &col_buffers, std::list< viennacl::backend::mem_handle > const &element_buffers, std::list< vcl_size_t > const &row_elimination_num_list)
Definition: common.hpp:224
viennacl::compressed_matrix::handle
const handle_type & handle() const
Returns the OpenCL handle to the matrix entry array.
Definition: compressed_matrix.hpp:942
viennacl::compressed_matrix::handle1
const handle_type & handle1() const
Returns the OpenCL handle to the row index array.
Definition: compressed_matrix.hpp:936
NumericT
float NumericT
Definition: bisect.cpp:40
viennacl::linalg::detail::level_scheduling_setup_impl
void level_scheduling_setup_impl(viennacl::compressed_matrix< NumericT, AlignmentV > const &LU, viennacl::vector< NumericT > const &diagonal_LU, std::list< viennacl::backend::mem_handle > &row_index_arrays, std::list< viennacl::backend::mem_handle > &row_buffers, std::list< viennacl::backend::mem_handle > &col_buffers, std::list< viennacl::backend::mem_handle > &element_buffers, std::list< vcl_size_t > &row_elimination_num_list, bool setup_U)
Definition: common.hpp:51
viennacl::compressed_matrix::handle2
const handle_type & handle2() const
Returns the OpenCL handle to the column index array.
Definition: compressed_matrix.hpp:938
viennacl::vcl_size_t
std::size_t vcl_size_t
Definition: forwards.h:75
common.hpp
Common routines for single-threaded or OpenMP-enabled execution on CPU.
viennacl::vector< NumericT >
viennacl::row
vector_expression< const matrix_base< NumericT, F >, const unsigned int, op_row > row(const matrix_base< NumericT, F > &A, unsigned int i)
Definition: matrix.hpp:910
viennacl::traits::context
viennacl::context context(T const &t)
Returns an ID for the currently active memory domain of an object.
Definition: context.hpp:40
viennacl::backend::mem_handle
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:89
viennacl::linalg::detail::level_scheduling_setup_L
void level_scheduling_setup_L(viennacl::compressed_matrix< NumericT, AlignmentV > const &LU, viennacl::vector< NumericT > const &diagonal_LU, std::list< viennacl::backend::mem_handle > &row_index_arrays, std::list< viennacl::backend::mem_handle > &row_buffers, std::list< viennacl::backend::mem_handle > &col_buffers, std::list< viennacl::backend::mem_handle > &element_buffers, std::list< vcl_size_t > &row_elimination_num_list)
Definition: common.hpp:191
viennacl::compressed_matrix
A sparse square matrix in compressed sparse rows format.
Definition: compressed_matrix.hpp:559
viennacl::backend::memory_create
void memory_create(mem_handle &handle, vcl_size_t size_in_bytes, viennacl::context const &ctx, const void *host_ptr=NULL)
Creates an array of the specified size. If the second argument is provided, the buffer is initialized...
Definition: memory.hpp:87
viennacl::vector_base::handle
const handle_type & handle() const
Returns the memory handle.
Definition: vector_def.hpp:128
misc_operations.hpp
Implementations of miscellaneous operations.
memory.hpp
Main interface routines for memory management.