- 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. - Hi Uwe,

you wrote:

> I use the atlas-bindings for linear algebra and noticed that

Let's start from the end:

> 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...

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.