Does julia has a Package designed to / for sparse matrix

I mean in solving Numerical PDE the matrix is usually quiet big but sparse. the normal “/” seems to use LU decomposition which makes a saprse matrix become two dense and will outofmemory qiuckly.

If you have a sparse matrix, then \ uses UMFPACK for LU decomposition, so it should not make it dense. Do you have an example of this happening?

2 Likes

This is the code.Most part is building a sparse matrix. It is solving a three dimensional PDE and N is the how much part we divide the PDE in different dimension. When N=30 it warks well, but when N=100

Stacktrace:
 [1] umfpack_numeric!(::SuiteSparse.UMFPACK.UmfpackLU{Float64,Int64}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\SuiteSparse\src\umfpack.jl:258
 [2] #lu#1(::Bool, ::typeof(lu), ::SparseMatrixCSC{Float64,Int64}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\SuiteSparse\src\umfpack.jl:160
 [3] lu at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\SuiteSparse\src\umfpack.jl:154 [inlined]
 [4] \(::SparseMatrixCSC{Float64,Int64}, ::Array{Float64,2}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\SparseArrays\src\linalg.jl:1374
 [5] pdesolver(::Int64, ::Float64, ::Float64, ::Float64, ::Array{Float64,2}, ::Int64) at C:\Users\lyu\Documents\DC\dachuang-again\threeDimensionalProblem.jl:349
 [6] pdesolver() at C:\Users\lyu\Documents\DC\dachuang-again\threeDimensionalProblem.jl:13
 [7] top-level scope at REPL[7]:1
using SparseArrays, PyPlot,NPZ
using Flux, Flux.Data.MNIST, Statistics
using Flux: onehotbatch, onehot, crossentropy, throttle,Params
using Base.Iterators: repeated
using LinearAlgebra
using Sobol
using BSON:@save
using BSON:@load

using FastGaussQuadrature,SparseGrids
## K and L is used to define the schostic layer
function pdesolver(σ=1,d=pi/2,L_y=2π,L_z=2π,ϵ=zeros(3,3),λ=-1)
    N=30
    K=size(ϵ,1)
    L=size(ϵ,2)
    dx=1/N;dy=1/N;dz=1/N
    x=0:dx:1;y=0:dx:1;z=0:dx:1
    Nx=length(x);Ny=length(y);Nz=length(z)
    # define the stochastic layer here
    h(rr1,rr2,sy,sz)    =   (                     sin(2*rr1*pi*sy)*sin(2*rr2*pi*sz));
    h_y(rr1,rr2,sy,sz)  =   (0   +  2*rr1*pi*     cos(2*rr1*pi*sy)*sin(2*rr2*pi*sz));
    h_yy(rr1,rr2,sy,sz) =   (0   -  (2*rr2*pi)^2* sin(2*rr1*pi*sy)*sin(2*rr2*pi*sz));
    h_z(rr1,rr2,sy,sz)  =   (0   +  2*rr2*pi*     sin(2*rr1*pi*sy)*cos(2*rr2*pi*sz));
    h_zz(rr1,rr2,sy,sz) =   (0   -  (2*rr2*pi)^2* sin(2*rr1*pi*sy)*sin(2*rr2*pi*sz));
    ##The main finite difference code for u, u1, and u2
    EPh =0;EU= zeros(Nx,Ny,Nz);
    Datah = zeros(Ny,Nz);   Datah_y = zeros(Ny,Nz);  Datah_yy = zeros(Ny,Nz);   Datah_z = zeros(Ny,Nz); Datah_zz = zeros(Ny,Nz);
    sumindex=0   # size: Nz*1
    for rr1=1:K
        for rr2=1:L
            multiplier=float(rr1)^(λ)*float(rr2)^(λ)
            Datah     +=  h.(rr1,rr2,y,z')*ϵ[rr1,rr2]*multiplier;
            Datah_y   +=  h_y.(rr1,rr2,y,z')*ϵ[rr1,rr2]*multiplier;
            Datah_yy  +=  h_yy.(rr1,rr2,y,z')*ϵ[rr1,rr2]*multiplier;
            Datah_z   +=  h_z.(rr1,rr2,y,z')*ϵ[rr1,rr2]*multiplier;
            Datah_zz  +=  h_zz.(rr1,rr2,y,z')*ϵ[rr1,rr2]*multiplier;
            sumindex  +=  multiplier
        end
    end
    Datah=Datah/sumindex
    Datah_y=Datah_y/sumindex
    Datah_yy=Datah_yy/sumindex
    Datah_z=Datah_z/sumindex
    Datah_zz=Datah_zz/sumindex
    ## conputing, loop over each experimental result
    ## form the data structure and index for stiffness matrix and load vector
    dx2=dx*dx; dy2=dy*dy; dz2=dz*dz; d4xy=4*dx*dy; d4xz=4*dx*dz; σ2=σ*σ;
    DOFsp=(Nx-3)*Ny*Nz*15+Ny*Nz*10*2;
    DOFu=(Nx-1)*Ny*Nz;
    spI4Au=zeros(DOFsp,1); spJ4Au=zeros(DOFsp,1); spS4Au=zeros(DOFsp,1);
    spI4bu=zeros(DOFu,1);  spJ4bu=ones(DOFu,1);
    spS4bu = zeros(DOFu,1);  spS4bu1 = zeros(DOFu,1);  spS4bu2 = zeros(DOFu,1);
    ## define the weights for trapezoidal rule
    trapz  = 2*ones(Nz,1);  trapy  = 2*ones(1,Ny); trapx= 2*ones(1,Nx);trapwt=zeros(Nx,Ny,Nz)  # define the weights for trapezoidal rule
    trapz[1] = 1; trapz[Nz] = 1;
    trapy[1] = 1; trapy[Ny] = 1;
    trapx[1] = 1; trapx[Nx] = 1;
    for i=1:Nx
        for j=1:Ny
            for k=1:Nz
            trapwt[i,j,k]=trapx[i]*trapy[j]*trapz[k];
            end
        end
    end
    wt=reshape(trapwt,DOFu+Ny*Nz,1);
    ##range of index for u. For x, 2<=i<=N_x; For y, 1<=j<=N_y; For z, 1<=k<=N_z
    index=0; index4b=0;
    for i=3:Nx-1
        for j=2:Ny-1
            for k=2:Nz-1
            spI4Au[index+1:index+15,1]=    ones(15)*((i  -2)*Ny*Nz+(j  -1)*Nz+k);
            spJ4Au[index+1:index+15,1]= [  (i-1-2)*Ny*Nz+(j-1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k; (i-1-2)*Nz*Ny+(j+1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k-1; (i-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                           (i  -2)*Ny*Nz+(j-1-1)*Nz+k; (i  -2)*Ny*Nz+(j-1)*Nz+k; (i  -2)*Nz*Ny+(j+1-1)*Nz+k; (i  -2)*Ny*Nz+(j-1)*Nz+k-1; (i  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                           (i+1-2)*Ny*Nz+(j-1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k; (i+1-2)*Nz*Ny+(j+1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k-1; (i+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
            spI4bu[index4b+1 ,1]= (i  -2)*Ny*Nz+(j  -1)*Nz+k;
            index=index+15; index4b=index4b+1;
            end
        end
    end
    for j=1:Ny   #front grids,i=2
        for k=1:Nz
            spI4Au[index+1:index+10,1]=ones(10)*((2-2)*Ny*Nz+(j-1)*Nz+k);
            if j==1
                if k==1
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(Ny-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(j+1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (2  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                               (2+1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
                elseif k==Nz
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(Ny-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(j+1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k-1; (2  -2)*Ny*Nz+(j-1)*Nz+2;
                                            (2+1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+2;];
                else
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(Ny-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(j+1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k-1; (2  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (2+1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
                end

            elseif j==Ny
                if k==1
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(j-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(2-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (2  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(2-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
                elseif k==Nz
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(j-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(2-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k-1; (2  -2)*Ny*Nz+(j-1)*Nz+2;
                                            (2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(2-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+2;];
                else
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(j-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(2-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k-1; (2  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(2-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
                end
            else
                if k==1
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(j-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(j+1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (2  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
                elseif k==Nz
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(j-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(j+1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k-1; (2  -2)*Ny*Nz+(j-1)*Nz+2;
                                            (2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+2;];
                else
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(j-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(j+1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k-1; (2  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
                end
            end
            spI4bu[index4b+1,1]=(2-2)*Ny*Nz+(j-1)*Nz+k;
            index=index+10; index4b=index4b+1;
        end
    end
    for j=1:Ny    #back grids
        for k=1:Nz
            if j==1
                if k==1
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (Nx  -2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+k+1;];
                elseif k==Nz
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+2;
                                            (Nx  -2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+2;];
                else
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (Nx  -2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+k+1;];
                end
            elseif j==Ny
                if k==1
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(2-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (Nx  -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(2-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+k+1;];
                elseif k==Nz
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(2-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+2;
                                            (Nx  -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(2-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+2;];
                else
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(2-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (Nx  -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(2-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+k+1;];
                end
            else
                if k==1
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (Nx  -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+k+1;];
                elseif k==Nz
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+2;
                                            (Nx  -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+2;];
                else
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (Nx  -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+k+1;];
                end
            end
            spI4Au[index+1:index+10,1]=ones(10)*((Nx-2)*Ny*Nz+(j-1)*Nz+k);
            spI4bu[index4b+1,1]=(Nx-2)*Ny*Nz+(j-1)*Nz+k;
            index=index+10; index4b=index4b+1;
        end
    end
    for i=3:Nx-1 #rightest grids
        for j=1:Ny #k=Nz
            if j==1
                spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(Ny-1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz; (i-1-2)*Nz*Ny+(j+1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i  -2)*Ny*Nz+(Ny-1-1)*Nz+Nz; (i  -2)*Ny*Nz+(j-1)*Nz+Nz; (i  -2)*Nz*Ny+(j+1-1)*Nz+Nz; (i  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i  -2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i+1-2)*Ny*Nz+(Ny-1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz; (i+1-2)*Nz*Ny+(j+1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
            elseif j==Ny
                spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz; (i-1-2)*Nz*Ny+(1+1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i  -2)*Ny*Nz+(j-1-1)*Nz+Nz; (i  -2)*Ny*Nz+(j-1)*Nz+Nz; (i  -2)*Nz*Ny+(1+1-1)*Nz+Nz; (i  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i  -2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i+1-2)*Ny*Nz+(j-1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz; (i+1-2)*Nz*Ny+(1+1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
            else
                spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz; (i-1-2)*Nz*Ny+(j+1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i  -2)*Ny*Nz+(j-1-1)*Nz+Nz; (i  -2)*Ny*Nz+(j-1)*Nz+Nz; (i  -2)*Nz*Ny+(j+1-1)*Nz+Nz; (i  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i  -2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i+1-2)*Ny*Nz+(j-1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz; (i+1-2)*Nz*Ny+(j+1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
            end
            spI4Au[index+1:index+15]=((i-2)*Ny*Nz+(j-1)*Nz+Nz)*ones(15);
            spI4bu[index4b+1,1]=(i-2)*Ny*Nz+(j-1)*Nz+Nz;
            index=index+15; index4b=index4b+1;
        end
    end
    for i=3:Nx-1 #leftest grids
        for j=1:Ny #k=1
            if j==1
            spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(Ny-1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+1; (i-1-2)*Nz*Ny+(j+1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i  -2)*Ny*Nz+(Ny-1-1)*Nz+1; (i  -2)*Ny*Nz+(j-1)*Nz+1; (i  -2)*Nz*Ny+(j+1-1)*Nz+1; (i  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i  -2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i+1-2)*Ny*Nz+(Ny-1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+1; (i+1-2)*Nz*Ny+(j+1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
            elseif j==Ny
            spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+1; (i-1-2)*Nz*Ny+(1+1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i  -2)*Ny*Nz+(j-1-1)*Nz+1; (i  -2)*Ny*Nz+(j-1)*Nz+1; (i  -2)*Nz*Ny+(1+1-1)*Nz+1; (i  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i  -2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i+1-2)*Ny*Nz+(j-1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+1; (i+1-2)*Nz*Ny+(1+1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
            else
            spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+1; (i-1-2)*Nz*Ny+(j+1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i  -2)*Ny*Nz+(j-1-1)*Nz+1; (i  -2)*Ny*Nz+(j-1)*Nz+1; (i  -2)*Nz*Ny+(j+1-1)*Nz+1; (i  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i  -2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i+1-2)*Ny*Nz+(j-1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+1; (i+1-2)*Nz*Ny+(j+1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
            end
            spI4Au[index+1:index+15]=ones(15)*((i-2)*Ny*Nz+(j-1)*Nz+1);
            spI4bu[index4b+1,1]=(i-2)*Ny*Nz+(j-1)*Nz+1;
            index=index+15; index4b=index4b+1;
        end
    end
    for i=3:Nx-1 #upper grids j=1
        j=1;
        for k=2:Nz-1
            spI4Au[index+1:index+15,1]=((i-2)*Ny*Nz+(1-1)*Nz+k)*ones(15);
            spJ4Au[index+1:index+15,1]=[(i-1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k; (i-1-2)*Nz*Ny+(j+1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k-1; (i-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                    (i  -2)*Ny*Nz+(Ny-1-1)*Nz+k; (i  -2)*Ny*Nz+(j-1)*Nz+k; (i  -2)*Nz*Ny+(j+1-1)*Nz+k; (i  -2)*Ny*Nz+(j-1)*Nz+k-1; (i  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                    (i+1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k; (i+1-2)*Nz*Ny+(j+1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k-1; (i+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
            spI4bu[index4b+1,1]= (i-2)*Ny*Nz+(j-1)*Nz+k;
            index=index+15; index4b=index4b+1;
        end
    end
    for i=3:Nx-1 #upper grids j=Ny
        j=Ny;
        for k=2:Nz-1
            spI4Au[index+1:index+15,1]=ones(15)*((i-2)*Ny*Nz+(j-1)*Nz+k);
            spJ4Au[index+1:index+15,1]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k; (i-1-2)*Nz*Ny+(1+1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k-1; (i-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                    (i  -2)*Ny*Nz+(j-1-1)*Nz+k; (i  -2)*Ny*Nz+(j-1)*Nz+k; (i  -2)*Nz*Ny+(1+1-1)*Nz+k; (i  -2)*Ny*Nz+(j-1)*Nz+k-1; (i  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                    (i+1-2)*Ny*Nz+(j-1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k; (i+1-2)*Nz*Ny+(1+1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k-1; (i+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
            spI4bu[index4b+1,1]= (i-2)*Ny*Nz+(j-1)*Nz+k;
            index=index+15; index4b=index4b+1;
        end
    end
    DataG=zeros(Nx,Ny,Nz)
    for i=1:Nx
        for j=1:Ny
            for k=1:Nz
                DataG[i,j,k]=-exp(-(d-Datah[j,k])*(1-x[i]))#-exp(-((d-Datah[j,k])*(1-x[i])).^2)
            end
        end
    end
    #-2sin(d*x[i])*(sin(L_y*y[j])+1)*(sin(L_z*z[k])+1)-sin(d*x[i])*sin(L_y*y[j])*(sin(L_z*z[k])+1)-sin(d*x[i])*(sin(L_y*y[j])+1)*sin(L_z*z[k])
    ## assmble and solve for u
    index=0;index4b=0;
    for i=3:Nx-1
        for j=2:Ny-1
            for k=2:Nz-1
                A=σ2*((1-x[i])^2*(Datah_y[j,k]^2+Datah_z[j,k]^2)+1)/(d-Datah[j,k])^2+ σ2/(d-Datah[j,k])^2;
                B=σ2/L_y^2;
                C=σ2/L_z^2;
                D=-2σ2*(1/L_y)*(1-x[i])*Datah_y[j,k]/(d-Datah[j,k]);
                E=-2σ2*(1/L_z)*(1-x[i])*Datah_z[j,k]/(d-Datah[j,k]);
                F=-2*σ2*((1-x[i])*(Datah_y[j,k]^2+Datah_z[j,k]^2)/(d-Datah[j,k])^2+0.5*(1-x[i])*(Datah_yy[j,k]+Datah_zz[j,k])/(d-Datah[j,k]));
                spS4Au[index+1:index+15,1]=[D/d4xy;   A/dx2-F/(2*dx);             -D/d4xy;   E/d4xz; -E/d4xz;
                                        B/dy2;    -2*A/dx2-2*B/dy2-2*C/dz2-1; B/dy2;     C/dz2;   C/dz2;
                                       -D/d4xy;   A/dx2+F/(2*dx);             D/d4xy;    -E/d4xz; E/d4xz;];
                spS4bu[index4b+1,1]=DataG[i,j,k]#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
                index=index+15; index4b=index4b+1;
            end
        end
    end
    i=2;
    for j=1:Ny   #front grids,i=2
        for k=1:Nz
            A=σ2*((1-x[i])^2*(Datah_y[j,k]^2+Datah_z[j,k]^2)+1)/(d-Datah[j,k])^2+ σ2/(d-Datah[j,k])^2;
            B=σ2/L_y^2;
            C=σ2/L_z^2;
            D=-2σ2*(1/L_y)*(1-x[i])*Datah_y[j,k]/(d-Datah[j,k]);
            E=-2σ2*(1L_z)*(1-x[i])*Datah_z[j,k]/(d-Datah[j,k]);
            F=-2*σ2*((1-x[i])*(Datah_y[j,k]^2+Datah_z[j,k]^2)/(d-Datah[j,k])^2+0.5*(1-x[i])*(Datah_yy[j,k]+Datah_zz[j,k])/(d-Datah[j,k]));
            spS4Au[index+1:index+10,1]=[B/dy2; -2*A/dx2-2*B/dy2-2*C/dz2-1; B/dy2; C/dz2;    C/dz2;
            -D/d4xy;A/dx2+F/(2*dx);            D/d4xy; -E/d4xz; E/d4xz;];
            spS4bu[index4b+1,1]=DataG[i,j,k]#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
            index=index+10; index4b=index4b+1;
        end
    end
    i=Nx;
    for j=1:Ny    #back grids
        for k=1:Nz
            A=σ2*((1-x[i])^2*(Datah_y[j,k]^2+Datah_z[j,k]^2)+1)/(d-Datah[j,k])^2+ σ2/(d-Datah[j,k])^2;
            B=σ2/L_y^2;
            C=σ2/L_z^2;
            D=-σ2*(2/L_y)*(1-x[i])*Datah_y[j,k]/(d-Datah[j,k]);
            E=-σ2*(2/L_z)*(1-x[i])*Datah_z[j,k]/(d-Datah[j,k]);
            F=-2*σ2*((1-x[i])*(Datah_y[j,k]^2+Datah_z[j,k]^2)/(d-Datah[j,k])^2+0.5*(1-x[i])*(Datah_yy[j,k]+Datah_zz[j,k])/(d-Datah[j,k]));
            spS4Au[index+1:index+10,1]=[0;      2*A/dx2; 0;0; 0;
                B/dy2; -2*A/dx2-2*B/dy2-2*C/dz2-1; B/dy2; C/dz2; C/dz2;];
            spS4bu[index4b+1,1]=DataG[i,j,k]#-#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
            index=index+10; index4b=index4b+1;
        end
    end
    k=Nz;
    for i=3:Nx-1 #rightest grids
        for j=1:Ny #k=Nz
            A=σ2*((1-x[i])^2*(Datah_y[j,k]^2+Datah_z[j,k]^2)+1)/(d-Datah[j,k])^2+ σ2/(d-Datah[j,k])^2;
            B=σ2/L_y^2;
            C=σ2/L_z^2;
            D=-σ2*(2/L_y)*(1-x[i])*Datah_y[j,k]/(d-Datah[j,k]);
            E=-σ2*(2/L_z)*(1-x[i])*Datah_z[j,k]/(d-Datah[j,k]);
            F=-2*σ2*((1-x[i])*(Datah_y[j,k]^2+Datah_z[j,k]^2)/(d-Datah[j,k])^2+0.5*(1-x[i])*(Datah_yy[j,k]+Datah_zz[j,k])/(d-Datah[j,k]));
            spS4Au[index+1:index+15,1]=[D/d4xy;   A/dx2-F/(2*dx);             -D/d4xy;   E/d4xz; -E/d4xz;
                B/dy2;    -2*A/dx2-2*B/dy2-2*C/dz2-1; B/dy2;     C/dz2;   C/dz2;
                -D/d4xy;   A/dx2+F/(2*dx);             D/d4xy;    -E/d4xz; E/d4xz;];
            spS4bu[index4b+1,1]=DataG[i,j,k]#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
            index=index+15; index4b=index4b+1;
        end
    end
    for i=3:Nx-1 #leftest grids
        k=1;
        for j=1:Ny #k=1
                A=σ2*((1-x[i])^2*(Datah_y[j,k]^2+Datah_z[j,k]^2)+1)/(d-Datah[j,k])^2+ σ2/(d-Datah[j,k])^2;
                B=σ2/L_y^2;
                C=σ2/L_z^2;
                D=-σ2*(2/L_y)*(1-x[i])*Datah_y[j,k]/(d-Datah[j,k]);
                E=-σ2*(2/L_z)*(1-x[i])*Datah_z[j,k]/(d-Datah[j,k]);
                F=-2*σ2*((1-x[i])*(Datah_y[j,k]^2+Datah_z[j,k]^2)/(d-Datah[j,k])^2+0.5*(1-x[i])*(Datah_yy[j,k]+Datah_zz[j,k])/(d-Datah[j,k]));
                spS4Au[index+1:index+15,1]=[D/d4xy;   A/dx2-F/(2*dx);             -D/d4xy;   E/d4xz; -E/d4xz;
                B/dy2;    -2*A/dx2-2*B/dy2-2*C/dz2-1; B/dy2;     C/dz2;   C/dz2;
                -D/d4xy;   A/dx2+F/(2*dx);             D/d4xy;    -E/d4xz; E/d4xz;];
            spS4bu[index4b+1,1]=DataG[i,j,k]#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
            index=index+15; index4b=index4b+1;
        end
    end
    for i=3:Nx-1 #upper grids j=1
    j=1;
        for k=2:Nz-1
            A=σ2*((1-x[i])^2*(Datah_y[j,k]^2+Datah_z[j,k]^2)+1)/(d-Datah[j,k])^2+ σ2/(d-Datah[j,k])^2;
            B=σ2/L_y^2;
            C=σ2/L_z^2;
            D=-σ2*(2/L_y)*(1-x[i])*Datah_y[j,k]/(d-Datah[j,k]);
            E=-σ2*(2/L_z)*(1-x[i])*Datah_z[j,k]/(d-Datah[j,k]);
            F=-2*σ2*((1-x[i])*(Datah_y[j,k]^2+Datah_z[j,k]^2)/(d-Datah[j,k])^2+0.5*(1-x[i])*(Datah_yy[j,k]+Datah_zz[j,k])/(d-Datah[j,k]));
            spS4Au[index+1:index+15,1]=[D/d4xy; A/dx2-F/(2*dx); -D/d4xy; E/d4xz; -E/d4xz;
                B/dy2; -2*A/dx2-2*B/dy2-2*C/dz2-1; B/dy2; C/dz2; C/dz2;
                -D/d4xy;A/dx2+F/(2*dx); D/4*d4xy; -E/d4xz; E/d4xz;];
            spS4bu[index4b+1,1]=DataG[i,j,k]#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
            index=index+15; index4b=index4b+1;
        end
    end
    for i=3:Nx-1 #upper grids j=Ny
        j=Ny;
        for k=2:Nz-1
            A=σ2*((1-x[i])^2*(Datah_y[j,k]^2+Datah_z[j,k]^2)+1)/(d-Datah[j,k])^2+ σ2/(d-Datah[j,k])^2;
            B=σ2/L_y^2;
            C=σ2/L_z^2;
            D=-σ2*(2/L_y)*(1-x[i])*Datah_y[j,k]/(d-Datah[j,k]);
            E=-σ2*(2/L_z)*(1-x[i])*Datah_z[j,k]/(d-Datah[j,k]);
            F=-2*σ2*((1-x[i])*(Datah_y[j,k]^2+Datah_z[j,k]^2)/(d-Datah[j,k])^2+0.5*(1-x[i])*(Datah_yy[j,k]+Datah_zz[j,k])/(d-Datah[j,k]));
            spS4Au[index+1:index+15,1]=[D/d4xy;   A/dx2-F/(2*dx);             -D/d4xy;   E/d4xz; -E/d4xz;
                B/dy2;    -2*A/dx2-2*B/dy2-2*C/dz2-1; B/dy2;     C/dz2;   C/dz2;
                -D/d4xy;   A/dx2+F/(2*dx);             D/d4xy;    -E/d4xz; E/d4xz;];
            spS4bu[index4b+1,1]=DataG[i,j,k]#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
            index=index+15; index4b=index4b+1;
        end
    end
    Au =sparse(vec(spI4Au),vec(spJ4Au),vec(spS4Au),DOFu,DOFu);
    bu = sparse(vec(spI4bu),vec(spJ4bu),vec(spS4bu),DOFu,1);
    u=Au\Matrix(bu)
    u= [zeros(Ny*Nz,1); u];
    U=reshape(Matrix(u),Nx,Ny,Nz);
    Phtmp = (1/8)*sum(sum(sum((d*ones(Nx,Nx^2)-repeat(Datah,1,Nx)).*reshape(U,Nx,Nx^2).*reshape(trapwt,Nx,Nx^2)*dx*dy*dz)))
    return  Phtmp
end

A minimal example would be kind on your helpers. Have a read through PSA: make it easier to help you, in particular point 4.

3 Likes

I understand your meaning, but the problem is building a sparse matrix is the main part of the code. In order to perserve the feature of this matrix, I cannot change most part of the code. Only thing I can do is change the value in each index into a rand number like under, but it may also cause some SingularException(0) problems because of the / sometimes.

Actually, the only useful line is the last five line u = Au\bu, when the size of the matric increase it will Outofmemory

function pdesolver(σ=1,d=pi/2,L_y=2π,L_z=2π,ϵ=zeros(3,3),λ=-1)
    N=20
    dx=1/N;dy=1/N;dz=1/N
    x=0:dx:1;y=0:dx:1;z=0:dx:1
    Nx=length(x);Ny=length(y);Nz=length(z)
    ##The main finite difference code for u, u1, and u2
    EPh =0;EU= zeros(Nx,Ny,Nz);
    Datah = zeros(Ny,Nz);   Datah_y = zeros(Ny,Nz);  Datah_yy = zeros(Ny,Nz);   Datah_z = zeros(Ny,Nz); Datah_zz = zeros(Ny,Nz);
    ## conputing, loop over each experimental result
    ## form the data structure and index for stiffness matrix and load vector
    dx2=dx*dx; dy2=dy*dy; dz2=dz*dz; d4xy=4*dx*dy; d4xz=4*dx*dz; σ2=σ*σ;
    DOFsp=(Nx-3)*Ny*Nz*15+Ny*Nz*10*2;
    DOFu=(Nx-1)*Ny*Nz;
    spI4Au=zeros(DOFsp,1); spJ4Au=zeros(DOFsp,1); spS4Au=zeros(DOFsp,1);
    spI4bu=zeros(DOFu,1);  spJ4bu=ones(DOFu,1);
    spS4bu = zeros(DOFu,1);  spS4bu1 = zeros(DOFu,1);  spS4bu2 = zeros(DOFu,1);
    ##range of index for u. For x, 2<=i<=N_x; For y, 1<=j<=N_y; For z, 1<=k<=N_z
    index=0; index4b=0;
    for i=3:Nx-1
        for j=2:Ny-1
            for k=2:Nz-1
            spI4Au[index+1:index+15,1]=    ones(15)*((i  -2)*Ny*Nz+(j  -1)*Nz+k);
            spJ4Au[index+1:index+15,1]= [  (i-1-2)*Ny*Nz+(j-1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k; (i-1-2)*Nz*Ny+(j+1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k-1; (i-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                           (i  -2)*Ny*Nz+(j-1-1)*Nz+k; (i  -2)*Ny*Nz+(j-1)*Nz+k; (i  -2)*Nz*Ny+(j+1-1)*Nz+k; (i  -2)*Ny*Nz+(j-1)*Nz+k-1; (i  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                           (i+1-2)*Ny*Nz+(j-1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k; (i+1-2)*Nz*Ny+(j+1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k-1; (i+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
            spI4bu[index4b+1 ,1]= (i  -2)*Ny*Nz+(j  -1)*Nz+k;
            index=index+15; index4b=index4b+1;
            end
        end
    end
    for j=1:Ny   #front grids,i=2
        for k=1:Nz
            spI4Au[index+1:index+10,1]=ones(10)*((2-2)*Ny*Nz+(j-1)*Nz+k);
            if j==1
                if k==1
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(Ny-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(j+1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (2  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                               (2+1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
                elseif k==Nz
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(Ny-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(j+1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k-1; (2  -2)*Ny*Nz+(j-1)*Nz+2;
                                            (2+1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+2;];
                else
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(Ny-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(j+1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k-1; (2  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (2+1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
                end

            elseif j==Ny
                if k==1
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(j-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(2-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (2  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(2-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
                elseif k==Nz
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(j-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(2-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k-1; (2  -2)*Ny*Nz+(j-1)*Nz+2;
                                            (2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(2-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+2;];
                else
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(j-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(2-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k-1; (2  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(2-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
                end
            else
                if k==1
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(j-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(j+1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (2  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
                elseif k==Nz
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(j-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(j+1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k-1; (2  -2)*Ny*Nz+(j-1)*Nz+2;
                                            (2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+2;];
                else
                    spJ4Au[index+1:index+10,1]=[(2  -2)*Ny*Nz+(j-1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k; (2  -2)*Nz*Ny+(j+1-1)*Nz+k; (2  -2)*Ny*Nz+(j-1)*Nz+k-1; (2  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (2+1-2)*Ny*Nz+(j-1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k; (2+1-2)*Nz*Ny+(j+1-1)*Nz+k; (2+1-2)*Ny*Nz+(j-1)*Nz+k-1; (2+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
                end
            end
            spI4bu[index4b+1,1]=(2-2)*Ny*Nz+(j-1)*Nz+k;
            index=index+10; index4b=index4b+1;
        end
    end
    for j=1:Ny    #back grids
        for k=1:Nz
            if j==1
                if k==1
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (Nx  -2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+k+1;];
                elseif k==Nz
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+2;
                                            (Nx  -2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+2;];
                else
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (Nx  -2)*Ny*Nz+(Ny-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+k+1;];
                end
            elseif j==Ny
                if k==1
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(2-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (Nx  -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(2-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+k+1;];
                elseif k==Nz
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(2-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+2;
                                            (Nx  -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(2-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+2;];
                else
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(2-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (Nx  -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(2-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+k+1;];
                end
            else
                if k==1
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (Nx  -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+k+1;];
                elseif k==Nz
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+2;
                                            (Nx  -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+2;];
                else
                    spJ4Au[index+1:index+10,1]=[(Nx-1-2)*Ny*Nz+(j-1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k; (Nx-1-2)*Nz*Ny+(j+1-1)*Nz+k; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k-1; (Nx-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                            (Nx  -2)*Ny*Nz+(j-1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k; (Nx  -2)*Nz*Ny+(j+1-1)*Nz+k; (Nx  -2)*Ny*Nz+(j-1)*Nz+k-1; (Nx  -2)*Ny*Nz+(j-1)*Nz+k+1;];
                end
            end
            spI4Au[index+1:index+10,1]=ones(10)*((Nx-2)*Ny*Nz+(j-1)*Nz+k);
            spI4bu[index4b+1,1]=(Nx-2)*Ny*Nz+(j-1)*Nz+k;
            index=index+10; index4b=index4b+1;
        end
    end
    for i=3:Nx-1 #rightest grids
        for j=1:Ny #k=Nz
            if j==1
                spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(Ny-1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz; (i-1-2)*Nz*Ny+(j+1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i  -2)*Ny*Nz+(Ny-1-1)*Nz+Nz; (i  -2)*Ny*Nz+(j-1)*Nz+Nz; (i  -2)*Nz*Ny+(j+1-1)*Nz+Nz; (i  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i  -2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i+1-2)*Ny*Nz+(Ny-1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz; (i+1-2)*Nz*Ny+(j+1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
            elseif j==Ny
                spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz; (i-1-2)*Nz*Ny+(1+1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i  -2)*Ny*Nz+(j-1-1)*Nz+Nz; (i  -2)*Ny*Nz+(j-1)*Nz+Nz; (i  -2)*Nz*Ny+(1+1-1)*Nz+Nz; (i  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i  -2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i+1-2)*Ny*Nz+(j-1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz; (i+1-2)*Nz*Ny+(1+1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
            else
                spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz; (i-1-2)*Nz*Ny+(j+1-1)*Nz+Nz; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i  -2)*Ny*Nz+(j-1-1)*Nz+Nz; (i  -2)*Ny*Nz+(j-1)*Nz+Nz; (i  -2)*Nz*Ny+(j+1-1)*Nz+Nz; (i  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i  -2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i+1-2)*Ny*Nz+(j-1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz; (i+1-2)*Nz*Ny+(j+1-1)*Nz+Nz; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
            end
            spI4Au[index+1:index+15]=((i-2)*Ny*Nz+(j-1)*Nz+Nz)*ones(15);
            spI4bu[index4b+1,1]=(i-2)*Ny*Nz+(j-1)*Nz+Nz;
            index=index+15; index4b=index4b+1;
        end
    end
    for i=3:Nx-1 #leftest grids
        for j=1:Ny #k=1
            if j==1
            spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(Ny-1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+1; (i-1-2)*Nz*Ny+(j+1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i  -2)*Ny*Nz+(Ny-1-1)*Nz+1; (i  -2)*Ny*Nz+(j-1)*Nz+1; (i  -2)*Nz*Ny+(j+1-1)*Nz+1; (i  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i  -2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i+1-2)*Ny*Nz+(Ny-1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+1; (i+1-2)*Nz*Ny+(j+1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
            elseif j==Ny
            spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+1; (i-1-2)*Nz*Ny+(1+1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i  -2)*Ny*Nz+(j-1-1)*Nz+1; (i  -2)*Ny*Nz+(j-1)*Nz+1; (i  -2)*Nz*Ny+(1+1-1)*Nz+1; (i  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i  -2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i+1-2)*Ny*Nz+(j-1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+1; (i+1-2)*Nz*Ny+(1+1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
            else
            spJ4Au[index+1:index+15]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+1; (i-1-2)*Nz*Ny+(j+1-1)*Nz+1; (i-1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i-1-2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i  -2)*Ny*Nz+(j-1-1)*Nz+1; (i  -2)*Ny*Nz+(j-1)*Nz+1; (i  -2)*Nz*Ny+(j+1-1)*Nz+1; (i  -2)*Ny*Nz+(j-1)*Nz+Nz-1; (i  -2)*Ny*Nz+(j-1)*Nz+1+1;
                                      (i+1-2)*Ny*Nz+(j-1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+1; (i+1-2)*Nz*Ny+(j+1-1)*Nz+1; (i+1-2)*Ny*Nz+(j-1)*Nz+Nz-1; (i+1-2)*Ny*Nz+(j-1)*Nz+1+1;];
            end
            spI4Au[index+1:index+15]=ones(15)*((i-2)*Ny*Nz+(j-1)*Nz+1);
            spI4bu[index4b+1,1]=(i-2)*Ny*Nz+(j-1)*Nz+1;
            index=index+15; index4b=index4b+1;
        end
    end
    for i=3:Nx-1 #upper grids j=1
        j=1;
        for k=2:Nz-1
            spI4Au[index+1:index+15,1]=((i-2)*Ny*Nz+(1-1)*Nz+k)*ones(15);
            spJ4Au[index+1:index+15,1]=[(i-1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k; (i-1-2)*Nz*Ny+(j+1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k-1; (i-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                    (i  -2)*Ny*Nz+(Ny-1-1)*Nz+k; (i  -2)*Ny*Nz+(j-1)*Nz+k; (i  -2)*Nz*Ny+(j+1-1)*Nz+k; (i  -2)*Ny*Nz+(j-1)*Nz+k-1; (i  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                    (i+1-2)*Ny*Nz+(Ny-1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k; (i+1-2)*Nz*Ny+(j+1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k-1; (i+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
            spI4bu[index4b+1,1]= (i-2)*Ny*Nz+(j-1)*Nz+k;
            index=index+15; index4b=index4b+1;
        end
    end
    for i=3:Nx-1 #upper grids j=Ny
        j=Ny;
        for k=2:Nz-1
            spI4Au[index+1:index+15,1]=ones(15)*((i-2)*Ny*Nz+(j-1)*Nz+k);
            spJ4Au[index+1:index+15,1]=[(i-1-2)*Ny*Nz+(j-1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k; (i-1-2)*Nz*Ny+(1+1-1)*Nz+k; (i-1-2)*Ny*Nz+(j-1)*Nz+k-1; (i-1-2)*Ny*Nz+(j-1)*Nz+k+1;
                                    (i  -2)*Ny*Nz+(j-1-1)*Nz+k; (i  -2)*Ny*Nz+(j-1)*Nz+k; (i  -2)*Nz*Ny+(1+1-1)*Nz+k; (i  -2)*Ny*Nz+(j-1)*Nz+k-1; (i  -2)*Ny*Nz+(j-1)*Nz+k+1;
                                    (i+1-2)*Ny*Nz+(j-1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k; (i+1-2)*Nz*Ny+(1+1-1)*Nz+k; (i+1-2)*Ny*Nz+(j-1)*Nz+k-1; (i+1-2)*Ny*Nz+(j-1)*Nz+k+1;];
            spI4bu[index4b+1,1]= (i-2)*Ny*Nz+(j-1)*Nz+k;
            index=index+15; index4b=index4b+1;
        end
    end
    DataG=zeros(Nx,Ny,Nz)
    for i=1:Nx
        for j=1:Ny
            for k=1:Nz
                DataG[i,j,k]=-exp(-(d-Datah[j,k])*(1-x[i]))#-exp(-((d-Datah[j,k])*(1-x[i])).^2)
            end
        end
    end
    #-2sin(d*x[i])*(sin(L_y*y[j])+1)*(sin(L_z*z[k])+1)-sin(d*x[i])*sin(L_y*y[j])*(sin(L_z*z[k])+1)-sin(d*x[i])*(sin(L_y*y[j])+1)*sin(L_z*z[k])
    ## assmble and solve for u
    index=0;index4b=0;
    for i=3:Nx-1
        for j=2:Ny-1
            for k=2:Nz-1
                spS4Au[index+1:index+15,1]=rand(15);
                spS4bu[index4b+1,1]=1#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
                index=index+15; index4b=index4b+1;
            end
        end
    end
    i=2;
    for j=1:Ny   #front grids,i=2
        for k=1:Nz
            spS4Au[index+1:index+10,1]=rand(10);
            spS4bu[index4b+1,1]=1#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
            index=index+10; index4b=index4b+1;
        end
    end
    i=Nx;
    for j=1:Ny    #back grids
        for k=1:Nz
            spS4Au[index+1:index+10,1]=rand(10)
            spS4bu[index4b+1,1]=1#-#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
            index=index+10; index4b=index4b+1;
        end
    end
    k=Nz;
    for i=3:Nx-1 #rightest grids
        for j=1:Ny #k=Nz
            spS4Au[index+1:index+15,1]=rand(15)
            spS4bu[index4b+1,1]=1
            index=index+15; index4b=index4b+1;
        end
    end
    for i=3:Nx-1 #leftest grids
        k=1;
        for j=1:Ny #k=1
                spS4Au[index+1:index+15,1]=rand(15);
            spS4bu[index4b+1,1]=1#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
            index=index+15; index4b=index4b+1;
        end
    end
    for i=3:Nx-1 #upper grids j=1
    j=1;
        for k=2:Nz-1
            spS4Au[index+1:index+15,1]=rand(15);
            spS4bu[index4b+1,1]=1#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
            index=index+15; index4b=index4b+1;
        end
    end
    for i=3:Nx-1 #upper grids j=Ny
        j=Ny;
        for k=2:Nz-1
            spS4Au[index+1:index+15,1]=rand(15);
            spS4bu[index4b+1,1]=1#-exp(-λ*(d-Datah[j,k])*(1-x[i]));
            index=index+15; index4b=index4b+1;
        end
    end
    Au =sparse(vec(spI4Au),vec(spJ4Au),vec(spS4Au),DOFu,DOFu);
    bu = sparse(vec(spI4bu),vec(spJ4bu),vec(spS4bu),DOFu,1);
    u=Au\Matrix(bu)
    u= [zeros(Ny*Nz,1); u];
    U=reshape(Matrix(u),Nx,Ny,Nz);
    return U
end

``

I mean, for large enough systems you will run out of memory even though you use a sparse storage format. How large is the matrix? And how sparse?

A = 1020100×1020100 SparseMatrixCSC{Float64,Int64} with 15199490 stored entries
b = 1020100×1 SparseMatrixCSC{Float64,Int64} with 1020100 stored entries:
But A is a block diagonal matrix not a diagonal matrix.
A and b can be stored but canot divide, the problem seems happen in LU part

ERROR: OutOfMemoryError()
Stacktrace:
 [1] umfpack_numeric!(::SuiteSparse.UMFPACK.UmfpackLU{Float64,Int64}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\SuiteSparse\src\umfpack.jl:258
 [2] #lu#1(::Bool, ::typeof(lu), ::SparseMatrixCSC{Float64,Int64}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\SuiteSparse\src\umfpack.jl:160
 [3] lu at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\SuiteSparse\src\umfpack.jl:154 [inlined]
 [4] \(::SparseMatrixCSC{Float64,Int64}, ::Array{Float64,2}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\SparseArrays\src\linalg.jl:1374
 [5] pdesolver(::Int64, ::Float64, ::Float64, ::Float64, ::Array{Float64,2}, ::Int64) at C:\Users\lyu\Downloads\test.jl:257

b can be stored in a Vector. Dont need sparse vector for this.

You should realize that there is an intrinsic limitation of sparse-direct methods like this for matrices that come from large 3d meshes/grids — they scale poorly for 3d problems, and tend to run out of memory before the problems become large. You could try a different sparse-direct solver like MUMPS (https://github.com/JuliaSmoothOptimizers/MUMPS.jl or https://github.com/JuliaInv/MUMPSjInv.jl), which might suffice in the short term, but in the long term you will hit the same scaling issues.

Here is a nice slide on the scaling (by Per Persson) that I use in one of my classes:

Note that for an n \times n \times n grid, N=n^3 so the storage scales as O(n^4). For more information, see e.g. Tim Davis’ book.

In consequence, for large 3d PDEs, the sparse-matrix problems must ultimately be solved by iterative methods (e.g. GMRES, CG, etc.) rather than sparse-direct methods. There are a few packages for iterative solvers in Julia, e.g. IterativeSolvers.jl, but to use iterative solvers effectively you have to know a bit more, and you often want to find a problem-specific preconditioner.

9 Likes