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 600800. > > >> 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_{t1}  grad(F(X_{t1})) > > > 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.numanalysis?!? 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_{t1} + grad(X, G(X_{t1},Y_{t1}) (A+aI)Y_{t} = aY_{t1} + grad(Y, G(X_{t1},Y_{t1})
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_{t1} (Y_{t1}).
Now, what i don't remember (the paper does not speak about this) is the update of the stepsize <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 semiimplicit 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

