Boost logo

Ublas :

Subject: Re: [ublas] todo list
From: Matwey V. Kornilov (matwey.kornilov_at_[hidden])
Date: 2010-07-24 08:39:40


As far as I understand when we talk about implementation of such algorithms
in uBLAS we mean that API should be 'ublas-style'. For instance, when we
write inner_prod( v, u ) we get an object of this operation which looks like
a scalar for us. The question is what is the 'ublas-style' for SVD,
Cholesky, LU, QR, etc. algorithms? It is not clear for me.

Most of the algorithms are in-place. For instance, we have to alter our
matrix in order to get its SVD decomposition. The next example. Cholesky
decomposition is defined only for positive-defined matrices. When matrix is
not positive-defined there is a square root of negative number in the
algorithm. The good C++ style is to throw exception here and to leave the
matrix in initial state (i.e not altered). But known Cholesky implementation
assumes undefined matrix at output in such cases.

For instance the result of SVD decomposition is a triple of matrices: one
singular diagonal S matrix and two unitary matrices U and V. Imagine we want
to know only singular values to calculate conditioning number of matrix or
it's rank. What happens when we don't want to know U and V matrices? As
uBLAS user I expect that decomposition saves calculations and memory for
placing unitary matrices. What happens when we want to use SVD for matrix
inversion? In this case we have to know U and V and thus they should be
calculated and allocated.

The third question is the estimations of rounding errors.

I suppose there are a several dozens of questions left.

I tried to implement SVD and Cholesky for my purposes. The code is
available:
  http://curl.sai.msu.ru/~matwey/lsp
Subversion:
  svn://curl.sai.msu.ru/lsp
But IMO the API definitely is not the good one.

Marco Guazzone wrote:

> # SVD, Cholesky, LU, QR, Shur
> Do you think this should be written from scratch or maybe using
> existent solvers like the one provided by LAPACK (possibly using the
> boost::numeric::bindings lib)?