Loading ...
Sorry, an error occurred while loading the content.

atlas problem

Expand Messages
  • Uwe Schmitt
    Hi, I use the atlas-bindings for linear algebra and noticed that I have to use column_major ordering for correct results (lu_solve)... Is this is a bug ? I
    Message 1 of 2 , Jun 7, 2003
    • 0 Attachment
      Hi,

      I use the atlas-bindings for linear algebra and noticed that
      I have to use column_major ordering for correct results (lu_solve)...
      Is this is a bug ? I expected row_major should work as it
      is the storage organization used by LAPACK...

      Greetings, Uwe.
    • Kresimir Fresl
      Hi Uwe, ... Let s start from the end: I am afraid that you are wrong: row_major is not storage organization used by LAPACK. LAPACK is written in Fortran, and
      Message 2 of 2 , Jun 7, 2003
      • 0 Attachment
        Hi Uwe,

        you wrote:

        > I use the atlas-bindings for linear algebra and noticed that
        > I have to use column_major ordering for correct results (lu_solve)...
        > Is this is a bug ? I expected row_major should work as it
        > is the storage organization used by LAPACK...

        Let's start from the end:

        I am afraid that you are wrong: row_major is not storage
        organization used by LAPACK. LAPACK is written in
        Fortran, and Fortran matrices are column_major.

        But ATLAS is written in C, and it allows both row_major
        and column_major ordering. And atlas bindings are
        smart enough ;o) to recognize the ordering.

        Now, regarding wrong results: without your code or, at
        least, your linear system(s), I can only guess. Namely,
        both orderings work well for me (see the example at
        the bottom).

        The only reason I can think of is the problem which is
        explained in the leading comments for `gesv()' in clapack.hpp:

        // [comments from `clapack_dgesv.c':]
        /* clapack_xgesv computes the solution to a system of linear equations
        * A * X = B,
        * where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
        */
        // [but ATLAS FAQ says:]
        /* What's the deal with the RHS in the row-major factorization/solves?
        * Most users are confused by the row major factorization and related
        * solves. The right-hand side vectors are probably the biggest source
        * of confusion. The RHS array does not represent a matrix in the
        * mathematical sense, it is instead a pasting together of the various
        * RHS into one array for calling convenience. As such, RHS vectors are
        * always stored contiguously, regardless of the row/col major that is
        * chosen. This means that ldb/ldx is always independent of NRHS, and
        * dependant on N, regardless of the row/col major setting.
        */
        // That is, it seems that, if B is row-major, it should be NRHS-by-N,
        // and RHS vectors should be its rows, not columns.

        But if you didn't read this comment, and if you did the obvious
        thing (`obvious' -- ATLAS FAQ notwithstanding ;o) -- i.e. if your
        A is 3x3, and your B is 3x1 (and not 1x3) -- there should be a
        run-time error:

        Parameter 8 to routine clapack_dgesv was incorrect
        ldb must be >= MAX(N,1): ldb=1 N=3

        Only if your B is 3x3 or 3x4 etc., that is, if you have more
        than 3 rhs `vectors', you will get wrong results. Is this the
        case?

        If not, can you post your program (and your matrices)?

        Regards,

        fres


        Simple example:
        ===========================================
        #include <cstddef>
        #include <iostream>
        #include <boost/numeric/bindings/atlas/cblas.hpp>
        #include <boost/numeric/bindings/atlas/clapack.hpp>
        #include <boost/numeric/bindings/traits/ublas_matrix.hpp>
        #include <boost/numeric/ublas/io.hpp>

        namespace ublas = boost::numeric::ublas;
        namespace atlas = boost::numeric::bindings::atlas;

        using std::size_t;
        using std::cout;
        using std::endl;

        #ifndef F_ROW_MAJOR
        typedef ublas::matrix<double, ublas::column_major> m_t;
        #else
        typedef ublas::matrix<double, ublas::row_major> m_t;
        #endif

        int main() {

        size_t n = 3, nrhs = 1;

        m_t a (n, n);
        a(0,0) = 1.; a(0,1) = 1.; a(0,2) = 1.;
        a(1,0) = 2.; a(1,1) = 3.; a(1,2) = 1.;
        a(2,0) = 1.; a(2,1) = -1.; a(2,2) = -1.;

        #ifndef ROW_MAJOR
        m_t b (n, nrhs);
        b(0,0) = 4.; b(1,0) = 9.; b(2,0) = -2.;
        #else
        m_t b (nrhs, n);
        b(0,0) = 4.; b(0,1) = 9.; b(0,2) = -2.;
        #endif

        cout << "A: " << a << endl;
        cout << "B: " << b << endl;

        atlas::lu_solve (a, b);
        cout << "X: " << b << endl;
        }
        ===============================================

        Without -DROW_MAJOR output is:

        A: [3,3]((1,1,1),(2,3,1),(1,-1,-1))
        B: [3,1]((4),(9),(-2))
        X: [3,1]((1),(2),(1))

        and with -DROW_MAJOR:

        A: [3,3]((1,1,1),(2,3,1),(1,-1,-1))
        B: [1,3]((4,9,-2))
        X: [1,3]((1,2,1))

        Example (with solution) is taken from Dorn & McCracken: Numerical
        Methods with Fortran IV Case Studies (Wiley, 1972) -- an old book ;o),
        but the only one with solved examples which I can find now.
      Your message has been successfully submitted and would be delivered to recipients shortly.