Mercurial > repos > public > sbplib_julia

comparison diffOp.jl @ 87:38733e84ef1a patch_based_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up bounds checking in stencil. change to using uview in DiffOp
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 25 Jan 2019 13:40:15 +0100
parents 7f72e7e14659
children 170e5447bc19
comparison
equal deleted inserted replaced
86:34fd86e9d0b9 87:38733e84ef1a
46 u .= L.a * u 46 u .= L.a * u
47 return nothing 47 return nothing
48 end 48 end
49 49
50 # u = L*v 50 # u = L*v
51 using UnsafeArrays
51 @inline function apply!(L::Laplace{2}, u::AbstractVector, v::AbstractVector) 52 @inline function apply!(L::Laplace{2}, u::AbstractVector, v::AbstractVector)
52 u .= 0*u 53 u .= 0*u # Fix this?
53 h = Grid.spacings(L.grid) 54 h = Grid.spacings(L.grid)
54 55
55 li = LinearIndices(L.grid.numberOfPointsPerDim) 56 li = LinearIndices(L.grid.numberOfPointsPerDim)
56 n_x, n_y = L.grid.numberOfPointsPerDim 57 n_x, n_y = L.grid.numberOfPointsPerDim
57 58
58 59
59 # For each x 60 # For each x
60 temp = zeros(eltype(u), n_y) 61 temp = zeros(eltype(u), n_y)
61 for i ∈ 1:n_x 62 for i ∈ 1:n_x
62 63 @inbounds indices = uview(li,i,:)
63 @inbounds v_i = view(v, li[i,:]) 64 @inbounds apply!(L.op, temp, uview(v, indices), h[2])
64 @inbounds apply!(L.op, temp, v_i, h[2]) 65 for i ∈ eachindex(indices)
65 66 @inbounds u[indices[i]] += temp[i]
66 @inbounds u[li[i,:]] += temp 67 end
67 end 68 end
68 69
69 # For each y 70 # For each y
70 temp = zeros(eltype(u), n_x) 71 temp = zeros(eltype(u), n_x)
71 for i ∈ 1:n_y 72 for i ∈ 1:n_y
72 @inbounds v_i = view(v, li[:,i]) 73 @inbounds indices = uview(li,:,i)
73 @inbounds apply!(L.op, temp, v_i, h[1]) 74 @inbounds apply!(L.op, temp, uview(v, indices), h[1])
74 75 for i ∈ eachindex(indices)
75 @inbounds u[li[:,i]] += temp 76 @inbounds u[indices[i]] += temp[i]
77 end
76 end 78 end
77 79
78 u .= L.a*u 80 u .= L.a*u
79 81
80 return nothing 82 return nothing