1

I'm revamping an old application that use Numerical Recipes' dmatrix quite extensively. Since one of the reasons I'm working on the application is because its code is about to be opened, I want to replace all of the Numerical Recipes code with code that can be freely distributed.

dmatrix is a function that returns a matrix of doubles. The called supplies the lower and upper bound for each index, like so:

double **mat = dmatrix(1,3,1,3);

mat now has 3 rows, from 1 to 3, and 3 columns, from 1 to 3, so that mat[1][1] is the first element and mat[3][3] is the last.

I looked at various C++ matrix implementations, none of them allowed me to specify the lower bound of each dimension. Is there something I can use, or do I have to write yet another matrix class for this?

zmbq
  • 38,013
  • 14
  • 101
  • 171
  • 1
    Never seen already existing implementation (but I hope there should be one). If you want to write a new implementation by yourself, you can inherit it from an existing implementation to reduce the cost of testing. – Mohit Jain Sep 23 '14 at 10:01
  • Writing a similar class in C++ is relatively simple, I'll have to test it for performance, though, because it'll have to be as fast as a `double**` eventually. – zmbq Sep 23 '14 at 10:03
  • 1
    Yuck :-/ Can you change the rest of the code to use lower-bound zero instead? Under the covers the old code must be over-allocating to make the indexes work on plain double pointers. – Rup Sep 23 '14 at 10:06
  • I don't want to do that. There are a lot of calculations in the old code, all of them rely on 1-based matrices. Change the code to start at 0 would require *a lot* of testing. – zmbq Sep 23 '14 at 10:10

1 Answers1

1

I believe that you can easily make a wrapper of some other matrix implementation to add the lower bound feature. Example (untested):

class Matrix {
    OtherMatrix m;
    int lowerX, lowerY;
public:

    Matrix(int lx, int hx, int ly, int hy) :
        m(hx-lx, hy-ly),
        lowerX(lx), lowerY(ly) { }

    MatrixCol operator[] (int x) {
        return {this, x};
    } 
};

class MatrixCol {
    friend class Matrix;
    Matrix* mm;
    int x;
public:
    double& operator[] (int y) {
        return mm->m[x - mm->lowerX, y - mm->lowerY];
    } 
};

This may require a little more robust implementation depending on your use case. But this is the basic idea, expand from it.

Guilherme Bernal
  • 8,183
  • 25
  • 43
  • Yes, that's more or less what I'm doing, only without another Matrix class, I'm using a std::vector of std::vector (which is basically what NR does, only with C arrays). – zmbq Sep 23 '14 at 11:58
  • 1
    @zmbq: Technically, NR doesn't use C arrays as those are zero-indexed. They're relying on rather dubious pointer abuse. – MSalters Sep 23 '14 at 12:04
  • Exactly. And it is this pointer abuse I'm trying to mimic in C++. I'll replace the other NR calculations, but other code is based on this same pointer abuse. – zmbq Sep 23 '14 at 12:31