1
$\begingroup$

I just learned about the RCM. I am trying to solve a problem that is a result of fluid dynamics and chemistry so I have a very large sparse matrix. I also learned reducing the bandwidth would accelerate the solver.

How do I know that the Matrix reordered matrix A in the system Ax=b will yield the same x vector as in the original A ?

I played a bit with Octave and Matlab to test it, but I don't see that is the case. So, I must understood something wrong, or did something wrong. Can someone here shed some light on the subject?

Well, it seems the answer is YES, but I don't know how to do it.

Here is an example done with octave:

Ny=900; Nx=900; ex = ones(Nx,1); Axx = spdiags([ex -2*ex ex], -1:1, Nx, Nx); Ix = speye(Nx); ey = ones(Ny,1); Ayy = spdiags([ey -2*ey ey], -1:1, Ny, Ny); Iy = speye(Ny); A = kron(Axx,Iy)+kron(Ix,Ayy); C=symrcm (A); %this is the actual RCM % just to visualize things  spy (A) spy (A(C,C)) b=ones(900,1); % store the permutated matrix Amut= (A(C,C)); % solve and compare x=A\b; xmut=Amut\b; x==xmut % results in avector full of zeros, e.g. False. 

So how do I permutate my xmut to be x again ? it seems the individual entries are correct, but not in the right order.

It seems that the values in the vectors are equal if I do that:

octave:136>r1=x(C) % permutate x ... octave:138> xmut(1:10) ans =    -2.0039   -3.5078   -3.5078   -4.7093   -6.3179   -4.7093   -5.7014   -8.6281   -8.6281   -5.7014  octave:139> r1(1:10) ans =    -2.0039   -3.5078   -3.5078   -4.7093   -6.3179   -4.7093   -5.7014   -8.6281   -8.6281   -5.7014  octave:140> isequal(r1(1:10),xmut(1:10)) ans = 0 

SO I think the answer is YES, I still need to find how to compare the vectors...

So the final answer: octave:144> r1-xmut

yields only numbers like -1.3323e-15, -4.4409e-16 which practically can be treated as zero because of the floating points limitations. Which means my vectors are the same. Which is what I wanted to find out.

1 Answers 1

2

It looks like you managed to convert the non-permuted to permuted with the line

r1 = x(C) % permutate x ... 

so I believe all you need now is a way to change permuted to non-permuted. One way to do this is

x(C) = A(C,C) \ b(C) 

which just immediately writes the correct values to the solution array x. Let me explain how this works. Say that the first entry is in C is 10 (so the 10th entry will be 1). Then the 10th value of x is set to be the first value in A(C,C) \ b(C), which is what we want. This is a result of how slicing is implemented in MATLAB and Octave (I only tested this in MATLAB, but Octave should do the same thing).