How can I speed up the multiplication of matrices?

Hello all,

I am having problems trying to improve the speed of my code. I wrote a code to solve the pressure and water saturation of a reservoir. I have to iterate up to 500 days (with increments of 1 day), solving the pressure implicitly and the water saturation explicitly. See the code please:

      clc; clear;
    tic
    warning('off', 'MATLAB:singularMatrix');
    %%P R O P E R T I E S   &   C O M P U T A T I O N    O F    M A T R I C E S
    % Calls the function inputvalueSN for the input files
    [reservoir, fluid, numerical, well, pressure, rperm, Sw] = inputProject2_2;

    % Computes the matrices T, B, Q, G
    Sw_hyst = zeros(size(reservoir.permx));
    [T, Tw, d12, d11, D, J, Q, Qwprod, G, gamma_w] = TDJQG2(reservoir, fluid, numerical, well, pressure, rperm, Sw, Sw_hyst);

    %%N U M E R I C A L   S O L U T I O N S
    % BEFORE 500 DAYS
    P_old = sparse(pressure.Po);
    P_new = sparse(zeros(size(pressure.Po)));
    Po_total = zeros(length(pressure.Po), numerical.final_t/numerical.dt + 1);
    Po_total(:,1) = pressure.Po;

      Sw_total = zeros(length(Sw), numerical.final_t / numerical.dt + 1);
      Sw_total(:,1) = Sw;
      t = 1;

      for i = 0:numerical.dt:15 %numerical.switch - 1   
          P_new = (T + J + D) \ (D * P_old + Q + G); % It computes the oil pressure
          Sw_new = (Sw_total(:,t) + d12 \(-d11 * (P_new - P_old)  + Qwprod...
              - Tw * (P_new - gamma_w .* reservoir.depth - pressure.Pc)));

          indexzero = (reservoir.depth == 0);
          Sw_new(indexzero) = 0; P_new(indexzero) = 0;

          [pressure, Sw_hyst2] = Pr_hyst(Sw_total(:,t), Sw_new, rperm, P_new, pressure, Sw_hyst); % Update capillary pressure and Pw
          Sw_hyst = Sw_hyst2;

          % Updating the matrices and vector dependent of pressure and saturation
          [T, Tw, d12, d11, D, J, Q, Qwprod, G, gamma_w] = TDJQG2(reservoir, fluid, numerical, well, pressure, rperm, Sw_new, Sw_hyst);

        P_old = P_new;
        Po_total(:, t+1) = P_new;
        Sw_total(:, t+1) = Sw_new;
        t = t + 1;

    end

However, the matrices T, J, D, Q and G are pressure dependent, i.e., I have to update them every time step within the loop. The function that computes them is: TDJQG2. The script is the following:

function [T, Tw, d12, d11, D, J, Q, Qwprod, G, gamma_w] = TDJQG2(reservoir, fluid, numerical, well, pressure, rperm, Swater, Sw_hyst)
% This function computes the transmisibility matrix T, accumulation matrix
% B, source vector Q, G and J

%%FVF of oil and water
fluid.Bo = fluid.Bo_sat * (1 - fluid.co .* (pressure.Po - fluid.Pb));
fluid.Bw = fluid.Bw_ref * (1 - fluid.cw .* (pressure.Pw - pressure.ref_Bw));

% fluid.Bw = fluid.Bw_ref * (1 - fluid.cw .* (pressure.Pw - pressure.ref - pressure.Pe));

%%Update densities of oleic phase
oleic_dens = (fluid.density_o + fluid.Rs * (0.0458171 /5.615)) ./ fluid.Bo;

%%TDJQG
Tw = zeros(6000);
To = zeros(6000);

D = zeros(6000);
d11 = zeros(6000);
d12 = zeros(6000);
d21 = zeros(6000);
d22 = zeros(6000);

Jw = zeros(6000);
Jo = zeros(6000);

Qw = zeros(6000, 1);
Qo = zeros(6000, 1);
Qwprod = zeros(6000, 1);

for L = 1:6000
    if rem(L, numerical.Nx) ~= 1 % If the block is NOT on left edge
        [tw, to] = inter4(L, L-1, reservoir, fluid, numerical, pressure, rperm, Swater, oleic_dens); % Computation of transmissibilities

          % Water Transmissibilities on the left
          Tw(L, L-1) = -tw;   Tw(L, L) = Tw(L, L) - Tw(L, L-1);
          % Oil Transmissibilities on the left
          To(L, L-1) = -to;   To(L, L) = To(L, L) - To(L, L-1);        
      end

      %%%%%%%%MORE CODE THAT WORKS %%%%%%%%%%%%%%

end

d12 = sparse(d12); d11 = sparse(d11); d22 = sparse(d22); D = sparse(D);
C = isnan(D);
D(C) = 0;

Tw = sparse(Tw * 6.33e-3); To = sparse(To * 6.33e-3);
T = ((-d22 * d12^-1) * Tw + To) * 6.33e-3;

Qw = sparse(Qw); Qo = sparse(Qo);
Q = (-d22 * d12^-1) * Qw + Qo;
Qwprod = sparse(Qwprod);

Jw = sparse(Jw); Jo = sparse(Jo);
J = ((-d22 * d12^-1) * Jw + Jo) * 6.33e-3;

gamma_o = ((fluid.density_o + (fluid.Rs * fluid.density_g *  (1/5.615))) ./ fluid.Bo)./144;
gamma_o = sparse(gamma_o);
gamma_w = (fluid.density_w ./ fluid.Bw) ./144;
gamma_w = sparse(gamma_w);

prePc = sparse(pressure.Pc);
depths = sparse(reservoir.depth);

G = ((-d22 * d12^-1) * (Tw * prePc) + diag(-d22 / d12) .*...
    (Tw * (gamma_w .* reservoir.depth))+ (To * (gamma_o .* depths)));
G = sparse(G);
end

My main issue is that the following lines consume 57% of all the computation time. (I clic "Run and Time" for 15 days).

G = ((-d22 * d12^-1) * (Tw * prePc) + diag(-d22 / d12) .*... (Tw * (gamma_w .* reservoir.depth))+ (To * (gamma_o .* depths)));

d12 = sparse(d12); d11 = sparse(d11); d22 = sparse(d22); D = sparse(D);

T = ((-d22 * d12^-1) * Tw + To) * 6.33e-3;

Q = (-d22 * d12^-1) * Qw + Qo;

J = ((-d22 * d12^-1) * Jw + Jo) * 6.33e-3;

Is there any way to improve the computation of those matrices mentioned above? d22, d12, d11 are diagonal matrices, Q only has 10 entries, and J only 10 entries. G is a column vector. I already tried the backslash \ operator to compute those value but I get the message that the matrices are baldy scalded. The solutions are correct, but the is poor in performance.

Thank you in advance.

PD: these are the speed performance results:

[Walter Roberson]

Please do not close questions that have an answer.

Answer 1

[James Tursa]

Q1: Why do you pre-allocate a bunch of 6000x6000 full matrices for d11, d12, and d22 only to downstream use sparse( ) on them and tell us they are diagonal? Why not just create the diagonals to begin with and work with that?

Q2: Are you sure things stay diagonal throughout your calculations, or do small off-diagonal elements creep in? This might make a big difference in how the linear algebra calculations are done for the sparse matrices.

Q3: Why do you do the (-d22 * d12^-1) calculation repeatedly? Why not do it only once?

Q4: Are any of your large calculations involving a mix of sparse and full matrices? This could cause calculations to be done in full and slow things down (because the 0's are explicitly multiplied).

[Jose Salazar]

Hello James,

Q1: I am now working with the diagonals only. I am using spdiags.

Q2: There are some values where the main diagonal is just full of zeros. Like 20%.

Q3: I have not thought about! It is indeed helpful.

Q4: Should I multiply everything in sparse?

Thank you!