how to make my own gaussian elimination spdiags

Question

Kiana Bahrami 2019-12-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/494243-how-to-make-my-own-gaussian-elimination-spdiags

Assignment6.m

Hi,

I'm trying to solve a sparse matrix using my own spdiag function instead of matlab's built in function. I have the following code to work with and the attached question file.

clear
clc
close all
n = 20;
dval = [-1 -1 4 -1 -1];
dloc = [-3 -1 0 1 3];
b = zeros(n,1);
b(1,1) = 1;
b(n,1) = -1;
nd = length(dloc);
Asp = zeros(n,nd);
% test
e = ones(n,1);
A = (spdiags(dval.*e,dloc,n,n))
for ind = 1:nd
    if dloc(ind)==0
        break
    end
end
for i = 1:ind
    Asp(1-dloc(i):n,i) = dval(i)
end
for i = ind+1:nd
    Asp(1:n-dloc(i),i) = dval(i)
end
[Asp,x] = sparse_GE(Asp,b,dloc)
[A1,x] = gauss_eli(A,b)
function [Asp,x] = sparse_GE(Asp,b,dloc)
[n,nd] = size(Asp);
x = zeros(n,1);
dmax = max(dloc);
dmin = min(dloc);
nd_new  = dmax-dmin+1;
dloc_new = dmin:dmax;
Asp_new = zeros(n,nd_new);
for i = 1:nd
    for j  = 1:nd_new
        if (dloc(i)==dloc_new(j))
           Asp_new(:,j) = Asp(:,i);
        end
    end
end
Asp = Asp_new;
dloc = dloc_new;
nd = nd_new;
% you can use the find function also
for ind = 1:nd_new
    if dloc_new(ind)==0
        break
    end
end
% Forward Elimination
for col = 1:n-1
    for i = ind-1:-1:1
        row = col-dloc(i);
        if row>n
            break
        end
        fac = Asp(row,i)/Asp(col,ind);
        Asp(row,i) = fac;
        for j = ind+1:nd
            m = i+j-ind;
            if (m>nd)
                break
            end
            Asp(row,m) = Asp(row,m) - fac*Asp(col,j);
        end
        b(row) = b(row) - fac*b(col);
    end
end
end
%% Gaussian Elimination Non Vectorized
function [A,x] = gauss_eli(A,b)
n = length(A); % we know it is a square system
% Forward Elimination
for col = 1:n-1
%    [A,b] = pivot(A,b,col,n);
    for row = col+1:n
        fac = A(row,col)/A(col,col);
        % non vectorized loop, to avoid cpu threading
        for k = col+1:n
            A(row,k) = A(row,k) - fac*A(col,k);
        end
        A(row,col) = fac;
        b(row) = b(row) - fac*b(col);
    end
end
% Back Substitution
for i = n:-1:1
    temp = 0;
    for j = i+1:n
        temp = temp+A(i,j)*x(j);
    end
    x(i,1) = (b(i)-temp)/A(i,i);
end
end