Vectorize matrix multiplication
I have a school project about running a SOR algorithm on Octave, but mine is very inefficient. So I have this snippet of code:
for ii=1:n r = 1/A(ii,ii); for jj=1:n if (ii!=jj) A(ii,jj) = A(ii,jj)*r; end; end; b(ii,1) = b(ii,1)*r; x(ii,1) = b(ii,1); end;
How can I vectorize this? My first attempt was this:
for ii=1:n r = 1/A(ii,ii); A(find(eye(length(A))!=1)) = A(find(eye(length(A))!=1))*r; b(ii,1) = b(ii,1)*r; x(ii,1) = b(ii,1); end;
But I'm not sure it helped much. Is there a better and/or more efficient way to do it?
You can totally avoid loops I believe. You have to see that you are taking as the reciprocal of diagonal elements of A, then why do you want to use a loop. Do it directly. That's the first step. Remove Inf and now you want to multiply r to corresponding non-diagonal elements of corresponding rows, right?
Therefore, use repmat to construct such a matrix which will have replicated elements along columns, since you are multiplying same r to (1,2), (1,3), ... , (1,n). But R has non-zero diagonal elements. Therefore make them zero. Now you will get your A except that the diagonal elements will be zero. Therefore, you just have to add them back from original A. That can be done by A=A.*R+A.*eye(size(A,1)).
Vectorization comes from experience and most importantly analyzing your code. Think at each step whether you want to use the loop, if not replace that step with the equivalent command, other code will follow (for example, I constructed a matrix R, whereas you were constructing individual elements r. So I only thought about converting r -> R and then the rest of the code will just fall into place).
The code is as follows:
R=1./(A.*eye(size(A,1))); %assuming matrix A is square and it does not contain 0 on the main diagonal R=R(~isinf(R)); R=R(:); R1=R; R=repmat(R,[1,size(A,2)]); R=R.*(true(size(A,1))-eye(size(A,1))); A=A.*R+A.*eye(size(A,1)); %code verified till here since A comes out to be the same b = b.*R1; x=b;