Demonstrate bug when using MKL.jl and MKLSparse.jl

calcpml3d.jl

"""
    calcpml3d(NGRID, NPML) -> (sx, sy, sz)

Calculate PML Parameters

# Input Arguments

- `NGRID`: Iterable containing `(Nx, Ny, Nz)`,  number of cells on grid
- `NPML`:  Iterable containing (NXLO, NXHI, NYLO, NYHI, NZLO, NZHI), the size of PML

# Return Values 

- `sx`, `sy`, `sz`: Complex matrices of size  `(Nx, Ny, Nz)`, resp., containing the 
  complex coordinate stretching parameters for each dimension of the grid.
"""
function calcpml3d(NGRID,NPML)

    #########################################################
    ## HANDLE INPUT ARGUMENTS
    #########################################################

    # DEFINE PML PARAMETERS
    amax = 4
    cmax = 1
    p = 3

    # EXTRACT GRID SIZE
    Nx, Ny, Nz = (n for n in NGRID)

    # EXTRACT PML SIZE
    NXLO, NXHI, NYLO, NYHI, NZLO, NZHI = (n for n in NPML)

    #########################################################
    ## CALCULATE PML PARAMETERS
    #########################################################

    # INITIALIZE PML PARAMETERS TO PROBLEM SPACE
    sx, sy, sz = (ones(ComplexF64, Nx, Ny, Nz) for _ in 1:3)

    # CALCULATE XLO PML
    for nx in 1:NXLO
        ax = 1 + (amax - 1) * (nx/NXLO)^p
        cx = cmax * sinpi(nx/(2*NXLO))^2;
        sx[NXLO - nx + 1, :, :] .= ax * (1 - im * 60cx)
    end

    # CALCULATE XHI PML
    for nx in 1:NXHI
        ax = 1 + (amax - 1) * (nx/NXHI)^p
        cx = cmax * sinpi(nx/(2*NXHI))^2
        sx[Nx - NXHI + nx, :, :] .= ax * (1 - im * 60cx)
    end

    # CALCULATE YLO PML
    for ny in 1:NYLO
        ay = 1 + (amax - 1) * (ny/NYLO)^p
        cy = cmax * sinpi(ny/(2*NYLO))^2
        sy[:, NYLO - ny + 1, :] .= ay * (1 - im * 60cy)
    end

    # CALCULATE YHI PML
    for ny in 1:NYHI
        ay = 1 + (amax - 1) * (ny/NYHI)^p
        cy = cmax * sinpi(ny/(2*NYHI))^2
        sy[:, Ny - NYHI + ny, :] .= ay * (1 - im * 60cy)
    end
    # CALCULATE ZLO PML
    for nz in 1:NZLO
        az = 1 + (amax - 1) * (nz/NZLO)^p
        cz = cmax * sinpi(nz/(2*NZLO))^2
        sz[:, :, NZLO - nz + 1] .= az * (1 - im * 60cz)
    end

    # CALCULATE ZHI PML
    for nz in 1:NZHI
        az = 1 + (amax - 1) * (nz/NZHI)^p
        cz = cmax * sinpi(nz/(2*NZHI))^2
        sz[:, :, Nz - NZHI + nz] .= az * (1 - im * 60cz)
    end

    return (sx, sy, sz)
end

fdfd3d.jl

"""
    fdfd3d(dev::NamedTuple, src::NamedTuple) -> dat::NamedTuple

 Three-dimensional FDFD for periodic structures

# INPUT ARGUMENTS
- `dev`: Device parameters, with the following fields:
    - `ER2` or
      `ER2xx`   `ER2xy`   `ER2xz`
      `ER2yx`   `ER2yy`   `ER2yz`     relative permittivity on 2x grid
      `ER2zx`   `ER2zy`   `ER2zz`
    - `UR2`: or
       `UR2xx`   `UR2xy`   `UR2xz`
       `UR2yx`   `UR2yy`   `UR2yz`     relative permeability on 2x grid
       `UR2zx`   `UR2zy`   `UR2zz`
    - `RES`: `[dx, dy]` grid resolution on Yee grid
    - `NPML`: `[NZLO, NZHI]` size of SCPML on Yee grid
- `src`: Source parameters with the following fields:
    - `λ₀`: free-space wavelength
    - `θ`: polar angle of incidence
    - `ϕ`: azimuthal angle of incidence
    - `pte`: complex amplitude of the TE polarization
    - `ptm`: complex amplitude of the TM polarization

# RETURN VALUES
- `dat`: Output data with the following fields:
    - `fx`, `fy`, `fz`:  simulated field
    - `m`, `n`: diffraction order indices
    - `RDE`:  diffraction efficiencies of reflected waves
    - `REF`: overall reflectance
    - `TDE`: diffraction efficiencies of transmitted waves
    - `TRN`: overall transmittance
    - `s11`: complex ref. coefficients of diffraction orders
    - `s21`: complex tran. coefficients of diffraction orders
    - `CON`: Power conservation parameter, equal to `REF + TRN`
"""
function fdfd3d(dev,src)

    # DEFINE SOLVER PARAMETERS
    tol   = 1e-6
    maxit = 15000

    # UTILITY FUNCTIONS
    diagonalize(x) = Diagonal(vec(x))

    #########################################################
    ## HANDLE INPUT ARGUMENTS
    #########################################################

    # GET SIZE OF GRID
    if haskey(dev, :ER2)
        Nx2, Ny2, Nz2 = size(dev.ER2)
    else
        Nx2, Ny2, Nz2 = size(dev.ER2xx)
    end

    # CALCULATE GRID
    Nx, Ny, Nz = (Nx2, Ny2, Nz2) .÷ 2
    dx, dy, dz = dev.RES
    Sx, Sy, Sz = (Nx, Ny, Nz) .* (dx, dy, dz)

    # GRID AXES
    xa2 = [1:Nx2;] * dx/2;  xa2 .-= mean(xa2)
    ya2 = [1:Ny2;] * dy/2;  ya2 .-= mean(ya2)
    za2 = [0.5:Nz2-0.5;] * dz/2
    #[Y2,X2,Z2] = meshgrid(ya2,xa2,za2)

    # WAVE NUMBER
    k0 = 2π / src.λ₀

    # SPECIAL MATRICES
    M  = Nx * Ny * Nz
    ZZ = spzeros(M, M)
    #I0  = speye(M,M);

    # PERMITTIVITY TENSOR ELEMENTS
    if haskey(dev, :ER2)
        ERxx = diagonalize(dev.ER2[2:2:Nx2, 1:2:Ny2,1:2:Nz2])
        ERxy = ZZ
        ERxz = ZZ
        ERyx = ZZ
        ERyy = diagonalize(dev.ER2[1:2:Nx2, 2:2:Ny2,1:2:Nz2])
        ERyz = ZZ
        ERzx = ZZ
        ERzy = ZZ
        ERzz = diagonalize(dev.ER2[1:2:Nx2, 1:2:Ny2,2:2:Nz2])
    else
        ERxx = diagonalize(dev.ER2xx[2:2:Nx2, 1:2:Ny2,1:2:Nz2])
        ERxy = diagonalize(dev.ER2xy[1:2:Nx2, 2:2:Ny2,1:2:Nz2])
        ERxz = diagonalize(dev.ER2xz[1:2:Nx2, 1:2:Ny2,2:2:Nz2])
        ERyx = diagonalize(dev.ER2yx[2:2:Nx2, 1:2:Ny2,1:2:Nz2])
        ERyy = diagonalize(dev.ER2yy[1:2:Nx2, 2:2:Ny2,1:2:Nz2])
        ERyz = diagonalize(dev.ER2yz[1:2:Nx2, 1:2:Ny2,2:2:Nz2])
        ERzx = diagonalize(dev.ER2zx[2:2:Nx2, 1:2:Ny2,1:2:Nz2])
        ERzy = diagonalize(dev.ER2zy[1:2:Nx2, 2:2:Ny2,1:2:Nz2])
        ERzz = diagonalize(dev.ER2zz[1:2:Nx2, 1:2:Ny2,2:2:Nz2])
    end

    # PERMEABILITY TENSOR ELEMENTS
    if haskey(dev, :UR2)
        URxx = diagonalize(dev.UR2[2:2:Nx2, 1:2:Ny2,1:2:Nz2])
        URxy = ZZ
        URxz = ZZ
        URyx = ZZ
        URyy = diagonalize(dev.UR2[1:2:Nx2, 2:2:Ny2,1:2:Nz2])
        URyz = ZZ
        URzx = ZZ
        URzy = ZZ
        URzz = diagonalize(dev.UR2[1:2:Nx2, 1:2:Ny2,2:2:Nz2])
    else
        URxx = diagonalize(dev.UR2xx[2:2:Nx2, 1:2:Ny2,1:2:Nz2])
        URxy = diagonalize(dev.UR2xy[1:2:Nx2, 2:2:Ny2,1:2:Nz2])
        URxz = diagonalize(dev.UR2xz[1:2:Nx2, 1:2:Ny2,2:2:Nz2])
        URyx = diagonalize(dev.UR2yx[2:2:Nx2, 1:2:Ny2,1:2:Nz2])
        URyy = diagonalize(dev.UR2yy[1:2:Nx2, 2:2:Ny2,1:2:Nz2])
        URyz = diagonalize(dev.UR2yz[1:2:Nx2, 1:2:Ny2,2:2:Nz2])
        URzx = diagonalize(dev.UR2zx[2:2:Nx2, 1:2:Ny2,1:2:Nz2])
        URzy = diagonalize(dev.UR2zy[1:2:Nx2, 2:2:Ny2,1:2:Nz2])
        URzz = diagonalize(dev.UR2zz[1:2:Nx2, 1:2:Ny2,2:2:Nz2])
    end

    #########################################################
    ## PERFORM FDFD ANALYS
    #########################################################

    # EXTRACT MATERIAL PROPERTIES IN EXTERNAL REGIONS
    ur_ref = URxx[1, 1]
    ur_trn = URxx[M, M]
    er_ref = ERxx[1, 1]
    er_trn = ERxx[M, M]
    n_ref  = sqrt(ur_ref * er_ref)
    n_trn  = sqrt(ur_trn * er_trn)

    # CALCULATE PML PARAMETERS
    sx2, sy2, sz2 = calcpml3d([Nx2 Ny2 Nz2], 2 .* [0, 0, 0, 0, dev.NPML...])

    sx_ey = diagonalize(1 ./ sx2[1:2:Nx2,2:2:Ny2,1:2:Nz2])
    sx_ez = diagonalize(1 ./ sx2[1:2:Nx2,1:2:Ny2,2:2:Nz2])
    sy_ex = diagonalize(1 ./ sy2[2:2:Nx2,1:2:Ny2,1:2:Nz2])
    sy_ez = diagonalize(1 ./ sy2[1:2:Nx2, 1:2:Ny2,2:2:Nz2])
    sz_ex = diagonalize(1 ./ sz2[2:2:Nx2, 1:2:Ny2,1:2:Nz2])
    sz_ey = diagonalize(1 ./ sz2[1:2:Nx2, 2:2:Ny2,1:2:Nz2])

    sx_hy = diagonalize(1 ./ sx2[2:2:Nx2, 1:2:Ny2,2:2:Nz2])
    sx_hz = diagonalize(1 ./ sx2[2:2:Nx2, 2:2:Ny2,1:2:Nz2])
    sy_hx = diagonalize(1 ./ sy2[1:2:Nx2, 2:2:Ny2,2:2:Nz2])
    sy_hz = diagonalize(1 ./ sy2[2:2:Nx2, 2:2:Ny2,1:2:Nz2])
    sz_hx = diagonalize(1 ./ sz2[1:2:Nx2, 2:2:Ny2,2:2:Nz2])
    sz_hy = diagonalize(1 ./ sz2[2:2:Nx2, 1:2:Ny2,2:2:Nz2])

    # CALCULATE WAVE VECTORS
    kinc   = k0 * n_ref .* [sin(src.θ) * cos(src.ϕ),
                            sin(src.θ) * sin(src.ϕ),
                            cos(src.θ)]
    m = [-(Nx ÷ 2):((Nx-1) ÷2);]
    n = [-(Ny ÷ 2):((Ny-1) ÷ 2);]
    kx = kinc[1] .- m * 2π / Sx
    ky = kinc[2] .- n * 2π / Sy
    kz_ref = sqrt.(complex(((k0*n_ref)^2 .- kx.^2 .- ky.^2)))
    kz_trn = sqrt.(complex(((k0*n_trn)^2 .- kx.^2 .- ky.^2)))

    # BUILD DERIVATIVE MATRICES
    NS  = [Nx, Ny, Nz]
    RES = [dx, dy, dz]
    BC  = [1, 1, 0]
    DEX, DEY, DEZ, DHX, DHY, DHZ = yeeder3d(NS, k0*RES, BC, kinc/k0)

    # CALCULATE INTERPOLATION MATRICES (assumes src.θ == 0)
    iszero(src.θ) || error("Nonzero θ=$θ not supported for RX, RY, and RZ calculation")
    RX = (0.5 * k0 * dx) .* abs.(DEX)
    RY = (0.5 * k0 * dy) .* abs.(DEY)
    RZ = (0.5 * k0 * dz) .* abs.(DEZ)

    # FORM MATERIALS TENSORS
    ER = [ERxx         RX*RY'*ERxy  RX*RZ'*ERxz
          RY*RX'*ERyx  ERyy         RY*RZ'*ERyz 
          RZ*RX'*ERzx  RZ*RY'*ERzy  ERzz]
    UR = [URxx         RX*RY'*URxy  RX*RZ'*URxz
          RY*RX'*URyx  URyy         RY*RZ'*URyz 
          RZ*RX'*URzx  RZ*RY'*URzy  URzz]

    # BUILD THE WAVE MATRIX A
    CE = [ZZ         -sz_hx*DEZ  sy_hx*DEY 
          sz_hy*DEZ   ZZ        -sx_hy*DEX
         -sy_hz*DEY   sx_hz*DEX  ZZ]
    CH = [ZZ         -sz_ex*DHZ  sy_ex*DHY 
          sz_ey*DHZ   ZZ        -sx_ey*DHX
         -sy_ez*DHY   sx_ez*DHX  ZZ]
    A = CH * (UR \ CE) - ER

    # CALCULATE POLARIZATION VECTOR P
    az = [0.0, 0.0, 1.0]

    if abs(src.θ) < 1e-6
        ate = [0.0, 1.0, 0.0]
    else
        ate = az × kinc
        ate /= norm(ate)
    end
    atm = kinc × ate
    atm /= norm(atm)
    P = src.pte * ate + src.ptm * atm

    # CALCULATE SOURCE FIELD fsrc
    fphase = [cis(-(kinc[1]*x + kinc[2]*y + kinc[3]*z)) for y in ya2, x in xa2, z in za2]
    fx = P[1] * fphase[2:2:Nx2, 1:2:Ny2, 1:2:Nz2]
    fy = P[2] * fphase[1:2:Nx2, 2:2:Ny2, 1:2:Nz2]
    fz = P[3] * fphase[1:2:Nx2, 1:2:Ny2, 2:2:Nz2]
    fsrc = [vec(fx); vec(fy); vec(fz)]

    # CALCULATE SCATERED FIELD MASKING MATRIX
    nz = dev.NPML[1] + 3
    Q = zeros(Nx, Ny, Nz)
    Q[:, :, 1:nz] .= 1
    Q = Diagonal(repeat(vec(Q), 3))
    #Q = diagonalize(Q)
    #Q = [Q ZZ ZZ; ZZ Q ZZ; ZZ ZZ Q]


    # CALCULATE SOURCE VECTOR B
    QAAQ = Q*A - A*Q
    b = QAAQ * fsrc

    @show typeof.((QAAQ, fsrc, b))
    @show size.((QAAQ, fsrc, b))
    @show nnz(QAAQ)
    @show maximum(abs ∘ imag, QAAQ)
    @show maximum(abs, QAAQ - real(QAAQ))
    @show maximum(abs, b - real(QAAQ) * fsrc)


    #=
    # COMPUTE JACOBI PRECONDITIONER
    Pl = Diagonal(inv.(Array(diag(A))))

    # SOLVE FOR FIELD
    Agpu = CuSparseMatrixCSC(A)
    bgpu = CuArray(b)
    Plgpu = CuArray(Pl)

    f, hist = idrs(A, b; Pl=Pl, log=true, verbose=true, maxiter=2000, abstol=1e-8, reltol=1e-6) 
    @time fgpu, histgpu = idrs(Agpu, bgpu; Pl=Plgpu, log=true, maxiter=2000, abstol=1e-8, reltol=1e-6) 

    #f = bicg(A,b,tol,maxit,[],[],f);
    #f = gather(f);

    # EXTRACT FIELD COMPONENTS
    fx = reshape(f[1:M], Nx, Ny, Nz)
    fy = reshape(f[M+1:2*M], Nx, Ny, Nz)
    fz = reshape(f[2*M+1:3*M], Nx, Ny, Nz)

    #########################################################
    ## ANALYZE REFLECTION AND TRANSMISSION
    #########################################################

    # EXTRACT FIELD FROM RECORD PLANES
    nz_ref = dev.NPML[1] + 2
    fx_ref = fx[:, :, nz_ref]
    fy_ref = fy[:, :, nz_ref]
    fz_ref = fz(:, :, nz_ref-1:nz_ref)

    nz_trn = Nz - dev.NPML[1]
    fx_trn = fx[:, :, nz_trn]
    fy_trn = fy[:, :, nz_trn]
    fz_trn = fz[:, :, nz_trn-1:nz_trn]

    # INTERPOLATE FIELDS TO ORIGIN OF YEE CELLS
    Ex_ref = copy(fx_ref)
    Ex_ref[1, :] = (fx_ref[Nx, :] + fx_ref[1, :]) / 2
    Ex_ref[2:Nx, :] = (fx_ref[1:Nx-1, :] + fx_ref[2:Nx, :]) / 2

    Ey_ref = copy(fy_ref)
    Ey_ref[:,1] = (fy_ref[:, Ny] + fy_ref[:, 1]) / 2
    Ey_ref[:, 2:Ny] = (fy_ref[:, 1:Ny-1] + fy_ref[:, 2:Ny]) / 2

    Ez_ref = (fz_ref[:,:,1] + fz_ref[:, :, 2]) / 2

    Ex_trn = copy(fx_trn)
    Ex_trn[1,:] = (fx_trn[Nx, :] + fx_trn[1, :]) / 2
    Ex_trn[2:Nx, :] = (fx_trn[1:Nx-1, :] + fx_trn[2:Nx, :]) / 2

    Ey_trn = copy(fy_trn)
    Ey_trn[:, 1]    = (fy_trn[:, Ny] + fy_trn[:, 1]) / 2
    Ey_trn[:, 2:Ny] = (fy_trn[:, 1:Ny-1] + fy_trn[:, 2:Ny]) / 2

    Ez_trn = (fz_trn[:, :, 1] + fz_trn[:, :, 2]) / 2

    # REMOVE PHASE TILT
    Ex_ref ./= fphase[2:2:Nx2, 1:2:Ny2, 1]
    Ey_ref ./= fphase[1:2:Nx2, 2:2:Ny2, 1]
    Ez_ref ./= fphase[1:2:Nx2, 1:2:Ny2, 2]

    Ex_trn ./= fphase[2:2:Nx2, 1:2:Ny2, 1]
    Ey_trn ./= fphase[1:2:Nx2, 2:2:Ny2, 1]
    Ez_trn ./= fphase[1:2:Nx2, 1:2:Ny2, 2]

    # CALCULATE AMPLITUDES OF DIFFRACTION ORDERS
    NxNy = Nx * Ny
    ax_ref = fftshift(fft(Ex_ref)) / NxNy
    ay_ref = fftshift(fft(Ey_ref)) / NxNy
    az_ref = fftshift(fft(Ez_ref)) / NxNy
    a_ref  = abs2.(ax_ref) + abs2.(ay_ref) + abs2.(az_ref)
    s11 = sqrt.(ax_ref.^2 + ay_ref.^2 + az_ref.^2)

    ax_trn = fftshift(fft(Ex_trn)) / NxNy
    ay_trn = fftshift(fft(Ey_trn)) / NxNy
    az_trn = fftshift(fft(Ez_trn)) / NxNy
    a_trn  = abs2.(ax_trn) + abs2.(ay_trn) + abs2.(az_trn)
    s21 = sqrt(ax_trn.^2 + ay_trn.^2 + az_trn.^2)

    # CALCULATE DIFFRACTION EFFICIENCIES
    RDE = a_ref .* real(kz_ref / kinc[3])
    TDE = a_trn .* real(ur_ref / ur_trn * kz_trn / kinc[3])

    # Calculate Overall Response
    REF = sum(RDE)
    TRN = sum(TDE)
    CON = REF + TRN;


    return (; fx, fy, fz, m, n, RDE, REF, TDE, TRN, s11, s21, CON)
=#
end

reproducer.jl

using LinearAlgebra
using SparseArrays
#using MKL, MKLSparse
using Statistics: mean

include("yeeder3d.jl")
include("calcpml3d.jl")
include("fdfd3d.jl")


# DEFINE UNITS
const meters = 1
const centimeters = 1e-2meters
const seconds = 1
const hertz = 1 / seconds
const gigahertz = 1e9 * hertz
const degrees = π / 180

# CONSTANTS
const c0 = 299792458 * meters/seconds

#########################################################
## DASHBOARD
#########################################################

# DEFINE SOURCE PARAMETERS
f1 = 4.5gigahertz
f2 = 5.5gigahertz
NFREQ = 500
FREQ = LinRange(f1, f2, NFREQ)
θ = ϕ = 0degrees
pte = 1 / sqrt(2)
ptm = im / sqrt(2)

# DEFINE GMRF PARAMETERS
er1 = 1.0
er2 = 2.0
erL = 3.0
erH = 5.0
a = 3.7centimeters
t = 1.5centimeters
frac = 0.5

# DEFINE GRID PARAMETERS
NRES = 20
NPML = [10, 10]
ermax = max(er1, er2, erL, erH)
nmax = sqrt(ermax)
lam_max  = c0 / minimum(FREQ)
SPACER = lam_max * [1, 1]

#########################################################
## CALCULATE OPTIMIZED GRID
#########################################################

# CALCULATE PRELIMINARY RESOLUTION
lam_min = c0 / maximum(FREQ)
dx = dy = dz = lam_min/nmax/NRES

# SNAP GRID TO CRITICAL DIMENSIONS
snap(w, Δ) = begin n = ceil(Int, w/Δ); δ = w/n; return (n, δ) end
nx, dx = snap(a, dx)
ny, dy = snap(a, dy)
nz, dz = snap(t, dz)

# CALCULATE GRID SIZE
Nx = ceil(Int, a / dx)
Sx = Nx * dx
Ny = ceil(Int, a / dy)
Sy = Ny * dy
Sz = SPACER[1] + t + SPACER[2]
Nz = NPML[1] + ceil(Int, Sz / dz) + NPML[2]
Sz = Nz * dz

# 2x GRID
Nx2, dx2 = (2Nx, dx/2)
Ny2, dy2 = (2Ny, dy/2)
Nz2, dz2 = (2Nz, dz/2)

# GRID AXES
xa = [1:Nx;] * dx; xa .-= mean(xa)
ya = [1:Ny;] * dy; ya .-= mean(ya)
za = [0.5:Nz-0.5;] * dz

xa2 = [1:Nx2;] * dx2; xa2 .-= mean(xa2)
ya2 = [1:Ny2;] * dy2; ya2 .-= mean(ya2)
za2 = [0.5:Nz2-0.5;] * dz2

# MESHGRIDS
#[Y,X,Z]    = meshgrid(ya,xa,za);
#[Y2,X2,Z2] = meshgrid(ya2,xa2,za2);

#########################################################
## BUILD MATERIALS ARRAYS
#########################################################

# INITIALIZE TO FREE SPACE
ER2 = ones(Nx2, Ny2, Nz2)
UR2 = ones(Nx2, Ny2, Nz2)

# BUILD GRATING LAYER
r = a * sqrt(frac / π)
r² = r * r
ER2 = [x^2 + y^2 ≤ r² ? erH : erL for y in ya2, x in xa2, z in za2]

# ADD SUPERSTRATE AND SUBSTRATE
nz1 = 2 * NPML[1] + round(Int, SPACER[1] / dz2) + 1
nz2 = nz1 + round(Int, t / dz2) - 1
ER2[:, :, 1:nz1-1]  .= er1
ER2[:, :, nz2+1:Nz2] .= er2;

#########################################################
## PERFORM FREQUENCY SWEEP
#########################################################

# dev and src
RES = [dx, dy, dz]
dev = (; NPML,  ER2, UR2, RES)

f0 = f1
λ₀ = c0 / f0
src = (; θ, ϕ, pte, ptm, λ₀)
    
# Call FDFD3D()
dat = fdfd3d(dev, src)
        

yeeder3d.jl

using LinearAlgebra, SparseArrays

"""
    yeeder3d(NS,RES,BC,kinc) -> (DEX,DEY,DEZ,DHX,DHY,DHZ)

Compute derivative matrices on a 3D Yee grid

# Input Arguments

- `NS`: `(Nx Ny Nz)` Grid size.
- `RES`: `(dx, dy, dz)` Grid resolution. For normalized grids contains `(k0*dx, k0*dy, k0*dz)`.
- `BC`:  `(xbc, ybc, zbc)` Boundary conditions, where 0(1) implies Dirichlet(periodic) boundary conditions.
- `kinc`: `[kx ky kz]` Incident wave vector (needed only for PBCs). For normalized grids contains `kinc/k0`.

# Return Values

- `DEX`: Derivative matrix wrt x for electric fields
- `DEY`: Derivative matrix wrt y for electric fields
- `DEZ`: Derivative matrix wrt z for electric fields
- `DHX`: Derivative matrix wrt x for magnetic fields
- `DHY`: Derivative matrix wrt y for magnetic fields
- `DHZ`: Derivative matrix wrt z for magnetic fields
"""
function yeeder3d(NS, RES, BC, kinc=zeros(3))

    #########################################################
    ## HANDLE INPUT ARGUMENTS
    #########################################################

    # EXTRACT GRID PARAMETERS
    Nx, Ny, Nz = NS
    dx, dy, dz = RES

    # DETERMINE MATRIX SIZE
    M = Nx * Ny * Nz

    # Determine if real or complex matrices are needed
    if (Nx == 1 && kinc[1] ≠ 0) || (Ny == 1 && kinc[2] ≠ 0) || (Nz == 1 && kinc[3] ≠ 0) || any(≠(0), BC)
        T = ComplexF64
    else
        T = Float64
    end

    #########################################################
    ## BUILD DEX
    #########################################################

    
    if Nx==1

        # HANDLE IF SIZE IS 1 CELL
        DEX = sparse(-im * kinc[1] * I, M, M)
    
    else

        # HANDLE ALL OTHER CASES

        d0 = fill(T(-1/dx), M) # Center Diagonal
        d1 = fill(T(1/dx), M) # Upper Diagonal
        d1[Nx+1:Nx:M] .= zero(T)
        d1 = d1[2:end]
        # Build Derivative Matrix with Dirichlet BCs
        DEX = spdiagm(M, M, 0 => d0, 1 => d1)
    
        # Incorporate Periodic Boundary Conditions
        if BC[1] == 1
            d1 = zeros(T, M - (Nx-1))
            d1[1:Nx:M] .= cis(-kinc[1] * (Nx*dx)) / dx
            DEX += spdiagm(M, M, 1-Nx => d1)
        end
    
    end

    #########################################################
    ## BUILD DEY
    #########################################################

    # HANDLE IF SIZE IS 1 CELL
    if Ny==1

        DEY = sparse(-im * kinc[2] * I, M, M)
    
    else
    
    # HANDLE ALL OTHER CASES
    
    d0 = fill(T(-1/dy), M) # Center Diagonal
    
    # Upper Diagonal
    d1 = [fill(T(1/dy), (Ny-1)*Nx); zeros(T, Nx)]
    d1 = repeat(d1, Nz-1)
    d1 = [zeros(T, Nx); d1; fill(T(1/dy), (Ny-1)*Nx)]
    d1 = d1[Nx+1:end]
    
    # Build Derivative Matrix with Dirichlet BCs
    DEY = spdiagm(M, M, 0 => d0, Nx => d1)
    
    # Incorporate Periodic Boundary Conditions
    if BC[2] == 1
        ph = cis(-kinc[2] * (Ny*dy)) / dy
        d1 = [fill(ph, Nx); zeros(T, (Ny-1)*Nx)]
        d1 = repeat(d1, Nz-1)
        d1 = [d1; fill(ph, (Ny-1)*Nx); zeros(T, Nx)]
        d1 = d1[1:M-Nx*(Ny-1)]
        DEY += spdiagm(M, M, -Nx*(Ny-1) => d1)
    end
    
end

    #########################################################
    ## BUILD DEZ
    #########################################################

    # HANDLE IF SIZE IS 1 CELL
    if Nz==1
        DEZ = sparse(-im * kinc[3] * I, M, M)
    else
        # HANDLE ALL OTHER CASES

        d0 = fill(T(1/dz), M)     # Center Diagonal

        # Build Derivative Matrix with Dirichlet BCs
        DEZ = spdiagm(M, M, 0 => -d0, Nx*Ny => d0[Nx*Ny+1:end])
    
        # Incorporate Periodic Boundary Conditions
        if BC[3] == 1
           d0  = fill(cis(-kinc[3] * (Nz*dz)) / dz, M - Nx*Ny*(Nz-1))
           DEZ += spdiagm(M, M, -Nx*Ny*(Nz-1) => d0)
        end
    
    end

    #########################################################
    ## BUILD DHX, DHY AND DHZ
    #########################################################

    DHX = sparse(-DEX')
    DHY = sparse(-DEY')
    DHZ = sparse(-DEZ')
    
    return (DEX, DEY, DEZ, DHX, DHY, DHZ)
end # function