'MultiVector' not updating in a loop

General Usage
1

Hello, I am using a for loop to insert values into a DenseMultiVector, a type that was created in the project I am continuing work on. I start by initializing an undef DenseMultiVector of the correct size, and then loop through by assigning each index of the DenseMultiVector with a new value. The new value, a sum, is being calculated correctly, but only the first value in each vector is being updated correctly. Any ideas of what the issue is?

Here is the code:

    comm = SerialComm{Int, Int, Int}()
    map = BlockMap(numVects, length, comm)
    data = Array{Float32}(undef, length, numVects)
    result = DenseMultiVector{}(map, data)

    A = scale!(A, a)
    B = scale!(B, b)

    for i = 1:numVectors(A)
        for j = 1:localLength(A)
            sum = A[i,j] + B[i,j]
            result[i,j] = sum
        end
    end
2

A couple of comments:

  1. scale! hasn’t been around for a long, long time (Julia 0.6 maybe?). Please update to a recent version of Julia.
  2. You are using several variables whose names are the same as built-in Julia functions, such as sum, map.
  3. Please provide a runnable minimal working example (MWE), item 4 in this post.
1 Like
3

It sounds like you might be trying to resurrect an old pre-1.0 project. To help you bring your code up to date with current Julia, you can install Julia 0.7 from here and run your code under that version. It will provide messages telling you what’s obsolete and what to replace the obsolete parts with.

4

The scale! function I am using here is not a built-in Julia function, it was written by the previous author of the project I am working on, specifically for type MultiVector.

using JuliaPetra
function update(a, A::MultiVector, b, B::MultiVector)
    if numVectors(A) == numVectors(B) && localLength(A) == localLength(B)
        myComm = SerialComm{Int, Int, Int}()
        myMap = BlockMap(numVectors(A), localLength(A), myComm)
        data = Array{Float32}(undef, localLength(A), numVectors(A))
        result = DenseMultiVector{}(myMap, data)

        A = scale!(A, a)
        B = scale!(B, b)

        for i = 1:numVectors(A)
             for j = 1:localLength(A)
                 result[i,j] = A[i,j] + B[i,j]
            end
         end

        else
            println("Error: MultiVectors must be the same size")
        end
        return result
    end

Is this more helpful?

*I would also like to note that when I printed out what A[i,j] + B[i,j] is in each iteration, the sum is correct which tells me the problem is in putting some of those values into the result MultiVector.

5

It would help a lot to see what the definition of DenseMultiVector is. Better yet would be a runnable MWE, to avoid misconceptions based on insufficient knowledge of your code (like my scale! comment).

Did you also try printing out the value of result[i,j] in the loop, say with a

@show result[i,j]

in the line right after the assignment?

6

Thank you for your help with this!
This is part of the project JuliaPetra (as I mentioned “using JuliaPetra”) so I’m not sure the best way to show you all the code, but here is the DenseMultiVector constructors:

function DenseMultiVector{Data}(map::BlockMap{GID, PID, LID}, numVecs::Integer, zeroOut=true) where {Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer}
     localLength = numMyElements(map)
     if zeroOut
         data = zeros(Data, (localLength, numVecs))
     else
         data = Array{Data, 2}(undef, localLength, numVecs)
     end
     DenseMultiVector{Data, GID, PID, LID}(data, numVecs, map)
end

 function DenseMultiVector(map::BlockMap{GID, PID, LID}, data::AbstractArray{Data, 2}) where {Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer}
    localLength = numMyElements(map)
    if size(data, 1) != localLength
         throw(InvalidArgumentError("Length of vectors does not match local length indicated by map"))
    end
    DenseMultiVector{Data, GID, PID, LID}(data, size(data, 2), map)
end

and then how I construct it to run the update function:

function update(a, A::MultiVector, b, B::MultiVector)
    if numVectors(A) == numVectors(B) && localLength(A) == localLength(B)
        myComm = SerialComm{Int, Int, Int}()
        myMap = BlockMap(numVectors(A), localLength(A), myComm)
        data = Array{Float32}(undef, localLength(A), numVectors(A))
        result = DenseMultiVector{}(myMap, data)

        A = scale!(A, a)
        B = scale!(B, b)

        for i = 1:numVectors(A)
            for j = 1:localLength(A)
                result[i,j] = A[i,j] + B[i,j]
            end
        end
     else
        println("Error: MultiVectors must be the same size")
    end
    return result
end
v1 = DenseMultiVector(curMap, ones(Float64, 2, 2))
v2 = DenseMultiVector(curMap, ones(Float64, 2, 2))
println("Update: ", update(2, v1, 3, v2))

What I expect is for the result to print [5.0 5.0; 5.0 5.0] but what is printed is [5.0 0.0; 5.0 0.0]. And when I do @show result[i,j] after it’s assignment, it prints: 5.0 0.0 5.0 0.0, just as is shown in the result.

7

Both of the methods you supplied for DenseMultiVector end with a call to what looks like a constructor method to the DenseMultiVector type, which would be defined in your package with a struct or mutable struct statement. Could you also please post that code? Also, since you are able to index into a DenseMultiVector there must be code that implements the getindex and setindex! methods for this type. Please post these as well.

1 Like
8

Here’s what I have for what you mentioned - thank you.

mutable struct DenseMultiVector{Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer} <: MultiVector{Data, GID, PID, LID}
    data::Array{Data, 2}
    numVectors::LID

    map::BlockMap{GID, PID, LID}
end


function Base.getindex(A::MultiVector, row::Integer, col::Integer)
    @boundscheck begin
        if !(1<=col<=numVectors(A))
            throw(BoundsError(A, (row, col)))
        end
    end

    lRow = lid(getMap(A), row)

    @boundscheck begin
        if lRow < 1
            throw(BoundsError(A, (row, col)))
        end
    end

    @inbounds value = getLocalArray(A)[lRow, col]
    value
end

function Base.getindex(A::MultiVector, i::Integer)
    if numVectors(A) != 1
        throw(ArgumentError("Can only use single index if there is just 1 vector"))
    end

    lRow = lid(getMap(A), i)

    @boundscheck begin
        if lRow < 1
            throw(BoundsError(A, I))
        end
    end

    @inbounds value = getLocalArray(A)[lRow, 1]
    value
end

function Base.setindex!(A::MultiVector, v, row::Integer, col::Integer)
    @boundscheck begin
        if !(1<=col<=numVectors(A))
            throw(BoundsError(A, (row, col)))
        end
    end

    lRow = lid(getMap(A), row)

    @boundscheck begin
        if lRow < 1
            throw(BoundsError(A, (row, col)))
        end
    end

    @inbounds getLocalArray(A)[lRow, 1] = v
    v
end

function Base.setindex!(A::MultiVector, v, i::Integer)
    if numVectors(A) != 0
        throw(ArgumentError("Can only use single index if there is just 1 vector"))
    end

    lRow = lid(getMap(A), i)

    @boundscheck begin
        if lRow < 1
            throw(BoundsError(A, I))
        end
    end

    @inbounds getLocalArray(A)[lRow, 1] = v
    v
end
9

Thanks. I can guess at the purpose of some of the functions referenced above, but that would be dangerous, especially when we’re on a bug hunt. So please also provide the codes that define numVectors, lid, getMap, getLocalArray, BlockMap, and localLength.

Note: This tedious back and forth is why it would be much easier to work with a complete MWE.

Note 2: I won’t be able to look at this further for a few hours.

10

With respect to the getindex/setindex!, methods, what is the purpose of doing manual bounds checks, and the call

?

Why not just rely on the automatic bounds checks? Doing manual checks seems like a risk of creating exactly the sort of weird bugs you are seeing.

11 
function numVectors end
function getLocalArray end

numVectors(mVect::DenseMultiVector) = mVect.numVectors

getMap(mVect::DenseMultiVector) = mVect.map

getLocalArray(mVect::DenseMultiVector) = mVect.data

@inline function lid(map::BlockMap{GID, PID, LID}, gid::Integer) where GID <: Integer where PID <: Integer where LID <: Integer
    data = map.data
        if (gid < data.minMyGID) || (gid > data.maxMyGID)
             LID(0)
        elseif data.linearMap
            LID(gid - data.minMyGID + 1)
        elseif gid >= data.myGlobalElements[1] && gid <= data.lastContiguousGID
            LID(gid - data.myGlobalElements[1] + 1)
        elseif haskey(data.lidHash, GID(gid))
            data.lidHash[gid]
        else
            LID(0)
        end
    end

struct BlockMap{GID <: Integer, PID <:Integer, LID <: Integer}
    data::BlockMapData{GID, PID, LID}

    function BlockMap{GID, PID, LID}(data::BlockMapData) where {GID <: Integer, PID <:Integer, LID <: Integer}
        new(data)
    end
end

**I also got rid of the bounds checks in getIndex! and setIndex! but the problem persists.

(Also, what would be the best way for me to share all the code needed to solve this issue?)

12

Thanks, I’ll look at this as soon as I’m able. If others want to jump in, they’re welcome as far as I’m concerned.

An easy way to share a large amount of code is a gist, or if you need to retain the directory structure of many files, just a standard Github repository.

13

I’m afraid that the BlockMap definition makes use of a so far undefined type BlockMapData which is also needed.

@fsmith001 : We could save some time by having you perform the same steps that I’ve been doing. As you provide each new block of code, I’ve been accumulating them in a single file, with the update function and its invocation positioned at the end of the file. I then try include-ing this file into a running Julia session, which quickly shows me what is still required for the code to compile and execute. If you could do this on your end and then post the file somewhere (such as in a gist or repository as I mentioned in my previous post), we could likely find your problem pretty quickly.

14

Okay, thank you! I did as you suggested, here is the link to the gist: trying to make the update() function work properly · GitHub

I think I see now where the problem is, but am not sure how to fix it. It seems that BlockMapData is not properly used based on how it is defined, however I do not know what exactly the issue is and how to use it correctly. Does this make sense?

15

Let’s continue the discussion on the gist.