source: branches/stable/TEMPLATES/matrix.h

// ========================================================= //
//                                                           //
//   File      : matrix.h                                    //
//   Purpose   : templates for matrices                      //
//                                                           //
//   Coded by Ralf Westram (coder@reallysoft.de) in Jul 20   //
//   http://www.arb-home.de/                                 //
//                                                           //
// ========================================================= //

#ifndef MATRIX_H
#define MATRIX_H

#ifndef ARB_ASSERT_H
#include <arb_assert.h>
#endif

template<typename T>
class symmetric_matrix : virtual Noncopyable {
    // Symmetrical Matrix (upper triangular matrix)

    size_t     Size;
    T **m;       // m[i][j]  i <= j !!!!

public:

    explicit symmetric_matrix(size_t Size_)
        : Size(Size_)
    {
        size_t elements = Size*(Size+1)/2;
        size_t headsize = Size * sizeof(*m);
        size_t datasize = elements * sizeof(*(m[0]));

        m    = (T**)ARB_calloc<char>(headsize+datasize);
        m[0] = (T*)(((char*)m)+headsize);

        for (size_t i=1; i<Size; i++) {
            m[i] = m[i-1]+i;
        }
    }

    ~symmetric_matrix() {
        free(m);
    }

    void set(size_t i, size_t j, T val) { if (i>j) m[i][j] = val; else m[j][i] = val; };

    T fast_get(size_t i, size_t j) const { arb_assert(i>=j); return m[i][j]; };
    T get(size_t i, size_t j) const { if (i>j) return m[i][j]; else return m[j][i]; };

    T get_max_value() const { // O(n*2)
        // Warning: ignores matrix diagonal

        T max = m[0][0]; // take any
        for (size_t i=0; i<Size; i++) {
            for (size_t j=0; j<i; j++) {
                if (m[i][j] > max) max = m[i][j];
            }
        }
        return max;
    }

    size_t size() const { return Size; }
};


#else
#error matrix.h included twice
#endif // MATRIX_H