## atlas problem

Expand Messages
• 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.
• 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

/* 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?

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.