search

OpenMendel / SnpArrays.jl

Compressed storage for SNP data
https://openmendel.github.io/SnpArrays.jl/latest
Other
44 stars 9 forks source link

A*x not correct if size(A, 1) % 4096 is 1, 2, or 3 #90

Closed biona001 closed 3 years ago

biona001 commented 3 years ago

MWE 

using Revise
using SnpArrays
using LinearAlgebra
using Random
using LoopVectorization
using MendelIHT
using Test

# any n between 4097 and 4099 doesn't work!
n = 4097
p = 10000
x = simulate_random_snparray(undef, n, p) #defined in my package `MendelIHT.jl`

A = SnpLinAlg{Float64}(x, model=ADDITIVE_MODEL, impute=false, center=false, scale=false)
b = ones(p)
c = A * b
ctrue = convert(Matrix{Float64}, A) * b
@test all(c .≈ ctrue) # fails!

The above fails. It seems Mrem at this line is the culprit

julia> [c ctrue]
4097×2 Matrix{Float64}:
 4930.0  4930.0
 4929.0  4929.0
 4929.0  4929.0
 5017.0  5017.0
 5044.0  5044.0
 5082.0  5082.0
 5112.0  5112.0
 4932.0  4932.0
 5055.0  5055.0
 4871.0  4871.0
 5009.0  5009.0
 5031.0  5031.0
 4960.0  4960.0
 4988.0  4988.0
 5001.0  5001.0
 5011.0  5011.0
 4969.0  4969.0
    ⋮
 4970.0  4970.0
 5001.0  5001.0
 4959.0  4959.0
 5065.0  5065.0
 5018.0  5018.0
 5037.0  5037.0
 5024.0  5024.0
 4987.0  4987.0
 4941.0  4941.0
 5016.0  5016.0
 5094.0  5094.0
 4901.0  4901.0
 5071.0  5071.0
 5045.0  5045.0
 5040.0  5040.0
 4984.0  4984.0
    0.0  4920.0