Trying to handle a large matrix of nx101 dimensions where n is unknown at start of loop, and I only need the previous vector to calculate the current vector

So, I am writing a code for calculating the temperature curve along a 1 dimensional rod as a function of time until steady state, and the numerical method I am using only needs the temperature curve from the previous time step (which is stored in the matrix), as well as the previous iteration within the current time step (which is thrown out after calculation), but I need to keep the values from all previous time steps even though they are no longer relevant to the calculation. The matrices that I am storing get extremely large, and therefore Matlab starts to slow down, I was wondering if there were any tricks I could use.

[John BG]

1.- USE DISK INSTEAD OF RAM:

If the 2 initial answers loading variables to RAM don't work for you, because you RAM limitation, then store/retrieve data on disk.

There are multiple options:

• save to/read from .mat file
• save to/read from .txt file
• spreadsheet will take longer read write

2.- THRESHOLD TOO SMALL:

Also, a difference threshold between temperatures of 1e-6?

isn't it a bit too strict?

you may have already reached a steady state, you don't know, the gap between temperatures may be 1e-5 or 1e-3, and you assume your simulation is crashing when perhaps a threshold of 0.01 Celsius would be ok.

If not asking for too much, if using disk instead of RAM your memory please let us know how you do it

John

[e-mail]

[Jacob Bresler]

It turns out the issue with the speed the code was running at was due to the expansion of the matrices within the T and N cell arrays which was solved by the GrowData2 function that was suggested by John D'Errico. However the code does get dangerously close to taking up all of the available RAM on my computerter and I was wondering if you could clarify how to save to/read from a .mat file, when I googled it the page wasn't very clear. As to changing the threshold, the small value is necessary since the change between time steps is actually relatively small even when it isn't at steady state (it takes on the order of 400-500 hours for this particular equation to reach steady state and the time step is 1 second)

[Jacob Bresler]

I modified the code so that the largest matrix (the one that holds all of the Temperature Vectors) gets moved from RAM to the Disk after it finishes for each for loop, this dramatically reduces the RAM necessary, I used a .mat file.

    %Jacob Bresler
    tic
    clear
    clc
    h_prime=.05;
    sigma_0_prime=2.7e-9;
    T_0=300;
    T_1=400;
    T_inf=200;
    rho=7870;
    c_p=447;
    k=80.2;
    L=10;
    h=.1; %x step
    K=1; %time step
    x=0:h:L;
    sx=floor(L/h);
    denom1=2*k*K+k*K*(h^2)*h_prime+rho*c_p*(h^2); %denominator for left half of bar, not time dependent
    beta=[.01,.1,1,10];
    N=cell(1,4);
    t_ss=cell(4,3);
    filename='Temp_1s';
    M=matfile(filename);
    M.Properties.Writable=true;
    s_T=zeros(1,4);
    T=cell(1,4);
    for b=1:4
    T_s=growdata2;
    T_i=T_0:(T_1-T_0)*h/L:T_1;
    T_s(T_i);
    N_s=growdata2;
    T_diff_t=1;
    T_diff_n=1;
    t=0;
    a=0;
    n=0;
    T_loop=T_i;
    T_t_old=T_i;
        while max(T_diff_t)/K>1e-6 %loop checking for steady state, corrects for time step size
            a=a+1;
            t=t+K;
            t_idx=floor(t/K)+1;
            while max(T_diff_n)>1e-6 %loop checking for convergence within time step
                n=n+1;
                T_n_old=T_loop(1,:); %vector for T(J+1)^n
                for i=2:sx/2+1
                    num1=k*K*(T_loop(i-1)+T_n_old(i+1)+(h^2)*h_prime*T_inf)+(h^2)*rho*c_p*T_t_old(i);
                    T_loop(1,i)=num1/denom1;
                end
                for i=sx/2+2:sx
                    num2=k*K*(T_loop(i-1)+T_n_old(i+1)+(h^2)*(h_prime*T_inf+sigma_0_prime*exp(-beta(b)*t)*(T_inf^4)))+(h^2)*rho*c_p*T_t_old(i);
                    denom2=k*K*(2+(h^2)*(h_prime+sigma_0_prime*exp(-beta(b)*t)*(T_n_old(i)^3)))+(h^2)*rho*c_p;
                    T_loop(1,i)=num2/denom2;
                end
                T_diff_n=abs(T_loop-T_n_old);
            end
            T_diff_n=1;
            T_diff_t=abs(T_loop-T_t_old);
            N_s(n);
            n=0;
            T_t_old=T_loop;
            T_s(T_loop);
            [a,b]
        end
        T_temp=T_s();
        M.T(1:t_idx,1+(sx+1)*(b-1):(sx+1)*b)=T_temp;
        s_T(b)=t_idx;
        clear T_temp
        N{1,b}=N_s();
        T_diff_t=1;
        f_h=floor(t/3600);
        r_h=round(t/3600);
        if abs((t/3600-r_h))<1e-6
            t_ss{b,1}=r_h;
        else
            t_ss{b,1}=f_h;
        end
        m=(t/3600-floor(t/3600))*60;
        f_m=floor(m);
        r_m=round(m);
        if abs((m-r_m))<1e-6
            t_ss{b,2}=r_m;
        else
            t_ss{b,2}=f_m;
        end
        s=(m-t_ss{b,2})*60;
        r_s=round(s);
        if abs(s-r_s)<1e-6
            t_ss{b,3}=r_s;
        else
            t_ss{b,3}=s;
        end
        t=0;
        a=0;
    end
    toc

Answer 1

[John D'Errico]

Easy, peasy (with the right tools.)

As an example, I'll simulate 101 parallel random walks in 1 dimension, stopping only when one of them goes outside of a preset bound, but I'll store all the walks in one array. At the very end, I'll unpack the function handle the data was stored in.

tic
% initialization step
funh = growdata2;
n = 101;
walklimit = 1000;
vec = zeros(1,n);
funh(vec);
theMoonIsBlue = true;
steps = 0;
while theMoonIsBlue
  steps = steps + 1;
  % Each step uses only the previous step.
  vec = vec + randn(1,n);
  funh(vec)
  if any(abs(vec) > walklimit)
    theMoonIsBlue = false;
  end
end
allwalkdata = funh();
toc

Elapsed time is 9.192757 seconds.

9 seconds is way better than the same computation would have taken had I grown an array dynamically, adding one new row to the array at every step.

steps
steps =
    189995

So almost 200k steps in the random walk.

Where did these walks terminate?

allwalkdata(end,:)
ans =
Columns 1 through 11
    -636.58      -541.39      -1000.9      -436.61       356.79      -157.98      -161.44       478.41      -65.325      -624.95      -795.86
Columns 12 through 22
    -569.76      -169.58       37.705      -373.35       -193.1      -316.91      -69.306       -430.6      -140.92       247.08      -143.95
Columns 23 through 33
    -339.97      -59.332       173.08       288.98      -78.927      -546.23        26.17       81.455       25.087       399.69       83.686
Columns 34 through 44
     520.92      -158.04        690.7       430.35       837.11       604.51       -459.4       369.44       47.953       17.497       265.34
Columns 45 through 55
    -73.708       171.93      -502.52      -436.02       213.63      -327.65        439.4      -595.36       -689.5      -609.44      -487.53
Columns 56 through 66
     56.602      -341.29       391.56      -345.06       104.16          324       407.17       396.58      -186.36       320.52       127.96
Columns 67 through 77
    -126.18       515.95       326.52      -160.58      -291.39      -671.77       375.42       161.05      -215.69      -709.72       -169.1
Columns 78 through 88
     442.65       206.96      -366.62       296.44       -68.06       85.838       180.54       93.865        343.9       36.668       196.27
Columns 89 through 99
     3.2942       253.73       503.08       482.23      -498.44       149.07      -615.69       144.99       55.706        501.4      -208.14
Columns 100 through 101
     -306.1      -342.34

You can find the growdata tools on the file exchange .

[Jacob Bresler]

that looks like it might work, will check it out and respond

[Jacob Bresler]

I haven't run it long enough to see for sure, but it looks like it is running into the same issue, it is just storing that huge matrix somewhere else but still in RAM, which is starting to get crowded, on previous runs when it hadn't even gotten close to convergence for the first of the four beta values, the matlab process was exceeding 1GB in size just from storing the matrix. I need a method for removing the matrix from RAM but not losing the data.

[John D'Errico]

Of course it stores the data. That is what you said you wanted to do!

Those tools avoid the problem of growing a matrix in memory, because they don't reallocate the array at every step. So these tools try to avoid the nasty quadratic time penalty of continuously reallocating an array in memory.

You said that you wanted to store the array. If it is so huge that you cannot store it in RAM (you did not emphasize that in your original question) then you will need to write out to a file on disk. Frequent disk writes won't be fast either though.

You need to accept that if you need to keep that much data around, then you need some place to put it, and you need to accept a time penalty.

[Jacob Bresler]

I ended up finding a solution (I save it to the HDD and then clear the matrix). Takes a bit longer since writing to the HDD is much slower than RAM but it saves a lot of RAM. Just in case you run into an issue sometime in the future that needs a solution. here is the code you would need lets say you want to store a matrix x in a matfile as x

M=matfile('filename')
%and then just 
M.x=x;
clear x

Answer 2

[dpb]

The primary one is to preallocate the array and populate rather than augmenting it. That is, start by

N=1000;                  % initial storage size
n=101;                   % solution vector size
output=zeros(N,n);       % preallocate
iter=0;                  % iteration counter
while ~steadystate       % the logic to determine solution end
  if mod(iter,N)==0
    output=[output;zeros(N,n)];    % allocate additional
  end
  iter=iter+1;           % increment counter
  output(iter,:)=solution(output(iter-1,:));  % new value
end

This should still be nearly invariant w/ iteration altho it will take a little longer when augmenting the array it won't be nearly as much overhead as the ploy of

output=[output;solution(previousoutput)];

that reallocates every pass thru the loop.

[Jacob Bresler]

preallocation doesn't help in this case because I have absolutely no Idea what the magnitude of n will be before I start the code. as to the other option, does that remove it from RAM? because my issue is that the code itself is taking up too much RAM and is starting to slow down due to that fact.

[dpb]

I don't understand...why is the size indeterminate? Are you using a variable mesh algorithm, perhaps?

It's still so that preallocation will help but you do need to be able to bound the array size, yes.

If you're running out of RAM and are running into swapping or the like, that's another issue that has to do with you'd do better to archive computed results to disk every so often and only keep in memory what is required for the immediate purpose.

Pseudo code showing the implementation now might help to explain your problem and would, if well done, let folks here see where there might be savings/optimizations.

[Jacob Bresler]

the T array is where the temperature vectors are stored as a 101xn matrix, and the size of the matrices within that array grows with every loop until the error criteria (max(T_diff_t)<1e6) is met

%Jacob Bresler
  tic
  clear
  clc
  h_prime=.05;
  sigma_0_prime=2.7e-9;
  T_0=300;
  T_1=400;
  T_inf=200;
  rho=7870;
  c_p=447;
  k=80.2;
  L=10;
  h=.1; %x step
  K=1; %time step
  x=0:h:L;
  sx=floor(L/h);
  denom1=2*k*K+k*K*(h^2)*h_prime+rho*c_p*(h^2); %denominator for left half of bar, not time dependent
  beta=[.01,.1,1,10];
  T=cell(1,4); %this is the storage point for the Temperature vectors
  for j=1:4
      T{1,j}(1,:)=T_0+(x/L).*(T_1-T_0);
  end
  t_ss=cell(4,3);
  N=cell(1,4);
  T_diff_t=1;
  T_diff_n=1;
  t=0;
  a=0;
  n=0;
  for b=1:4
      while max(T_diff_t)/K>1e-6 %loop checking for steady state, corrects for time step size
          a=a+1;
          t=t+K;
          t_idx=floor(t/K)+1;
          T_loop=T{1,b}(t_idx-1,:); %vector for T(j+1)^n+1
          T_t_old=T{1,b}(t_idx-1,:); %vector for T(j)
          while max(T_diff_n)>1e-6 %loop checking for convergence within time step
              n=n+1;
              T_n_old=T_loop(1,:); %vector for T(J+1)^n
              for i=2:sx/2+1
                  num1=k*K*(T_loop(i-1)+T_n_old(i+1)+(h^2)*h_prime*T_inf)+(h^2)*rho*c_p*T_t_old(i);
                  T_loop(1,i)=num1/denom1;
              end
              for i=sx/2+2:sx
                  num2=k*K*(T_loop(i-1)+T_n_old(i+1)+(h^2)*(h_prime*T_inf+sigma_0_prime*exp(-beta(b)*t)*(T_inf^4)))+(h^2)*rho*c_p*T_t_old(i);
                  denom2=k*K*(2+(h^2)*(h_prime+sigma_0_prime*exp(-beta(b)*t)*(T_n_old(i)^3)))+(h^2)*rho*c_p;
                  T_loop(1,i)=num2/denom2;
              end
              T_diff_n=abs(T_loop-T_n_old);
          end
          T{1,b}(t_idx,:)=T_loop;
          T_diff_n=1;
          T_diff_t=abs(T_loop-T_t_old);
          N{1,b}(a)=n;
          n=0;
          [a,b]
      end
      T_diff_t=1;
      f_h=floor(t/3600);
      r_h=round(t/3600);
      if abs((t/3600-r_h))<1e-6
          t_ss{b,1}=r_h;
      else
          t_ss{b,1}=f_h;
      end
      m=(t/3600-floor(t/3600))*60;
      f_m=floor(m);
      r_m=round(m);
      if abs((m-r_m))<1e-6
          t_ss{b,2}=r_m;
      else
          t_ss{b,2}=f_m;
      end
      t_ss{b,3}=(m-t_ss{b,2})*60;%seconds (decimal)
      t=0;
      a=0;
  end
  toc

Answer 3

[Jacob Bresler]

I added a tic toc to the loop checking for convergence in time, and it is showing that the code does indeed slow down at a considerable rate as loops even with the addition of the growdata2 function. It may be running faster than previously, but it still slows down at a considerable rate. Is there a way to fully remove the superfluous vectors from RAM but not lose the values?

Edit: Ok, sorry, I found a leftover line of code that was actually causing it to still create the matrix on the fly within my code, I fixed it and it is looking better now