Re: how do i solve sparse linear equation like A*x=b

Yan Yaqin,
See the example below. Best regards

On Thursday 02 September 2004 11:52, Yan Yaqin wrote:
> hi,
> I'm really appreciate for your help, but there's another thing.
> One parameter for lu_substitute is vector_expression, i'm not familiar with
> this type, Would you like to provide me a simple real number demo?
>
> Thanks and best regards!

> ls -l
total 8
lrwxrwxrwx  1 leopardi users   25 2004-09-03 00:24 boost 
-> ../boost/boost-1_31/boost
-rw-r--r--  1 leopardi users 1457 2004-09-03 00:41 sample_solve.cpp
-rw-r--r--  1 leopardi users 1458 2004-09-03 00:40 sample_solve.cpp~

> g++ -I. sample_solve.cpp -o sample_solve

> ./sample_solve
A: [2,2]((0,1),(2,0))
B: [2,2]((1,0),(0,-1))
X: [2,2]((0,-0.5),(1,0))

> cat sample_solve.cpp
#include <fstream>

#include <boost/numeric/ublas/config.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/vector_sparse.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/operation_sparse.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/lu.hpp>

namespace ublas = boost::numeric::ublas;

typedef ublas::row_major                           orientation_t;
typedef ublas::compressed_matrix< double, orientation_t >
                                                   matrix_t;
typedef matrix_t::size_type                        matrix_index_t;
typedef ublas::permutation_matrix<matrix_index_t>  permutation_t;

// Geometric quotient
matrix_index_t
sample_solve (matrix_t& X, const matrix_t& A, const matrix_t& B)
{
  // We solve A*X == B

  matrix_t LU = A;
  permutation_t pvector(A.size1());
  matrix_index_t singular = ublas::lu_factorize(LU, pvector);

  if (! singular)
  {
    X = B;
    ublas::lu_substitute(LU, pvector, X);
  }
  return singular;
}

main()
{
  matrix_t A(2,2);
  matrix_t B(2,2);
  matrix_t X(2,2);

  A(0,1) =     1;
  A(1,0) = 2;
  std::cout << "A: " << A << std::endl;

  B(0,0) = 1;
  B(1,1) =    -1;
  std::cout << "B: " << B << std::endl;

  sample_solve(X,A,B);
  std::cout << "X: " << X << std::endl;
}
>

_______________________________________________
This list is archived at http://www.osl.iu.edu/MailArchives/mtl-devel/