Efficient Band Matrix Computations

Listing 1: Template class bandStorage

#ifndef BANDSTOR_H
#define BANDSTOR_H

#include <iostream.h>
#include <math.h>
#include <vector.h>

inline int MAX(int i, int j) { return (i>j ? i : j);}
inline int MIN(int i, int j) { return (i<j ? i : j);}

template <class T> class bandStorage
{
private:
    int n;
    int upperMostBand;
    int lowerMostBand;

// data for the matrix

    vector< vector <T> > a;

public:
    bandStorage() : n(0),
        upperMostBand(0),
        lowerMostBand(0),
        a(0) {}
    bandStorage(const int);
    bandStorage& operator=(const bandStorage&);
    inline int order() const {return n;};
    inline int& upperBandWidth(void) {return upperMostBand;}
    inline int& lowerBandWidth(void) {return lowerMostBand;}
    vector<T>& diag(const int);
    const vector<T>& diag(const int) const;
    T& operator()(const int, const int);    
    vector<T> operator*(const vector<T>&) const;
};

template <class T>
    bandStorage<T>::bandStorage(const int order)
{
    n = order;
    upperMostBand = 0;
    lowerMostBand = 0;
    a = vector< vector <T> >(2*n-1);
}

template <class T>
    vector<T>& bandStorage<T>::diag(const int location)
{
    return a[location+n-1];
}

template <class T>
    const vector<T>& bandStorage<T>::diag(const int location) const
{
    return a[location+n-1];
}

template <class T>
T& bandStorage<T>::operator()(const int i, const int j)
{
    static T zero = T(0.0);
     if (j <= (i+upperMostBand) && i <= (j-lowerMostBand))
        if (i <= j) // upper triangle
            return a[j-i+n-1].operator[](i);
        else
            return a[j-i+n-1].operator[](j);
    else
        return zero;
}

template <class T>
vector<T> bandStorage<T>::operator*(const vector<T>& x) const
{
    vector<T> y(n, T(0.0));
    vector<T> currentBand(n);
    int i, j;
    for (i=lowerMostBand; i<=upperMostBand; i++) {
        currentBand = diag(i);
        if (i <= 0)
            for (j=-i; j<n; j++)
                y[j] += currentBand[j+i]*x[j+i];
        else
            for (j=0; j<n-i; j++)
                y[j] += currentBand[j]*x[j+i];
    }
    return y;
}

#endif
/* End of File */