I have a function that requires a 3-element array and a 2-element array as inputs. Since its possible that each one would be an array block, my function looks like the following:

Code: Select all

template <typename Derived>
void solve_quadratic(const ArrayBase<Derived> &c, ArrayBase<Derived> &r, int &nre) {

When I call the function with an array created using

ArrayXd c(3,1);
ArrayXd r(2,1);

It works as expected;

However, when I try to call it using a fixed-size array, it fails.

Array<double,3,1> c;
Array<double,2,1> r;

Shouldn't ArrayBase be allowing for either one of these? On a more general note, I would expect it to be better to tend towards fixed-size arrays if possible, is this not the case?

Here's an example code, uncommenting the second call to solve_quadratic results in a compile-time error. ( gcc version 4.2.1 (Apple Inc. build 5666) )

Code: Select all

#include <Eigen/Dense>

using namespace Eigen;

template <typename Derived>
void solve_quadratic(const ArrayBase<Derived> &c, ArrayBase<Derived> &r, int &nre) {
   /**
   * Solves for x in
   * c[0]*x**2 + c[1]*x + c[2] = 0
   *
   *     If the quadratic has two real solutions x1 and x2 they are returned as
   *     r[0] = x1
   *     r[1] = x2
   *     nre  = 2
   *
   *     If the quadratic has no real solutions x +- i*y they are returned as
   *     r[0] = x
   *     r[1] = y
   *     nre  = 0
   *
   * Arguments:
   *     c_in - The coefficients c as given in the equation above
   *     r    - On output, returns the roots of r
   *     nre  - On output, returns the number of real roots in r
   **/
   double q, A, B;

    q = pow(c(1),2) - 4.0*c(0)*c(2);
    A = -0.5*c(1)/c(0);
    B = 0.5*c(0);

    if (q >= 0) {
        nre = 2;
        r(0) = A + B * sqrt(q);
        r(1) = A - B * sqrt(q);
    } else {
        nre = 0;
        r(0) = A;
        r(1) = B * sqrt(abs(q));
    }
}




int main( int argc, char* argv[]) {
   using namespace std;
    ArrayXd c(3,1);
    ArrayXd r(2,1);
    Array3d c2;
    Array2d r2;
    int nre;
   
    c << 0.5, 0, -20.0;
    c2 << 0.5, 0.0, -20.0;
    solve_quadratic(c,r,nre);
    //solve_quadratic(c2,r2,nre);  // COMPILE TIME ERROR
    /**
    test.cpp: In function 'int main(int, char**)':
    test.cpp:58: error: no matching function for call to 'solve_quadratic(Eigen::Array3d&, Eigen::Array2d&, int&)'
    **/   
    cout << "Quadratic has " << nre  << " real roots" << endl;
    cout << r << endl;
   return 0;
}

Manhigh wrote:Shouldn't ArrayBase be allowing for either one of these? On a more general note, I would expect it to be better to tend towards fixed-size arrays if possible, is this not the case?

ArrayBase does allow for fixed-size arrays. The problem is that the size is encoded in the template parameter. Array3f inherits from ArrayBase<Array3f> and Array2f inherits from ArrayBase<Array2f> (this is called CRTP = Curiously Recurring Template Pattern). In your function, both the array arguments are of type ArrayBase<Derived>, so they have to be of the same type. So you should change the prototype to

Code: Select all

template <typename Derived1, typename Derived2>
void solve_quadratic(const ArrayBase<Derived1> &c, ArrayBase<Derived2> &r, int &nre)


And hopefully it will work then. I don't understand your second question; if this does not answer it then please ask again.

Thanks, that was indeed the problem.

The second part of the question was more about the efficiency of dynamic vs fixed arrays. I would suspect that fixed size arrays are fast, but if you're not reallocating the dynamically sized array, I'm not sure. It's something that I can play around with as my project matures.

fixed sizes also allow for explicit unrolling and so much faster performance for very small matrices (e.g., 2, 3, 4).