Reverse Cuthill McKee Ordering and Solution of systems of Linear Equations

Question

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.

Answer 1

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).