b2linear_algebra_utils.H

//------------------------------------------------------------------------
// b2linear_algebra_utils.H --
//
//     A framework for generic linear algebraic (matrix) computations.
//
// written by Mathias Doreille
//
// Copyright (c) 2004-2012 SMR Engineering & Development SA
//                         2502 Bienne, Switzerland
//
// All Rights Reserved.  Proprietary source code.  The contents of
// this file may not be disclosed to third parties, copied or
// duplicated in any form, in whole or in part, without the prior
// written permission of SMR.
//------------------------------------------------------------------------
 
#ifndef __B2LINEAR_ALGEBRA_UTILS_H__
#define __B2LINEAR_ALGEBRA_UTILS_H__
 
#include <cassert>
 
#include "b2linear_algebra_def.H"
#include "b2ppconfig.h"
#include "utils/b2logging.H"
 
namespace b2000 {
 
using std::size_t;
 
namespace b2linalg {
 
// ====================================
// temporary buffer for sparse operation
// ====================================
 
#ifdef HAVE_LIB_TBB
#define thread_local_storage __thread
#else
#define thread_local_storage
#endif
 
template <typename T>
class TemporaryBuffer {
public:
    static T* get(size_t size) {
        if (buffer_size < size) {
            while (buffer_size < size) { buffer_size = buffer_size ? buffer_size * 2 : 1; }
            delete[] buffer;
            buffer = new T[buffer_size];
            std::fill_n(buffer, buffer_size, initialisation_value);
        }
        return buffer;
    }
 
private:
    static thread_local_storage size_t buffer_size;
    static thread_local_storage T* buffer;
    static T initialisation_value;
};
 
template <typename T>
thread_local_storage size_t TemporaryBuffer<T>::buffer_size = 0;
 
template <typename T>
thread_local_storage T* TemporaryBuffer<T>::buffer = 0;
 
template <typename T>
T TemporaryBuffer<T>::initialisation_value = T(0);
 
template <>
size_t TemporaryBuffer<size_t>::initialisation_value;
 
// ====================================
// Slice for sub Vector/Matrix
// ====================================
 
class Interval {
public:
    Interval(size_t start_, size_t end_) : start(start_), end(end_) {}
 
    bool is_null() const { return this == &null; }
 
    bool is_all() const { return this == &all; }
 
    size_t size() const { return end - start; }
 
    size_t operator[](size_t i) const { return start + i; }
 
    size_t get_start() const { return start; }
 
    size_t get_end() const { return end; }
 
    std::pair<bool, size_t> find(size_t i) const {
        if (i >= start && i < end) { return std::pair<bool, size_t>(true, i - start); }
        return std::pair<bool, size_t>(false, 0);
    }
 
    static Interval null;
    static Interval all;
 
protected:
    MVFRIEND;
    friend class Index;
    size_t start;
    size_t end;
};
 
class Slice {
public:
    Slice(size_t start_, size_t end_ = 0, size_t step_ = 1)
        : start(start_), end(end_ == 0 ? start_ + 1 : end_), step(step_) {}
 
    bool is_null() const { return this == &null; }
 
    bool is_all() const { return this == &all; }
 
    size_t size() const { return (end - start) / step; }
 
    size_t increment() const { return step; }
 
    size_t operator[](size_t i) const { return start + i * step; }
 
    std::pair<bool, size_t> find(size_t i) const {
        if (i >= start && i < end) {
            size_t r = i - start;
            if ((r % step) == 0) { return std::pair<bool, size_t>(true, r / start); }
        }
        return std::pair<bool, size_t>(false, 0);
    }
 
    static Slice null;
    static Slice all;
 
protected:
    MVFRIEND;
    size_t start;
    size_t end;
    size_t step;
};
 
class Index : public std::vector<size_t> {
public:
    Index(size_t s = 0) : std::vector<size_t>(s) {}
 
    Index(const Index::const_iterator& ibegin, const Index::const_iterator& iend)
        : std::vector<size_t>(ibegin, iend) {}
 
    Index(size_t start, size_t end) {
        reserve(end - start);
        while (start != end) { push_back(start++); }
    }
 
    bool is_null() const { return this == &null; }
 
    bool is_all() const { return this == &all; }
 
    void set_zero() { std::fill(begin(), end(), 0); }
 
    std::pair<bool, size_t> find(size_t i) const {
        std::vector<size_t>::const_iterator ii = std::lower_bound(begin(), end(), i);
        if (*ii == i) { return std::pair<bool, size_t>(true, ii - begin()); }
        return std::pair<bool, size_t>(false, 0);
    }
 
    bool sorted() const {
        for (size_t i = 1; i != size(); ++i) {
            if (operator[](i - 1) > operator[](i)) { return false; }
        }
        return true;
    }
 
    Index operator()(const Interval& i) { return Index(begin() + i.start, begin() + i.end); }
 
    Index make_dual() const {
        Index index(size());
        for (size_t i = 0; i != size(); ++i) {
            const size_t ii = operator[](i);
            assert(ii < size());
            index[ii] = i;
        }
        return index;
    }
 
    static Index null;
    static Index all;
 
    friend logging::Logger& operator<<(logging::Logger& l, const Index& v) {
        l << "index vector ";
        l.write(v.size(), &v[0], 1);
        return l;
    }
};
 
std::ostream& operator<<(std::ostream& out, const Index& index);
 
}  // namespace b2linalg
}  // namespace b2000
 
#endif