[ublas] lu factorize/substitute question

From: Prasad Ramanan (prasad.ramanan_at_[hidden])
Date: 2007-04-30 12:05:01

Hi,

I have a question regarding the lu_factorize/lu_substitute functions
in ublas. I actually wrote to this list with a similar question in Nov
2005 (and was helped a great deal by Gunter Winkler), so I apologize
in advance for a repetitive query, although my new question is (I
think) regarding a different problem.

I wrote the following code to solve Ax = rhs:

#include <iostream>
#include <fstream>
#include <string>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/lu.hpp>

### includes suggested by Gunter (thanks for your help!)
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/triangular.hpp>

using namespace boost::posix_time;
using namespace std;

void matrixTest() {
        namespace ublas = boost::numeric::ublas;

        const n = 10;
        ublas::matrix<double> A(n,n);
        ublas::vector<double> x(n);
        ublas::vector<double> rhs(n);
        for(int i = 0; i < n; i++) {
                A(i,i) = 1;
                rhs(i) = i;
        }

        ublas::permutation_matrix<double> P(n);
        ublas::lu_factorize(A,P);
        x = rhs;
        ublas::lu_substitute(A,P,rhs);
}

This code compiles, and I was asked by Gunter Winkler to add the
includes marked by my comment since at the time I was using boost
1-32. Otherwise, the code did not compile (compiler errors are
documented in one of my previous messages to this list; I am using
Visual Studio 7.1.3088).

I have since upgraded to boost 1-33.1 (downloaded from
SourceForge.net), but I still require the includes suggested by Gunter
in order to get the code to compile, although he had suggested that I
would no longer need them after upgrading to 1-33. Once the code is
compiled, I get a runtime error caused by the lu_substitute method
while running:

Assertion failed in file
c:\boost\include\boost-1_33_1\boost/numeric/ublas/lu.hpp at line 271:
detail::expression_type_check (prod
(triangular_adaptor<const_matrix_type, upper> (m), e), cv2)

I am not sure what is happening here, since the example is actually
quite simple in terms of linear algebra. Of course my code is
contained inside a much larger project that relies on other boost
libraries, but I am not sure that is relevant here.

Again, I apologize for a slightly repetitive query, and thanks in
advance for any help that you can provide.

Best,
Prasad