Но к сравнениям это отношения не имеет

Да, проблема в том, что матрицы по размеру разные. И в общем они временные. Сейчас внутри цикла создается матрица, заполняется нулями, потом с ней что-то происходит. И так каждую итерацию. Идея была - экономить память - создать перед циклом матрицу максимально требуемого размера и из нее брать view() и с ним работать - не работает. Думал может с view что-то не так. Проверил: внутри цикла создал матрицу, сделал ее view() - с ним дальше работал - проблем нет.

Из описания, всё же, не ясно, что конкретно делается. Если если минимальный пример и он не работает, то можем разобрать здесь 🙂

Пример огромный - очень большая функция, вот так работает: function reml_sweep_β(lmm, data::AbstractLMMDataBlocks, θ::Vector{T}) where T <: Number
   noerrors      = true
   n             = length(lmm.covstr.vcovblock)
   N             = length(lmm.data.yv)
   c             = (N - lmm.rankx)*log(2π)
   #---------------------------------------------------------------------------
   V⁻¹           = Vector{AbstractArray{T}}(undef, n)
   θ₁            = zero(T)
   θ₂            = zeros(T, lmm.rankx, lmm.rankx)
   θ₃            = zero(T)
   βm            = zeros(T, lmm.rankx)
   β             = Vector{T}(undef, lmm.rankx)
   #---------------------------------------------------------------------------
   local logdetθ₂::T
   akk           = Vector{T}(undef, lmm.covstr.maxn + lmm.rankx) #temp for sweep
   #Vm            = Matrix{T}(undef, lmm.covstr.maxn + lmm.rankx, lmm.covstr.maxn + lmm.rankx) #!!
   for i = 1:n
       q    = length(lmm.covstr.vcovblock[i])
       qswm = q + lmm.rankx
       Vpm  = Matrix{T}(undef, qswm, qswm) #!!
       Vp  = view(Vpm, 1:qswm, 1:qswm) #!!
       fill!(Vp, zero(T))
       V   = view(Vp, 1:q, 1:q)
       Vx  = view(Vp, 1:q, q+1:qswm)
       Vc  = view(Vp, q + 1:qswm, q + 1:qswm)
       copyto!(Vx, data.xv[i])
       vmatrix!(V, θ, lmm, i)
       #-----------------------------------------------------------------------
       try
           swr = sweepb!(fill!(view(akk, 1:qswm), zero(T)), Vp, 1:q; logdet = true)
           θ₁ += swr[2]
       catch
           noerrors = false
           lmmlog!(lmm, LMMLogMsg(:ERROR, "θ₁ not estimated during REML calculation, V isn't positive definite or |V| less zero."))
           return (1e100, nothing, nothing, 1e100, noerrors)
       end
       V⁻¹[i] = V
       #-----------------------------------------------------------------------
       subutri!(θ₂, Vc)
       #θ₂ -= UpperTriangular(view(Vp, q + 1:qswm, q + 1:qswm))
       mulαtβinc!(βm, Vx, data.yv[i])
       #-----------------------------------------------------------------------
   end
   θs₂ = Symmetric(θ₂)
   mul!(β, inv(θs₂), βm)
   for i = 1:n
       θ₃   += mulθ₃(data.yv[i], data.xv[i], β, V⁻¹[i])
   end
   try
       logdetθ₂ = logdet(θs₂)
   catch
       noerrors = false
       lmmlog!(lmm, LMMLogMsg(:ERROR, "logdet(θ₂) not estimated during REML calculation"))
       return (1e100, nothing, nothing, 1e100, noerrors)
   end
   return   θ₁ + logdetθ₂ + θ₃ + c, β, θs₂, θ₃, noerrors #REML, β, iC, θ₃, errors
end

вот так не работает: function reml_sweep_β(lmm, data::AbstractLMMDataBlocks, θ::Vector{T}) where T <: Number
   noerrors      = true
   n             = length(lmm.covstr.vcovblock)
   N             = length(lmm.data.yv)
   c             = (N - lmm.rankx)*log(2π)
   #---------------------------------------------------------------------------
   V⁻¹           = Vector{AbstractArray{T}}(undef, n)
   θ₁            = zero(T)
   θ₂            = zeros(T, lmm.rankx, lmm.rankx)
   θ₃            = zero(T)
   βm            = zeros(T, lmm.rankx)
   β             = Vector{T}(undef, lmm.rankx)
   #---------------------------------------------------------------------------
   local logdetθ₂::T
   akk           = Vector{T}(undef, lmm.covstr.maxn + lmm.rankx) #temp for sweep
   Vm            = Matrix{T}(undef, lmm.covstr.maxn + lmm.rankx, lmm.covstr.maxn + lmm.rankx) #!!
   for i = 1:n
       q    = length(lmm.covstr.vcovblock[i])
       qswm = q + lmm.rankx
       #Vpm  = Matrix{T}(undef, qswm, qswm) #!!
       Vp  = view(Vm, 1:qswm, 1:qswm) #!!
       fill!(Vp, zero(T))
       V   = view(Vp, 1:q, 1:q)
       Vx  = view(Vp, 1:q, q+1:qswm)
       Vc  = view(Vp, q + 1:qswm, q + 1:qswm)
       copyto!(Vx, data.xv[i])
       vmatrix!(V, θ, lmm, i)
       #-----------------------------------------------------------------------
       try
           swr = sweepb!(fill!(view(akk, 1:qswm), zero(T)), Vp, 1:q; logdet = true)
           θ₁ += swr[2]
       catch
           noerrors = false
           lmmlog!(lmm, LMMLogMsg(:ERROR, "θ₁ not estimated during REML calculation, V isn't positive definite or |V| less zero."))
           return (1e100, nothing, nothing, 1e100, noerrors)
       end
       V⁻¹[i] = V
       #-----------------------------------------------------------------------
       subutri!(θ₂, Vc)
       #θ₂ -= UpperTriangular(view(Vp, q + 1:qswm, q + 1:qswm))
       mulαtβinc!(βm, Vx, data.yv[i])
       #-----------------------------------------------------------------------
   end
   θs₂ = Symmetric(θ₂)
   mul!(β, inv(θs₂), βm)
   for i = 1:n
       θ₃   += mulθ₃(data.yv[i], data.xv[i], β, V⁻¹[i])
   end
   try
       logdetθ₂ = logdet(θs₂)
   catch
       noerrors = false
       lmmlog!(lmm, LMMLogMsg(:ERROR, "logdet(θ₂) not estimated during REML calculation"))
       return (1e100, nothing, nothing, 1e100, noerrors)
   end
   return   θ₁ + logdetθ₂ + θ₃ + c, β, θs₂, θ₃, noerrors #REML, β, iC, θ₃, errors
end

Пара замечаний:
1. Такие огромные штуки лучше через gist выкладывать.

2. А возможно от него куски отрезать, пока минимальный кусок не останется, в котором баг воспроизводится?

Ведь нет необходимости, чтобы какой-то вменяемый результат получился, нужно, чтобы ошибка воспроизводилась.

Ещё по коду:

        Vp  = view(Vm, 1:qswm, 1:qswm) #!!
        fill!(Vp, zero(T))
        V   = view(Vp, 1:q, 1:q)
        Vx  = view(Vp, 1:q, q+1:qswm)
        Vc  = view(Vp, q + 1:qswm, q + 1:qswm)

Есть begin/end чтобы не ошибиться в граничных индексах

сейчас попробую

https://gist.github.com/PharmCat/01f43a0f3795ede52a3a940c78a49e01#file-gistfile1-txt вот, к примеру

Я сходу ошибку не вижу. Но что бы сделал - распечатку подозрительной матрицы или контрольных сумм с неё на каждой итерации. Потенциально, что-то не то с индексами

еще уменьшил... индексов не осталось почти... при чем 1я различается, вторая нет...

V это view ведь? ну типа не может такое случиться, что Vr[1] присвоили

а потом когда V меняли, то Vr[1] тоже поменялась

А вообще Vr[1] == Vr[2] ?

я бы положил так, если бы сам писал: Vr[i] = deepcopy(V)

Все! я понял где я неправ.   Действительно Vr[i] = deepcopy(V), без этого не обойтись видимо... Большое всем спасибо!

Ну не знаю, не знаю.

Выглядит подозрительно.