luca
Posts:
11
Registered:
5/19/10
|
|
Re: Sparse linear system
Posted:
May 19, 2010 11:24 AM
|
|
On 19 Mag, 16:54, Nicolas Neuss <lastn...@kit.edu> wrote:
> luca <luca.frammol...@gmail.com> writes:
> > On 19 Mag, 12:07, luca <luca.frammol...@gmail.com> wrote:
> >> Hi,
>
> >> suppose i have a matrix M X N of reals. This matrix is sparse. Every
> >> row has only 3 non zero values (they are always -1, 2, -1). The M and
> >> N are on the order of 600-800.
>
> >> Is there a fast library (way) to solve a sparse linear system A X = B,
> >> where the matrix A is of the type defined below??
>
> >> I need a library writte in C/C++ that is possibly portable.
>
> >> Thank you,
> >> Luca
>
> > A is a NxN real sparse symmetric matrix. Sorry for the confusion. I
> > will give you guys more details. The problem i am facing is the
> > following. Given a NxN real sparse symmetric matrix A, solve the
> > following system of linear equation:
> > (A + aI)X_{t} = aX_{t-1} - grad(F(X_{t-1}))
>
> > where <a> is a real scalar (the step size), X_{t} is a vector of N
> > reals, F is differentiable function, I is the identity matrix.
> > One start with an arbitrary X_0 and <a> = some constant and every
> > iteration involves the solution of this system of linear equations and
> > the update of the scalar <a> using some schedule.
>
> > Since A is a big sparse matrix, A+aI is a big sparse matrix too. So i
> > could use a LU factoriztion of A+aI and than solve
> > the system in the usual way.
>
> > A does not change during the iterations, but the coefficient matrix is
> > A+aI, so i need to compute the LU of A+aI.
> > I would really like to avoid the LU factorization at every iteration,
> > but i don't see any alternative. Maybe there is a way
> > to compute LU of A and than, given aI, compute (in a fast way) the LU
> > of A+aI using the LU of A...
>
> > So, what i need is a free portable C/C++ code for computing a fast LU
> > factorization of a sparse matrix.
>
> > Thank you,
> > Luca
>
> If you have a symmetric tridiagonal matrix with constant coefficients,
> the decomposition will be of the form L D L^t with L having only one
> subdiagonal and ones on the diagonal and D is diagonal. Furthermore,
> the entries of D and L can be calculated recursively which is very
> cheap. Maybe you do not even need to store them.
>
> Nicolas- Nascondi testo citato
>
> - Mostra testo citato -
Ok, first thank you all for your patience. What happened to
sci.math.num-analysis?!? It is full of spam!
Second, i would like to give you more details. Since my exam of num
analysis was some years ago ;) maybe
you can help me understand better.
Suppose we have the function F(X,Y). We need the solution to the PDE:
grad(X, F(X,Y)) = 0
grad(Y, F(X,y)) = 0
where grad(X, F) (and grad(Y,F)) is the gradient of F with respect to
X (to Y). Ok, the paper i am reading suggest to solve this system in
the following way:
(A+aI)X_{t} = aX_{t-1} + grad(X, G(X_{t-1},Y_{t-1})
(A+aI)Y_{t} = aY_{t-1} + grad(Y, G(X_{t-1},Y_{t-1})
where G is a differentiable function.
A is the symmetric, positive but not definite, real, big, sparse
matrix, <a> is a real parameter (the step size) and so on (read my
previous messages).
It is and iterative method. Start with random X_0, Y_0 and stop when
X_t (Y_t) is not much different from X_{t-1} (Y_{t-1}).
Now, what i don't remember (the paper does not speak about this) is
the update of the step-size <a>. I don't know if <a> has to be updated
(in some way) at every iteration or if it is constant.
This makes a big difference. If <a> is constant, i would compute once
for all the LU of (A+aI) and than solve, at every iteration, a system
of linear
equations. Otherwise, if <a> is not constant, i need to do LU at every
iteration (and find a way to update it).
The paper call this method semi-implicit and cite this other paper:
http://www.cs.ucla.edu/~dt/papers/ijcv88/ijcv88.pdf
If you look at the last page (the appendix), you will see the same
equations above.
So, is <a> constant?
Is there a C/C++ free LU factorization library for spare big
matrices??
Thank you again.
Best regards,
Luca
|
|
