b2matrix_structured.H

//------------------------------------------------------------------------
// b2matrix_structured.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 __B2MATRIX_STRUCTURED_H__
#define __B2MATRIX_STRUCTURED_H__
 
#include "b2linear_algebra_def.H"
 
namespace b2000::b2linalg {
 
template <typename MTYPE>
struct Mstructured_constref {
    using base = Mstructured_constref;
    using const_base = Mstructured_constref;
    using copy = typename MTYPE::copy;
};
 
template <typename T, typename MTYPE>
class Matrix<T, Mstructured_constref<MTYPE>> {
public:
    Matrix() : m(Matrix<T, MTYPE>::null), s(0), s2(0), index(0), index2(0) {}
 
    Matrix(size_t size_, const size_t* index_, const Matrix<T, MTYPE>& m_)
        : m(m_), s(size_), s2(size_), index(index_), index2(0) {}
 
    Matrix(
          size_t size_1, size_t size_2, const size_t* index1_, const size_t* index2_,
          const Matrix<T, MTYPE>& m_)
        : m(m_), s(size_1), s2(size_2), index(index1_), index2(index2_) {}
 
    Matrix(const Matrix& m_) : m(m_.m), s(m_.s), s2(m_.s2), index(m_.index), index2(m_.index2) {}
 
    std::pair<size_t, size_t> size() const { return std::pair<size_t, size_t>(s, s2); }
 
    size_t size1() const { return s; }
 
    size_t size2() const { return s2; }
 
    size_t esize() const { return s * (s + 1) / 2; }
 
    T operator()(size_t i, size_t j) const {
        const size_t* ii = std::find(index, index + m.size1(), i);
        if (ii == index + m.size1()) { return 0; }
        const size_t* jj;
        if (index2) {
            jj = std::find(index2, index2 + m.size1(), j);
            if (jj == index2 + m.size1()) { return 0; }
        } else {
            jj = std::find(index, index + m.size1(), j);
            if (jj == index + m.size1()) { return 0; }
        }
        return m(ii - index, jj - index);
    }
 
private:
    const Matrix<T, MTYPE> m;
    size_t s;
    size_t s2;
    const size_t* index;
    const size_t* index2;
    MVFRIEND;
};
 
template <typename T, typename MTYPE>
Matrix<T, Mstructured_constref<typename MTYPE::const_base>> StructuredMatrix(
      size_t size, const Index& index, const Matrix<T, MTYPE>& m) {
    if (index.size() != m.size1()) {
        Exception() << "The index size differ from the matrix size." << THROW;
    }
    return Matrix<T, Mstructured_constref<typename MTYPE::const_base>>(size, &index[0], m);
}
 
template <typename T, typename MTYPE>
Matrix<T, Mstructured_constref<typename MTYPE::const_base>> StructuredMatrix(
      size_t size1, size_t size2, const Index& index1, const Index& index2,
      const Matrix<T, MTYPE>& m) {
    if (index1.size() != m.size1() && index2.size() != m.size2()) {
        Exception() << "The index size differ from the matrix size." << THROW;
    }
    return Matrix<T, Mstructured_constref<typename MTYPE::const_base>>(
          size1, size2, &index1[0], &index2[0], m);
}
 
}  // namespace b2000::b2linalg
 
#endif