Laplacian, Octave - rextester

Language: Layout:

function [lambda, V, A] = laplacian(varargin)

% Input/Output handling.
    if (nargin < 1 || nargin > 3)
      print_usage;
    endif

u = varargin{1};
    dim2 = size(u);
    if dim2(1) ~= 1
        error('BLOPEX:laplacian:WrongVectorOfGridPoints',...
            '%s','Number of grid points must be in a row vector.')
    end
    if dim2(2) > 3
        error('BLOPEX:laplacian:WrongNumberOfGridPoints',...
            '%s','Number of grid points must be a 1, 2, or 3')
    end
    dim=dim2(2); clear dim2;

uint = round(u);
    if max(uint~=u)
        warning('BLOPEX:laplacian:NonIntegerGridSize',...
            '%s','Grid sizes must be integers. Rounding...');
        u = uint; clear uint
    end
    if max(u<=0 )
        error('BLOPEX:laplacian:NonIntegerGridSize',...
            '%s','Grid sizes must be positive.');
    end

if nargin == 3
        m = varargin{3};
        B = varargin{2};
    elseif nargin == 2
        f = varargin{2};
        a = whos('regep','f');
        if sum(a.class(1:4)=='cell') == 4
            B = f;
            m = 0;
        elseif sum(a.class(1:4)=='doub') == 4
            if dim == 1
                B = {'DD'};
            elseif dim == 2
                B = {'DD' 'DD'};
            else
                B = {'DD' 'DD' 'DD'};
            end
            m = f;
        else
            error('BLOPEX:laplacian:InvalidClass',...
                '%s','Second input must be either class double or a cell array.');
        end
    else
        if dim == 1
            B = {'DD'};
        elseif dim == 2
            B = {'DD' 'DD'};
        else
            B = {'DD' 'DD' 'DD'};
        end
        m = 0;
    end

if max(size(m) - [1 1]) ~= 0
        error('BLOPEX:laplacian:WrongNumberOfEigenvalues',...
            '%s','The requested number of eigenvalues must be a scalar.');
    end

maxeigs = prod(u);
    mint = round(m);
    if mint ~= m || mint > maxeigs
        error('BLOPEX:laplacian:InvalidNumberOfEigs',...
            '%s','Number of eigenvalues output must be a nonnegative ',...
            'integer no bigger than number of grid points.');
    end
    m = mint;

bdryerr = 0;
    a = whos('regep','B');
    if sum(a.class(1:4)=='cell') ~= 4 || sum(a.size == [1 dim]) ~= 2
        bdryerr = 1;
    else
        BB = zeros(1, 2*dim);
        for i = 1:dim
            if (length(B{i}) == 1)
                if B{i} == 'P'
                    BB(i) = 3;
                    BB(i + dim) = 3;
                else
                    bdryerr = 1;
                end
            elseif (length(B{i}) == 2)
                if B{i}(1) == 'D'
                    BB(i) = 1;
                elseif B{i}(1) == 'N'
                    BB(i) = 2;
                else
                    bdryerr = 1;
                end
                if B{i}(2) == 'D'
                    BB(i + dim) = 1;
                elseif B{i}(2) == 'N'
                    BB(i + dim) = 2;
                else
                    bdryerr = 1;
                end
            else
                bdryerr = 1;
            end
        end
    end

if bdryerr == 1
        error('BLOPEX:laplacian:InvalidBdryConds',...
            '%s','Boundary conditions must be a cell of length 3 for 3D, 2',...
            ' for 2D, 1 for 1D, with values ''DD'', ''DN'', ''ND'', ''NN''',...
            ', or ''P''.');
    end

% Set the component matrices. SPDIAGS converts int8 into double anyway.
    e1 = ones(u(1),1); %e1 = ones(u(1),1,'int8');
    if dim > 1
        e2 = ones(u(2),1);
    end
    if dim > 2
        e3 = ones(u(3),1);
    end

% Calculate m smallest exact eigenvalues.
    if m > 0
        if (BB(1) == 1) && (BB(1+dim) == 1)
            a1 = pi/2/(u(1)+1);
            N = (1:u(1))';
        elseif (BB(1) == 2) && (BB(1+dim) == 2)
            a1 = pi/2/u(1);
            N = (0:(u(1)-1))';
        elseif ((BB(1) == 1) && (BB(1+dim) == 2)) || ((BB(1) == 2)...
                && (BB(1+dim) == 1))
            a1 = pi/4/(u(1)+0.5);
            N = 2*(1:u(1))' - 1;
        else
            a1 = pi/u(1);
            N = floor((1:u(1))/2)';
        end
        
        lambda1 = 4*sin(a1*N).^2;
        
        if dim > 1
            if (BB(2) == 1) && (BB(2+dim) == 1)
                a2 = pi/2/(u(2)+1);
                N = (1:u(2))';
            elseif (BB(2) == 2) && (BB(2+dim) == 2)
                a2 = pi/2/u(2);
                N = (0:(u(2)-1))';
            elseif ((BB(2) == 1) && (BB(2+dim) == 2)) || ((BB(2) == 2)...
                    && (BB(2+dim) == 1))
                a2 = pi/4/(u(2)+0.5);
                N = 2*(1:u(2))' - 1;
            else
                a2 = pi/u(2);
                N = floor((1:u(2))/2)';
            end
            lambda2 = 4*sin(a2*N).^2;
        else
            lambda2 = 0;
        end
        
        if dim > 2
            if (BB(3) == 1) && (BB(6) == 1)
                a3 = pi/2/(u(3)+1);
                N = (1:u(3))';
            elseif (BB(3) == 2) && (BB(6) == 2)
                a3 = pi/2/u(3);
                N = (0:(u(3)-1))';
            elseif ((BB(3) == 1) && (BB(6) == 2)) || ((BB(3) == 2)...
                    && (BB(6) == 1))
                a3 = pi/4/(u(3)+0.5);
                N = 2*(1:u(3))' - 1;
            else
                a3 = pi/u(3);
                N = floor((1:u(3))/2)';
            end
            lambda3 = 4*sin(a3*N).^2;
        else
            lambda3 = 0;
        end
        
        if dim == 1
            lambda = lambda1;
        elseif dim == 2
            lambda = kron(e2,lambda1) + kron(lambda2, e1);
        else
            lambda = kron(e3,kron(e2,lambda1)) + kron(e3,kron(lambda2,e1))...
                + kron(lambda3,kron(e2,e1));
        end
        [lambda, p] = sort(lambda);
        if m < maxeigs - 0.1
            w = lambda(m+1);
        else
            w = inf;
        end
        lambda = lambda(1:m);
        p = p(1:m)';
    else
        lambda = [];
    end

V = []; 
    if nargout > 1 && m > 0 % Calculate eigenvectors if specified in output.
        
        p1 = mod(p-1,u(1))+1;
        
        if (BB(1) == 1) && (BB(1+dim) == 1)
            V1 = sin(kron((1:u(1))'*(pi/(u(1)+1)),p1))*(2/(u(1)+1))^0.5;
        elseif (BB(1) == 2) && (BB(1+dim) == 2)
            V1 = cos(kron((0.5:1:u(1)-0.5)'*(pi/u(1)),p1-1))*(2/u(1))^0.5;
            V1(:,p1==1) = 1/u(1)^0.5;
        elseif ((BB(1) == 1) && (BB(1+dim) == 2))
            V1 = sin(kron((1:u(1))'*(pi/2/(u(1)+0.5)),2*p1 - 1))...
                *(2/(u(1)+0.5))^0.5;
        elseif ((BB(1) == 2) && (BB(1+dim) == 1))
            V1 = cos(kron((0.5:1:u(1)-0.5)'*(pi/2/(u(1)+0.5)),2*p1 - 1))...
                *(2/(u(1)+0.5))^0.5;
        else
            V1 = zeros(u(1),m);
            a = (0.5:1:u(1)-0.5)';
            V1(:,mod(p1,2)==1) = cos(a*(pi/u(1)*(p1(mod(p1,2)==1)-1)))...
                *(2/u(1))^0.5;
            pp = p1(mod(p1,2)==0);
            if ~isempty(pp)
                V1(:,mod(p1,2)==0) = sin(a*(pi/u(1)*p1(mod(p1,2)==0)))...
                    *(2/u(1))^0.5;
            end
            V1(:,p1==1) = 1/u(1)^0.5;
            if mod(u(1),2) == 0
                V1(:,p1==u(1)) = V1(:,p1==u(1))/2^0.5;
            end
        end
        
        
        if dim > 1
            p2 = mod(p-p1,u(1)*u(2));
            p3 = (p - p2 - p1)/(u(1)*u(2)) + 1;
            p2 = p2/u(1) + 1;
            if (BB(2) == 1) && (BB(2+dim) == 1)
                V2 = sin(kron((1:u(2))'*(pi/(u(2)+1)),p2))*(2/(u(2)+1))^0.5;
            elseif (BB(2) == 2) && (BB(2+dim) == 2)
                V2 = cos(kron((0.5:1:u(2)-0.5)'*(pi/u(2)),p2-1))*(2/u(2))^0.5;
                V2(:,p2==1) = 1/u(2)^0.5;
            elseif ((BB(2) == 1) && (BB(2+dim) == 2))
                V2 = sin(kron((1:u(2))'*(pi/2/(u(2)+0.5)),2*p2 - 1))...
                    *(2/(u(2)+0.5))^0.5;
            elseif ((BB(2) == 2) && (BB(2+dim) == 1))
                V2 = cos(kron((0.5:1:u(2)-0.5)'*(pi/2/(u(2)+0.5)),2*p2 - 1))...
                    *(2/(u(2)+0.5))^0.5;
            else
                V2 = zeros(u(2),m);
                a = (0.5:1:u(2)-0.5)';
                V2(:,mod(p2,2)==1) = cos(a*(pi/u(2)*(p2(mod(p2,2)==1)-1)))...
                    *(2/u(2))^0.5;
                pp = p2(mod(p2,2)==0);
                if ~isempty(pp)
                    V2(:,mod(p2,2)==0) = sin(a*(pi/u(2)*p2(mod(p2,2)==0)))...
                        *(2/u(2))^0.5;
                end
                V2(:,p2==1) = 1/u(2)^0.5;
                if mod(u(2),2) == 0
                    V2(:,p2==u(2)) = V2(:,p2==u(2))/2^0.5;
                end
            end
        else
            V2 = ones(1,m);
        end
        
        if dim > 2
            if (BB(3) == 1) && (BB(3+dim) == 1)
                V3 = sin(kron((1:u(3))'*(pi/(u(3)+1)),p3))*(2/(u(3)+1))^0.5;
            elseif (BB(3) == 2) && (BB(3+dim) == 2)
                V3 = cos(kron((0.5:1:u(3)-0.5)'*(pi/u(3)),p3-1))*(2/u(3))^0.5;
                V3(:,p3==1) = 1/u(3)^0.5;
            elseif ((BB(3) == 1) && (BB(3+dim) == 2))
                V3 = sin(kron((1:u(3))'*(pi/2/(u(3)+0.5)),2*p3 - 1))...
                    *(2/(u(3)+0.5))^0.5;
            elseif ((BB(3) == 2) && (BB(3+dim) == 1))
                V3 = cos(kron((0.5:1:u(3)-0.5)'*(pi/2/(u(3)+0.5)),2*p3 - 1))...
                    *(2/(u(3)+0.5))^0.5;
            else
                V3 = zeros(u(3),m);
                a = (0.5:1:u(3)-0.5)';
                V3(:,mod(p3,2)==1) = cos(a*(pi/u(3)*(p3(mod(p3,2)==1)-1)))...
                    *(2/u(3))^0.5;
                pp = p1(mod(p3,2)==0);
                if ~isempty(pp)
                    V3(:,mod(p3,2)==0) = sin(a*(pi/u(3)*p3(mod(p3,2)==0)))...
                        *(2/u(3))^0.5;
                end
                V3(:,p3==1) = 1/u(3)^0.5;
                if mod(u(3),2) == 0
                    V3(:,p3==u(3)) = V3(:,p3==u(3))/2^0.5;
                end
                
            end
        else
            V3 = ones(1,m);
        end
        
        if dim == 1
            V = V1;
        elseif dim == 2
            V = kron(e2,V1).*kron(V2,e1);
        else
            V = kron(e3, kron(e2, V1)).*kron(e3, kron(V2, e1))...
                .*kron(kron(V3,e2),e1);
        end
        
        if m ~= 0
            if abs(lambda(m) - w) < maxeigs*eps('double')
                sprintf('\n%s','Warning: (m+1)th eigenvalue is  nearly equal',...
                    ' to mth.')
                
            end
        end
        
    end

A = [];
    if nargout > 2 %also calulate the matrix if specified in the output
        
        % Set the component matrices. SPDIAGS converts int8 into double anyway.
        %    e1 = ones(u(1),1); %e1 = ones(u(1),1,'int8');
        D1x = spdiags([-e1 2*e1 -e1], [-1 0 1], u(1),u(1));
        if dim > 1
            %        e2 = ones(u(2),1);
            D1y = spdiags([-e2 2*e2 -e2], [-1 0 1], u(2),u(2));
        end
        if dim > 2
            %        e3 = ones(u(3),1);
            D1z = spdiags([-e3 2*e3 -e3], [-1 0 1], u(3),u(3));
        end
        
        
        % Set boundary conditions if other than Dirichlet.
        for i = 1:dim
            if BB(i) == 2
                eval(['D1' char(119 + i) '(1,1) = 1;'])
            elseif BB(i) == 3
                eval(['D1' char(119 + i) '(1,' num2str(u(i)) ') = D1'...
                    char(119 + i) '(1,' num2str(u(i)) ') -1;']);
                eval(['D1' char(119 + i) '(' num2str(u(i)) ',1) = D1'...
                    char(119 + i) '(' num2str(u(i)) ',1) -1;']);
            end
            
            if BB(i+dim) == 2
                eval(['D1' char(119 + i)...
                    '(',num2str(u(i)),',',num2str(u(i)),') = 1;'])
            end
        end
        
        % Form A using tensor products of lower dimensional Laplacians
        Ix = speye(u(1));
        if dim == 1
            A = D1x;
        elseif dim == 2
            Iy = speye(u(2));
            A = kron(Iy,D1x) + kron(D1y,Ix);
        elseif dim == 3
            Iy = speye(u(2));
            Iz = speye(u(3));
            A = kron(Iz, kron(Iy, D1x)) + kron(Iz, kron(D1y, Ix))...
                + kron(kron(D1z,Iy),Ix);
        end
    end

disp('  ')
    if ~isempty(V)
        a = whos('regep','V');
        disp(['The eigenvectors take ' num2str(a.bytes) ' bytes'])
    end
    if  ~isempty(A)
        a = whos('regexp','A');
        disp(['The Laplacian matrix takes ' num2str(a.bytes) ' bytes'])
    end
    disp('  ')
endfunction
clf;
% A=delsq(numgrid('S',m));
[~,~,A]=laplacian(30,30);
%[lamda,A] = laplacian([50,20,20],{'DD' 'NN' 'P'});

spy(A);
print -dpng name.png;

[ + ] Show input
Absolute running time: 1.02 sec, cpu time: 1.67 sec, memory peak: 39 Mb, absolute service time: 1,07 sec