Skip to content

Fast matrix apply #73

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
dlg0 wants to merge 2 commits into
base: reference
Choose a base branch
from
Open

Fast matrix apply #73

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
44da516
Adding SG to FG mapping
stefan-schnake Feb 1, 2022
0687339
Changed backward_euler to reflect faster evaluation
stefan-schnake Feb 1, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion asgard_run_pde.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -536,7 +536,7 @@
else
figs.mass = figure('Name','mass(t)','Units','normalized','Position',[0.7,0.5,0.2,0.3]);
end
plot(outputs.mass_t/outputs.mass_t(1));
plot(mass_t/mass_t(1));
xtitle='timestep';
ytitle='mass/mass(t=0)';
ylim([0,2]);
Expand Down
76 changes: 76 additions & 0 deletions sg_to_fg_mapping_2d.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,76 @@
function [perm,iperm,pvec] = sg_to_fg_mapping_2d(pde,opts,A_data)
% Maps the sparse_grid index to standard tensor product full_grid index.
% This to allow multiplication by the full-grid kronecker sum stiffness
% matrix.

% perm -> SG-type index to FG-standard index
% pvec -> logical vector showing which elements are lit. Used for SG->FG
% conversion.
% iperm -> inverse mapping

% This reindexing allows us to use the standard sum of kroncker products
% stiffness matrix in lieu of the construction used in apply_A. Using the
% conversion to full-grid and back is also much faster.

% Works for adaptiviity

% Only works in 2D for now, but can be easily changed for higher diminsions

% Function assumes A_data.element_global_row_index is 1:sparse_dof.
% --Does not check for this

% USAGE: If the matrix A and vector f are in the SG-index, then
% we can compute A*f=f1 using A_F by the sequence
% f_F = zeros(size(A_F,1),1);
% f_F(pvec) = f(perm(pvec));
% Af_F = A_F*f_F;
% f1 = Af_F(iperm);
% where A_F is the standard full-grid sum of kronecker products.
% Specifically, A_F is generated by the code:
% A_F = 0;
% for i=1:numel(pde.terms)
% A_F = A_F + kron(pde.terms{i}.terms_1D{1}.mat,...
% pde.terms{i}.terms_1D{2}.mat);
% end
%
% The SG-index stiffness matrix and be computed from the full grid matrix
% by using:
% A = A_F(iperm,iperm);

%Checking for only 2 dimensions
assert(numel(pde.dimensions) == 2);

deg = opts.deg;

num_sparse_ele = numel(A_data.element_local_index_D{1}); %Number of sparse elements
iperm = zeros(num_sparse_ele*deg^2,1);

%Get FG_dof for each direction
max_x_dof = size(pde.terms{1}.terms_1D{1}.mat,1);
max_y_dof = size(pde.terms{1}.terms_1D{2}.mat,1);
perm = nan(max_x_dof*max_y_dof,1);

idx = 1; %Keep track of index
for i=1:num_sparse_ele
%First compute the x dim FG dof
dof_x = (A_data.element_local_index_D{1}(i)-1)*deg + (1:deg);
%Now y dim FG dof
dof_y = (A_data.element_local_index_D{2}(i)-1)*deg + (1:deg);

%Mapping is (i,j) -> (i-1)*max_y_dof+j;

for deg_x=1:deg
for deg_y=1:deg
iperm(idx) = (dof_x(deg_x)-1)*max_y_dof+dof_y(deg_y);
idx = idx + 1;
end
end
end


perm(iperm) = 1:numel(iperm);
pvec = ~isnan(perm);


end

172 changes: 19 additions & 153 deletions time_advance.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -251,162 +251,28 @@
s0 = source_vector(pde,opts,hash_table,t+dt);
bc0 = boundary_condition_vector(pde,opts,hash_table,t+dt);

if opts.time_independent_A || opts.time_independent_build_A
persistent A;
%Get SG <-> FG conversion
[~,iperm,~] = sg_to_fg_mapping_2d(pde,opts,A_data);

%Construct full-grid A and ALHS
A_F = 0;
for i=1:numel(pde.terms)
A_F = A_F + kron(pde.terms{i}.terms_1D{1}.mat,...
pde.terms{i}.terms_1D{2}.mat);
end

if opts.time_independent_A % relies on inv() so is no good for poorly conditioned problems
persistent AA_inv;
persistent ALHS_inv;

if isempty(AA_inv)
if applyLHS
[~,A,ALHS] = apply_A(pde,opts,A_data,f0,deg);
I = speye(numel(diag(A)));
ALHS_inv = inv(ALHS);
AA = I - dt*(ALHS_inv * A);
b = f0 + dt*(ALHS_inv * (s0 + bc0));
else
[~,A] = apply_A(pde,opts,A_data,f0,deg);
I = speye(numel(diag(A)));
AA = I - dt*A;
b = f0 + dt*(s0 + bc0);
end

if numel(AA(:,1)) <= 4096
condAA = condest(AA);
if ~opts.quiet; disp([' condest(AA) : ', num2str(condAA,'%.1e')]); end

if condAA > 1e6
disp(['WARNING: Using time_independent_A=true for poorly conditioned system not recommended']);
disp(['WARNING: cond(A) = ', num2str(condAA,'%.1e')]);
end
end

AA_inv = inv(AA);
f1 = AA_inv * b;
else
if applyLHS
b = f0 + dt*(ALHS_inv * (s0 + bc0));
else
b = f0 + dt*(s0 + bc0);
end
f1 = AA_inv * b;
end

else % use the backslash operator instead
if applyLHS
if opts.time_independent_build_A
if isempty(A);[~,A,ALHS] = apply_A(pde,opts,A_data,f0,deg);end
else
[~,A,ALHS] = apply_A(pde,opts,A_data,f0,deg);
end
I = speye(numel(diag(A)));
AA = I - dt*(ALHS \ A);
b = f0 + dt*(ALHS \ (s0 + bc0));
else
if opts.time_independent_build_A
if isempty(A);[~,A] = apply_A(pde,opts,A_data,f0,deg);end
else
[~,A] = apply_A(pde,opts,A_data,f0,deg);
end
I = speye(numel(diag(A)));
AA = I - dt*A;
b = f0 + dt*(s0 + bc0);
end

% Direct solve
rescale = false;
if rescale
% [AA_rescaled,diag_scaling] = rescale2(AA);
% AA_thresholded = sparsify(AA_rescaled,1e-5);
[P,R,C] = equilibrate(AA);
AA_rescaled = R*P*AA*C;
b_rescaled = R*P*b;
f1 = AA_rescaled \ b_rescaled;
f1 = C*f1;
else
f1 = AA \ b;
A = A_F(iperm,iperm);
ALHS = speye(size(A));
if applyLHS
ALHS_F = 0;
for i=1:numel(pde.termsLHS)
ALHS_F = ALHS_F + kron(pde.termsLHS{i}.terms_1D{1}.mat,...
pde.terms{i}.terms_1D{2}.mat);
end

% % Pre-kron solve
% num_terms = numel(pde.terms);
% num_dims = numel(pde.dimensions);
% for d=1:num_dims
% dim_mat_list{d} = zeros(size(pde.terms{1}.terms_1D{1}.mat));
% for t=1:num_terms
% dim_mat_list{d} = dim_mat_list{d} + pde.terms{t}.terms_1D{d}.mat;
% end
% end
% for d=1:num_dims
% A = dim_mat_list{d};
% I = speye(numel(diag(A)));
% AA_d = I - dt*A;
% dim_mat_inv_list{d} = inv(AA_d);
% end
% use_kronmultd = true;
% if use_kronmultd
% f1a = kron_multd(num_dims,dim_mat_inv_list,b);
% else
% f1a = kron_multd_full(num_dims,dim_mat_inv_list,b);
% end

% % Iterative solve
% restart = [];
% tol=1e-6;
% maxit=1000;
% tic;
% [f10,flag,relres,iter,resvec] = gmres(AA,b,restart,tol,maxit);
% t_gmres = toc;
% figure(67)
% semilogy(resvec);
%
% % Direct solve - LU approach to reducing repeated work
% [L,U,P] = lu(AA);
% n=size(AA,1);
% ip = P*reshape(1:n,n,1); % generate permutation vector
% err = norm(AA(ip,:)-L*U,1); % err should be small
% tol = 1e-9;
% isok = (err <= tol * norm(AA,1) ); % just a check
% disp(['isok for LU: ', num2str(isok)]);
% % to solve a linear system A * x = b, we have P * A * x = P*b
% % then from LU factorization, we have (L * U) * x = (P*b), so x = U \ (L \ ( P*b))
% f1_LU = U \ (L \ (P*b));
%
% % Direct solve - QR reduced rank
% n=size(AA,1);
% [Q,R,P] = qr(AA); % A*P = Q*R
% % where R is upper triangular, Q is orthogonal, Q’*Q is identity, P is column permutation
% err = norm( AA*P - Q*R,1);
% tol=1e-9;
% isok = (err <= tol * norm(AA,1)); % just a check
% disp(['isok for QR: ', num2str(isok)]);
% % to solve A * x = b, we have A * P * (P’*x) = b, (Q*R) * (P’*x) = b
% % y = R\(Q’*b), P*y = x
% tol=1e-9;
% is_illcond = abs(R(n,n)) <= tol * abs(R(1,1));
% if(is_illcond)
% disp('is_illcond == true');
% bhat = Q'*b;
% y = zeros(n,1);
% yhat = R(1:(n-1), 1:(n-1)) \ bhat(1:(n-1)); % solve with smaller system
% y(1:(n-1)) = yhat(1:(n-1));
% y(n) = 0; % force last component to be zero
% f1_QR = P * y;
% else
% disp('is_illcond == false');
% f1_QR = P * (R \ (Q'*b));
% end
%
% disp(['f1-f10: ',num2str(norm(f1-f10)/norm(f1))]);
% disp(['f1-f1_LU: ',num2str(norm(f1-f1_LU)/norm(f1))]);
% disp(['f1-f1_QR: ',num2str(norm(f1-f1_QR)/norm(f1))]);
% disp(['direct runtime: ', num2str(t_direct)]);
% disp(['gmres runtime: ', num2str(t_gmres)]);
%
% f1 = f1_QR;

ALHS = ALHS_F(iperm,iperm);
end

f1 = (ALHS-dt*A)\(ALHS*f0 + dt*(s0+bc0));

end


Expand Down