http://matrixprogramming.com/TAUCS/
Thank you for your help.
I have noted (it seems that I have missed it before) that TAUCS has the
libraries for ATLAS as well as for the reference BLAS. Please use only
the libraries from ATLAS. This way it is much faster - see comments in
the text.
Happy New Year,
Evgenii
Actually it should be exactly the same but as I have said taucs_linsolve
has a bug. If you are ready invest time in debugging it, I could give
you a hint. My guess is that this line in taucs_linsolve
int opt_factor = 1;
is wrong. It forces taucs_linsolve to go through the factorization both
times. Try to change it to
int opt_factor = 0;
recompile TAUCS and see what happens. Yet, this is only a guess, it well
might be other bugs as well.
taucs_linsolve actually uses the low-level function similar to the code
that I have presented. In my practice I am working with TAUCS with
low-level functions only and I do not use taucs_linsolve.
> in the code of test_taucs2.cpp i have some questions
>
> Question 1:
> in this example
> A=[1 0.5 0 0
> 0 1 0.5 0
> 0 0 1 0.5
> 0 0 0 1 ] (right?)
Actually not. The matrix to factorize is symmetric, so the original
matrix is
A=[1 0.5 0 0
0.5 1 0.5 0
0 0.5 1 0.5
0 0 0.5 1 ]
Otherwise the result (0 2 0 4) would be wrong. The form that you have
written just saves memory, as one keeps only half of a symmetric matrix.
Actually according to CCS (Compressed Column Storage) employed in TAUCS
one keeps the low triangle, that is
A=[1
0.5 1
0 0.5 1
0 0 0.5 1 ]
See for example
http://www.cs.utexas.edu/users/inderjit/Resources/sparse_matrices
but here the matrix is unsymmetric.
> then
> Could you please give the value of each variable in every step ? I
> don't konw how to check the value of A or prem or imvperm .so didn't
> understant the effect of each function
> // 1) Reordering
> taucs_ccs_order(&A, &perm, &invperm, "metis");
A is the original matrix
perm and invperm are two permutation vector. The memory for them is
allocate in this function automatically. You can release this memory at
the end of your program by free(perm) and free(invperm).
"metis" is an reordering scheme
See slide 18 in
http://modelreduction.com/doc/teaching/eurosime/lecture2.pdf
about reordering and comparison of different reordering schemes.
> Aod = taucs_ccs_permute_symmetrically(&A, perm, invperm);
This function uses the permutation vectors from above to permute the
matrix.
> taucs_vec_permute(dim, TAUCS_DOUBLE, b, bod, perm);
As the matrix was permuted one must also permute the right hand side.
> // 2) Factoring
> F = taucs_ccs_factor_llt_mf(Aod);
>
> // 3) Back substitution and reodering the solution back
> taucs_supernodal_solve_llt(F, xod, bod);
This is a solution for the permuted system
> taucs_vec_ipermute(dim, TAUCS_DOUBLE, xod, x, perm);
so it is necessary to permute it back to the original ordering.
The reordering step is crucial for direct solvers. Theoretically you can
solve the original system without permutation but it will take forever
and presumably memory of your computer will not be enough.
You will find more information about these function in the TAUCS manual.
But if you have more questions, please.
> Question 2:
>
> when we creat matrix A we use
> A.flags = (TAUCS_DOUBLE | TAUCS_SYMMETRIC | TAUCS_LOWER);
> Dose it means A is double ,and summetric and lower triangular?
Yes, provided that TAUCS_LOWER means that we keep the lower triangular
part of the symmetric matrix. This means that the original matrix is not
lower triangular. The latter is just a storage scheme. This again comes
from the TAUCS manual.
> why we
> use "|" rather than "&"
Because we would like to preserve in A.flags the bits from all three
macros. Print TAUCS_DOUBLE, TAUCS_SYMMETRIC, TAUCS_LOWER, and A.flags.
Then change |" to "&" and again print A.flags. You will see the difference.
> These Questions may seem stupid ,but it is useful to me,for I am not
> very famliar with programming.
No problem. Feel free to ask any questions. Yet, at the same time it is
also good start reading the manual.
If you aren't building a software library yourself and all you care
about is solving the problem, maybe you should use Matlab (which is
not free) or Octave (which is free both as in open-source and as in no
cost) instead. Both of these programs support sparse matrices and
sparse matrix factorizations, and both are much easier to use for
those who aren't familiar with programming.
mfh
I would agree. If the goal is to learn numerics, environments like
Matlab are a good starting point. I should say that I like Mathematica
more, but this is a matter of tastes. From the free software there is
also SciPy.
On the other hand, if you plan to go to high dimensional problems or
high efficiency is a must, then C or C++ are good choices. As far as I
know Matlab does not support symmetric storage for sparse matrices.
Mathematica also does not support it (at least this was the case in
Mathematica 5). And without the symmetric storage the work with
symmetric sparse matrices will not be efficient. Yet, it is necessary to
invest time in learning.
It depends on your matrix and on the particular application. If you
are very concerned about getting all the performance you can, then it
would be worth your time to try a few sparse factorization routines
and see which one is faster for typical problems of yours. I can
think of three such routines: Taucs, UMFPACK (which happens to be in
Matlab, so you can try it out there), and SuperLU. There are many
more. You should also try different orderings of the sparse matrix.
If you have a copy of Matlab, you can experiment there with different
orderings to see which one makes the factorization go faster, as
Matlab implements several different orderings. If you don't care
about small constant factors in performance, just take the solver
that's easiest to use. In either case it's worth looking at a few
different factorization routines.
Also, if your matrix comes from an optimization problem and you are
forming it as B * B' then it is probably symmetric. You may be able
to get some performance benefit from using a factorization that can
exploit symmetry. I would recommend a symmetric indefinite
factorization unless you know for sure that B * B' is positive
definite.
Best of luck!
mfh
A couple of more comments.
TAUCS is for symmetric positive definite matrices. If this is not the
case for you, just forget about TAUCS. TAUCS has a solver for indefinite
symmetric matrices but it is not multifrontal - as a result it is very slow.
If you have symmetric indefinite or unsymmetric, I would recommend you
MUMPS. It is much faster than UMFPACK. The problem with UMFPACK is that
it does not use METIS.
I should say that I have not tried SuperLU.
Evgenii