matrix manipulation of a large diagonally symmetric matrix to reduce computational time

Hi there, Given a large square matrix (maybe up to 20k X 20k). I would like to reduce computational time for calculating special functions (elliptic integrals). When I image up this matrix it is diagonally symmetric as shown below. This diagonal symmetric behaviour is as expected I've shown 3 cases below.

What is the most efficient approach to reduce computational time for a matrix with this symmetry? Perhaps remove half the matrix then do computation then replicate using repmat? If this can be done it will greatly reduce the computational time. Any suggestions will be very helpful

Thanks

case 1

case 2

case 3

[James Tursa]

If you decide to create a new variable with only the lower (or upper) triangle to work with, that means data copies in both directions. While this can be done at the m-file level, you may end up with large intermediate matrices in the process. This will slow things down and potentially create memory problems given your stated 20K x 20K size. These large intermediate variables can be avoided with a fairly simple mex routine. Do you have a C compiler available?

[Udit Gupta]

Are you doing an element by element operation or a matrix operation?

[Joseph]

Hi James, Yes I have a C compiler. I dont have much experience with using c/c++ should I look into that?

[Joseph]

Hi Udit, Below is what I'm doing to operate on the upper triangular half of the matrix where NT is the length (1 < NT <=20k+). I can perhaps set the preallocating matrix to sparse. I still have allocate before looping,..

    tic
k2ij=zeros(NT,NT);
MNTxNT=zeros(NT,NT);
toc
tic
for i=1:NT-1
    for j=i+1:NT
        k2ij(i,j)=4*ri(i)*ri(j)/((ri(i)+ri(j)).^2+(zi(i)-zi(j)).^2);
        MNTxNT(i,j) = sqrt((ri(i)*ri(j))/(k2ij(i,j))); %  need this later
    end
end
MNTxNT(logical(eye(size(MNTxNT)))) = 0;
k2ij(k2ij(:)==0) = []; % remove zero elements
MNTxNT(MNTxNT(:)==0) = []; % remove zero elements
toc
  end

tic
[K,ENTxNT] = ellipke(k2ij); % ellipke(k2ij,eps*1e14) <- helps abit speed wise
toc

[James Tursa]

Since you are creating the matrices from scratch with the zeros function, I have some suggestions:

1) Use the uninit function from the FEX to pre-allocate ... it is much faster than the zeros function for large matrices. You can find uninit here:

http://www.mathworks.com/matlabcentral/fileexchange/31362-uninit-create-an-uninitialized-variable--like-zeros-but-faster-

2) Do the operations in-place in a mex routine, which will prevent large intermediate copies. (I can do this coding for you)

I can work on the mex routine to do this over the next few days, but I have a question about the last three lines.

MNTxNT(logical(eye(size(MNTxNT)))) = 0;

If the diagonal starts as 0's because of the zeros pre-allocation, and you don't set the diagonal to anything in your loops (j index starts at i+1), why do you go through the trouble of setting it to 0's explicitly? It should already be 0's. Also, even if you insist on doing this, a much more efficient way of setting the diagonal to 0's would have been:

MNTxNT(1:NT+1:end) = 0;

And then these lines:

k2ij(k2ij(:)==0) = []; % remove zero elements
MNTxNT(MNTxNT(:)==0) = []; % remove zero elements

You are physically removing the 0 elements, which will turn the matrices into vectors (i.e., extracting the triangle). Is this what you want for downstream programming beyond the ellipke function? Or are you going to turn them to full again later on? If you really want them this way, then better to just build them that way from the beginning in the mex routine. Let me know how you really want this so I can start coding.

Basically, what explicitly do you want as variables after all of this is finished? Back to full square symmetic matrices? Knowing this will determine how the mex routines are coded.

[Joseph]

Thanks James, you are right I don't need the line that removes the zeros..

Actually my final result is thankfully only scalar so there is no need to reconstruct the matrix. I simply take the result and multiply the final scalar by 2 to account for the symmetry. I'll now try allocating using the method you suggested above. It maybe just easier for to pass you function .m file.

[Joseph]

I should add that my matrix size can range between 1x1 - 20k+20k+. Since I'm running the function through optimization.