From 85132716d335a05bb35f8bb48c65d1c6dce934d3 Mon Sep 17 00:00:00 2001 From: "jan.reininghaus" Date: Mon, 24 Jun 2013 09:48:27 +0000 Subject: new block_spectral_sequence_reduction git-svn-id: https://phat.googlecode.com/svn/trunk@134 8e3bb3c2-eed4-f18f-5264-0b6c94e6926d --- .../algorithms/block_spectral_sequence_reduction.h | 79 ++++++++++++++++++++++ src/benchmark.cpp | 7 +- src/phat.cpp | 7 +- src/self_test.cpp | 14 +++- 4 files changed, 101 insertions(+), 6 deletions(-) create mode 100644 include/phat/algorithms/block_spectral_sequence_reduction.h diff --git a/include/phat/algorithms/block_spectral_sequence_reduction.h b/include/phat/algorithms/block_spectral_sequence_reduction.h new file mode 100644 index 0000000..398236e --- /dev/null +++ b/include/phat/algorithms/block_spectral_sequence_reduction.h @@ -0,0 +1,79 @@ +/* Copyright 2013 IST Austria + Contributed by: Jan Reininghaus + + This file is part of PHAT. + + PHAT is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + PHAT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with PHAT. If not, see . */ + +#pragma once + +#include +#include + +namespace phat { + class block_spectral_sequence_reduction { + public: + template< typename Representation > + void operator () ( boundary_matrix< Representation >& boundary_matrix ) { + + const index nr_columns = boundary_matrix.get_num_cols(); + std::vector< index > lowest_one_lookup( nr_columns, -1 ); + + //const index num_stripes = sqrt( nr_columns ); + const index num_stripes = omp_get_max_threads(); + + index block_size = ( nr_columns % num_stripes == 0 ) ? nr_columns / num_stripes : block_size = nr_columns / num_stripes + 1; + + std::vector< std::vector< index > > unreduced_cols_cur_pass( num_stripes ); + std::vector< std::vector< index > > unreduced_cols_next_pass( num_stripes ); + + for( index cur_dim = boundary_matrix.get_max_dim(); cur_dim >= 1 ; cur_dim-- ) { + #pragma omp parallel for schedule( guided, 1 ) + for( index cur_stripe = 0; cur_stripe < num_stripes; cur_stripe++ ) { + index col_begin = cur_stripe * block_size; + index col_end = std::min( (cur_stripe+1) * block_size, nr_columns ); + for( index cur_col = col_begin; cur_col < col_end; cur_col++ ) + if( boundary_matrix.get_dim( cur_col ) == cur_dim && boundary_matrix.get_max_index( cur_col ) != -1 ) + unreduced_cols_cur_pass[ cur_stripe ].push_back( cur_col ); + } + for( index cur_pass = 0; cur_pass < num_stripes; cur_pass++ ) { + boundary_matrix.sync(); + #pragma omp parallel for schedule( guided, 1 ) + for( int cur_stripe = 0; cur_stripe < num_stripes; cur_stripe++ ) { + index row_begin = (cur_stripe - cur_pass) * block_size; + index row_end = row_begin + block_size; + unreduced_cols_next_pass[ cur_stripe ].clear(); + for( index idx = 0; idx < unreduced_cols_cur_pass[ cur_stripe ].size(); idx++ ) { + index cur_col = unreduced_cols_cur_pass[ cur_stripe ][ idx ]; + index lowest_one = boundary_matrix.get_max_index( cur_col ); + while( lowest_one != -1 && lowest_one >= row_begin && lowest_one < row_end && lowest_one_lookup[ lowest_one ] != -1 ) { + boundary_matrix.add_to( lowest_one_lookup[ lowest_one ], cur_col ); + lowest_one = boundary_matrix.get_max_index( cur_col ); + } + if( lowest_one != -1 ) { + if( lowest_one >= row_begin && lowest_one < row_end ) { + lowest_one_lookup[ lowest_one ] = cur_col; + boundary_matrix.clear( lowest_one ); + } else { + unreduced_cols_next_pass[ cur_stripe ].push_back( cur_col ); + } + } + } + unreduced_cols_next_pass[ cur_stripe ].swap( unreduced_cols_cur_pass[ cur_stripe ] ); + } + } + } + } + }; +} diff --git a/src/benchmark.cpp b/src/benchmark.cpp index 59b5f32..ca2497c 100644 --- a/src/benchmark.cpp +++ b/src/benchmark.cpp @@ -29,6 +29,7 @@ #include #include #include +#include #include @@ -37,7 +38,7 @@ enum Representation_type {VECTOR_VECTOR, VECTOR_SET, SPARSE_PIVOT_COLUMN, FULL_PIVOT_COLUMN, BIT_TREE_PIVOT_COLUMN, VECTOR_LIST}; -enum Algorithm_type {STANDARD, TWIST, ROW, CHUNK, CHUNK_SEQUENTIAL}; +enum Algorithm_type {STANDARD, TWIST, ROW, CHUNK, CHUNK_SEQUENTIAL, BLOCK_SPECTRAL_SEQUENCE}; enum Ansatz_type {PRIMAL, DUAL}; void print_help() { @@ -51,7 +52,7 @@ void print_help() { std::cerr << "--dualize -- use only dualization approach" << std::endl; std::cerr << "--primal -- use only primal approach" << std::endl; std::cerr << "--vector_vector, --vector_set, --vector_list, --full_pivot_column, --sparse_pivot_column, --bit_tree_pivot_column -- use only a subset of representation data structures for boundary matrices" << std::endl; - std::cerr << "--standard, --twist, --chunk, --chunk_sequential, --row -- use only a subset of reduction algorithms" << std::endl; + std::cerr << "--standard, --twist, --chunk, --chunk_sequential, --block_spectral_sequence, --row -- use only a subset of reduction algorithms" << std::endl; } void print_help_and_exit() { @@ -80,6 +81,7 @@ void parse_command_line( int argc, char** argv, bool& use_binary, std::vector< R else if( argument == "--twist" ) algorithms.push_back( TWIST ); else if( argument == "--row" ) algorithms.push_back( ROW ); else if( argument == "--chunk_sequential" ) algorithms.push_back( CHUNK_SEQUENTIAL ); + else if( argument == "--block_spectral_sequence" ) algorithms.push_back( BLOCK_SPECTRAL_SEQUENCE ); else if( argument == "--chunk" ) algorithms.push_back( CHUNK ); else if( argument == "--primal" ) ansaetze.push_back( PRIMAL ); else if( argument == "--dual" ) ansaetze.push_back( DUAL ); @@ -153,6 +155,7 @@ void benchmark( std::string input_filename, bool use_binary, Ansatz_type ansatz case TWIST: std::cout << " twist,"; benchmark< phat::Representation, phat::twist_reduction >( input_filename, use_binary, ansatz ); break; \ case ROW: std::cout << " row,"; benchmark< phat::Representation, phat::row_reduction >( input_filename, use_binary, ansatz ); break; \ case CHUNK: std::cout << " chunk,"; benchmark< phat::Representation, phat::chunk_reduction >( input_filename, use_binary, ansatz ); break; \ + case BLOCK_SPECTRAL_SEQUENCE: std::cout << " block spectral sequence,"; benchmark< phat::Representation, phat::block_spectral_sequence_reduction >( input_filename, use_binary, ansatz ); break; \ case CHUNK_SEQUENTIAL: std::cout << " chunk_sequential,"; \ int num_threads = omp_get_max_threads(); \ omp_set_num_threads( 1 ); \ diff --git a/src/phat.cpp b/src/phat.cpp index c19aea8..1cc5ebb 100644 --- a/src/phat.cpp +++ b/src/phat.cpp @@ -29,11 +29,12 @@ #include #include #include +#include #include enum Representation_type {VECTOR_VECTOR, VECTOR_SET, SPARSE_PIVOT_COLUMN, FULL_PIVOT_COLUMN, BIT_TREE_PIVOT_COLUMN, VECTOR_LIST}; -enum Algorithm_type {STANDARD, TWIST, ROW, CHUNK, CHUNK_SEQUENTIAL }; +enum Algorithm_type {STANDARD, TWIST, ROW, CHUNK, CHUNK_SEQUENTIAL, BLOCK_SPECTRAL_SEQUENCE }; void print_help() { std::cerr << "Usage: " << "phat " << "[options] input_filename output_filename" << std::endl; @@ -46,7 +47,7 @@ void print_help() { std::cerr << "--verbose -- verbose output" << std::endl; std::cerr << "--dualize -- use dualization approach" << std::endl; std::cerr << "--vector_vector, --vector_set, --vector_list, --full_pivot_column, --sparse_pivot_column, --bit_tree_pivot_column -- selects a representation data structure for boundary matrices (default is '--bit_tree_pivot_column')" << std::endl; - std::cerr << "--standard, --twist, --chunk, --chunk_sequential, --row -- selects a reduction algorithm (default is '--twist')" << std::endl; + std::cerr << "--standard, --twist, --chunk, --chunk_sequential, --block_spectral_sequence, --row -- selects a reduction algorithm (default is '--twist')" << std::endl; } void print_help_and_exit() { @@ -79,6 +80,7 @@ void parse_command_line( int argc, char** argv, bool& use_binary, Representation else if( option == "--row" ) algorithm = ROW; else if( option == "--chunk" ) algorithm = CHUNK; else if( option == "--chunk_sequential" ) algorithm = CHUNK_SEQUENTIAL; + else if( option == "--block_spectral_sequence" ) algorithm = BLOCK_SPECTRAL_SEQUENCE; else if( option == "--verbose" ) verbose = true; else if( option == "--help" ) print_help_and_exit(); else print_help_and_exit(); @@ -150,6 +152,7 @@ void compute_pairing( std::string input_filename, std::string output_filename, b case STANDARD: compute_pairing< phat::Representation, phat::standard_reduction> ( input_filename, output_filename, use_binary, verbose, dualize ); break; \ case TWIST: compute_pairing< phat::Representation, phat::twist_reduction> ( input_filename, output_filename, use_binary, verbose, dualize ); break; \ case ROW: compute_pairing< phat::Representation, phat::row_reduction >( input_filename, output_filename, use_binary, verbose, dualize ); break; \ + case BLOCK_SPECTRAL_SEQUENCE: compute_pairing< phat::Representation, phat::block_spectral_sequence_reduction >( input_filename, output_filename, use_binary, verbose, dualize ); break; \ case CHUNK: compute_pairing< phat::Representation, phat::chunk_reduction >( input_filename, output_filename, use_binary, verbose, dualize ); break; \ case CHUNK_SEQUENTIAL: int num_threads = omp_get_max_threads(); \ omp_set_num_threads( 1 ); \ diff --git a/src/self_test.cpp b/src/self_test.cpp index e3aaa32..cfecd0e 100644 --- a/src/self_test.cpp +++ b/src/self_test.cpp @@ -29,6 +29,7 @@ #include #include #include +#include int main( int argc, char** argv ) { @@ -132,6 +133,11 @@ int main( int argc, char** argv ) phat::boundary_matrix< BitTree > row_boundary_matrix = boundary_matrix; phat::compute_persistence_pairs< phat::row_reduction >( row_pairs, row_boundary_matrix ); + std::cout << "Running Block spectral sequence - BitTree ..." << std::endl; + phat::persistence_pairs ss_pairs; + phat::boundary_matrix< BitTree > ss_boundary_matrix = boundary_matrix; + phat::compute_persistence_pairs< phat::block_spectral_sequence_reduction >( ss_pairs, ss_boundary_matrix ); + if( twist_pairs != std_pairs ) { std::cerr << "Error: twist and standard differ!" << std::endl; error = true; @@ -144,8 +150,12 @@ int main( int argc, char** argv ) std::cerr << "Error: chunk and row differ!" << std::endl; error = true; } - if( row_pairs != twist_pairs ) { - std::cerr << "Error: row and twist differ!" << std::endl; + if( row_pairs != ss_pairs ) { + std::cerr << "Error: row and block spectral sequence differ!" << std::endl; + error = true; + } + if( ss_pairs != twist_pairs ) { + std::cerr << "Error: block spectral sequence and twist differ!" << std::endl; error = true; } -- cgit v1.2.3