Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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.

share|improve this question
add comment

2 Answers

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

share|improve this answer
add comment
up vote 0 down vote accepted

Well,I'm just answering so the question is not left unanswered.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.