Thread Subject: eigs for larger FEM matrices, out of memory

I have two FEM (finite element) matrices A and M of the size around 1 milion. Both are symetric and positive definite coming from 3D discetization of an elliptic equation.
The command

options.disp=0;
few_largest_eigenvalues=eigs(M,A,1,'LM',options);

unfortunately fails giving me back

Out of memory. Type HELP MEMORY for your options.

Error in ==> eigs>checkInputs/CHOLfactorB at 1047
                        [RB,idxB] = chol(B(ppB,ppB));

Error in ==> eigs>checkInputs at 877
            [RB,RBT,permB,BisHpd] = CHOLfactorB;

Error in ==> eigs at 96
[A,Amatrix,isrealprob,issymA,n,B,classAB,k,eigs_sigma,whch, ...

In fact, I would only need the larges eigenvalue or ist reasonable approximation. Now, I am just wondering whether, I should edit eigs fuction and modify it or whether there is an easier solution to it. My desktop has 8Gb memory, so with some tricks, it might hopefully just fit in. For 2D FEM discretization I can even compute eigenvalues for matrices of size 4 milion on the same machine.

Any suggestions,
Jan Valdman
University of Iceland

"Jan Valdman" <Jan.Valdman@gmail.com> wrote in message <hegkek$cc9$1@fred.mathworks.com>...

> In fact, I would only need the larges eigenvalue or ist reasonable approximation.

If an approximation will do, and if the matrices are reasonably concentrated around their diagonals, you might try to approximate the largest eignenvalue as

max(diag(M)./diag(A))

>My desktop has 8Gb memory, so with some tricks, it might hopefully just fit in.

It would take some very fancy tricks to give MATLAB access to the whole 8Gb. There's no way that I know of to go higher than 3Gb (using the \3Gb switch)

"Matt " <xys@whatever.com> wrote in message <heguqg$5ba$1@fred.mathworks.com>...
> "Jan Valdman" <Jan.Valdman@gmail.com> wrote in message <hegkek$cc9$1@fred.mathworks.com>...
>
> > In fact, I would only need the larges eigenvalue or ist reasonable approximation.
>
> If an approximation will do, and if the matrices are reasonably concentrated around their diagonals, you might try to approximate the largest eignenvalue as
>
> max(diag(M)./diag(A))
>
> >My desktop has 8Gb memory, so with some tricks, it might hopefully just fit in.
>
> It would take some very fancy tricks to give MATLAB access to the whole 8Gb. There's no way that I know of to go higher than 3Gb (using the \3Gb switch)

Do you think going to 64 bit make eigs(), svds(), etc. successful in this case?

Thanks Matt,

           M is a finite element mass matrix and A is a stifness matrix, they are indeed allocated around diagonal but there are offdiagonal terms as well (so that there is so called strictly diagonally dominancy). Therefore, the bound max(diag(M)./diag(A)) is very inacurate, unfortunatelly.

          In fact, I have 64 bit version of Matlab using all 8 Gb.

          Jan
          

"Jan Valdman" <Jan.Valdman@gmail.com> wrote in message <hegvs8$b7c$1@fred.mathworks.com>...
> Thanks Matt,
>
> M is a finite element mass matrix and A is a stifness matrix, they are indeed allocated around diagonal but there are offdiagonal terms as well (so that there is so called strictly diagonally dominancy). Therefore, the bound max(diag(M)./diag(A)) is very inacurate, unfortunatelly.
========

Then perhaps use the Gershgorin circle theorem. You can analyze the Gershgorin discs to get an upper bound on the eigenvalues of M and a lower bound on the eigenvalues of A, and then take their ratio. I can't say how tight a bound this will be (it depends on the strength of the strict diagonal dominance), but it will at least take into account the off-diagonal terms...

 

See this thread: just apply Krylov/Arnoldi iteration.
http://www.mathworks.com/matlabcentral/newsreader/view_thread/247873#639456

Replace the the fixed point assignment "x = A\x" by "x = A\(M*x)" (I believe).

Bruno

Dear Bruno,

        your advice helped me already. I just implemented
        
function lambda=krylov_iteration(M,A)
tol=1e-6;
xold=0;
x=zeros(length(A),1);
x(1)=1;
while norm(x-xold)>tol
    xold=x;
    x = A\(M*x);
    lambda=1/norm(x);
    x=lambda*x;
end
lambda=1/lambda;

which computes the maximal generalized eigenvalue of M x = lambda A x. For very small matrices, it converges and is faster than eigs function of Matlab. However, for larger matrices it is slower that eigs and for my matrix of size around 1 000 000 x 1 000 000 it again does not fit in the memory again:

??? Error using ==> mldivide
Out of memory. Type HELP MEMORY for your options.

Error in ==> krylov_iteration at 8
    x = A\(M*x);

I guess, I have to use an iterative method (CG, PCG) for solving x instead, I can do it.

One more question. Where do you plug in lambda, if you know a good approximation of it?

.
>
> One more question. Where do you plug in lambda, if you know a good approximation of it?

% Do something silimar to eigs(M,A,1,lbdguess)

tol=1e-6;
xold=0;
x=rand(length(A),1);
while norm(x-xold)>tol
    xold=x;
    for i=1:2
        x = (M - lbdguess*A) \ (A * x);
        if i==1
            s = sign(dot(x,xold));
        end
        a = norm(x);
        x = x/a;
    end
end

lambda = 1./(s*a)+lbdguess

% Bruno