Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP: typestable factorizations #9575

Merged
merged 1 commit into from
Jan 28, 2015
Merged

WIP: typestable factorizations #9575

merged 1 commit into from
Jan 28, 2015

Conversation

skariel
Copy link
Contributor

@skariel skariel commented Jan 3, 2015

There are several more factorization functions that can be made type-stable. If this is positively accepted I'l gladly change them too :)

@skariel skariel changed the title initial changes to qrfact typestable qrfact Jan 3, 2015
@andreasnoack
Copy link
Member

👍 However, we need to figure out if we want use types or instances for this, see #9452. I think I prefer the type version.

stagedfunction qrfact{T}(A::StridedMatrix{T}, ::Pivot)
S = typeof(one(T)/norm(one(T)))
S != T ? :(qrfact!(convert(AbstractMatrix{S},A),pivot)) : :(qrfact!(copy(A)),pivot)
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is a stagedfunction really necessary here? Couldn't you do:

copy_oftype{T}(A::AbstractArray{T}, ty::Type{T}) = copy(A)
copy_oftype{T,N}(A::AbstractArray{T,N}, ty) = convert(AbstractArray{ty,N}, A)
qrfact{T}(A::StridedMatrix{T}) = qrfact!(copy_oftype(A, typeof(one(T)/norm(one(T))))

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The solution we find better should be used for all the promotion methods in the linear algebra code.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since this is likely to be useful outside of Base, perhaps we could have:

copy{T}(::Type{T}, x::T) = copy(x)
copy{T}(::Type{T}, x) = convert(T, x)

Although the second method may be too broad, since convert to a new type is not guaranteed to copy everything that copy would.

@andreasnoack
Copy link
Member

As I read the discussion, we should use the type version. It would be great if you could also prepare a similar change for the Cholesky as it also has a pivoted version.

@ViralBShah
Copy link
Member

Before merging - it would be nice to squash the commits, and have an explanation on type stability in the commit message.

I faintly recollect a discussion on the mailing list too, but couldn't find it. If so, could someone who knows, link it here?

@skariel
Copy link
Contributor Author

skariel commented Jan 5, 2015

Sure, I'll do the following:

  1. use copy_oftype instead of stagedfunction
  2. type stability also for Cholesky factorization
  3. squash the commits into a single one

There is some discussion on this type instability in the docs, is this what you mean, see here:
http://julia.readthedocs.org/en/latest/stdlib/linalg/?highlight=qrfact#Base.qrfact
(see note in blue below...)

@andreasnoack
Copy link
Member

@skariel Great. I like @simonster's suggestion of extending copy to allow a type specification.

I think @ViralBShah is thinking of this discussion

@skariel skariel changed the title typestable qrfact WIP: typestable qrfact Jan 5, 2015
@skariel skariel changed the title WIP: typestable qrfact WIP: typestable factorizations Jan 5, 2015
@@ -538,8 +540,8 @@ Deprecated or removed
argument specifying the floating point type to which they apply. The old
behaviour and `[get/set/with]_bigfloat_rounding` functions are deprecated ([#5007])

* `cholpfact` and `qrpfact` are deprecated in favor of keyword arguments in
`cholfact(..., pivot=true)` and `qrfact(..., pivot=true)` ([#5330])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is in the Julia 0.3 release notes. I think we need a separate entry for these changes under Julia 0.4. Also it would be good to have deprecations in deprecated.jl.

@skariel
Copy link
Contributor Author

skariel commented Jan 6, 2015

lufact has a reversed default pivot argument ie pivot=true. So separating this into types means defining a type NoPivot and a cooresponging constant nopivot. Is this the way to go? the other options would be to use a different default, pivot=false like other functions do. But then this a a change in api that would require many changes all over the place. I actually started doing this but maybe I'll just revert.

What do you think?

@skariel
Copy link
Contributor Author

skariel commented Jan 6, 2015

there is no need actually for copy_oftype since there is no dependence of return type on value for this. So this is not a type instability. Only interesting are cases where return type depends on value of pivot. I'm changing this now

@simonster
Copy link
Member

@skariel I'm pretty sure you need copy_oftype because Julia still considers the branch that's not taken in type inference (even though it could be inferred to be dead; see #5560), and if S != T the two branches give different types, e.g.:

julia> f{T}(A::StridedMatrix{T}) = (S = typeof(one(T)/norm(one(T)));S != T ? convert(AbstractMatrix{S},A) : copy(A));

julia> code_typed(f, (Matrix{Int},))
1-element Array{Any,1}:
 :($(Expr(:lambda, Any[:A], Any[Any[:S,:_var0,:_var1,:_var2,:_var3,:_var4,:_var6,:_var7,:_var9,:_var10,:_var11,:_var12],Any[Any[:A,Array{Int64,2},0],Any[:S,Type{Float64},18],Any[:_var0,Int64,18],Any[:_var1,Bool,18],Any[:_var2,Type{Int64},18],Any[:_var3,Type{Float64},18],Any[:_var4,Int64,18],Any[:_var6,Int64,18],Any[:_var7,Int64,18],Any[:_var9,Int64,18],Any[:_var10,Int64,18],Any[:_var11,Int64,18],Any[:_var12,Array{Int64,2},18]],Any[]], :(begin  # none, line 1:
        _var0 = norm(1)::Int64
        S = typeof((top(box))(Float64,(top(div_float))((top(box))(Float64,(top(sitofp))(Float64,1)),(top(box))(Float64,(top(sitofp))(Float64,_var0::Int64)))))::Type{Float64}
        _var3 = S::Type{Float64}
        _var2 = T::Any
        _var1 = _var3::Type{Float64} == _var2::Type{Int64}::Bool
        unless (top(box))(Bool,(top(not_int))(_var1::Bool)) goto 0
        _var4 = (top(arraysize))(A::Array{Int64,2},1)::Int64
        _var7 = _var4::Int64
        _var6 = (top(arraysize))(A::Array{Int64,2},2)::Int64
        return copy!((top(ccall))(:jl_alloc_array_2d,$(Expr(:call1, :(top(apply_type)), :Array, Float64, 2)),$(Expr(:call1, :(top(tuple)), :Any, :Int, :Int)),Array{Float64,2},0,_var7::Int64,0,_var6::Int64,0)::Array{Float64,2},A::Array{Int64,2})::Array{Float64,2}
        0: 
        _var9 = (top(arraysize))(A::Array{Int64,2},1)::Int64
        _var11 = _var9::Int64
        _var10 = (top(arraysize))(A::Array{Int64,2},2)::Int64
        _var12 = (top(ccall))(:jl_alloc_array_2d,$(Expr(:call1, :(top(apply_type)), :Array, Int64, 2)),$(Expr(:call1, :(top(tuple)), :Any, :Int, :Int)),Array{Int64,2},0,_var11::Int64,0,_var10::Int64,0)::Array{Int64,2}
        return copy!(_var12::Array{Int64,2},1,A::Array{Int64,2},1,(top(arraylen))(A::Array{Int64,2})::Int64)::Array{Int64,2}
    end::Union(Array{Int64,2},Array{Float64,2})))))

@simonster
Copy link
Member

We could define !(::Pivot) = NoPivot() and use !pivot instead of nopivot.

@skariel
Copy link
Contributor Author

skariel commented Jan 6, 2015

yes, I see your point. I like the !pivot suggestion, lets use this

@skariel
Copy link
Contributor Author

skariel commented Jan 6, 2015

About copy_oftype -- it won't work for more elaborate cases I think, like this one for e.g.:

function cholfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0)
    S = promote_type(typeof(chol(one(T))),Float32)
    S <: BlasFloat && return cholfact!(convert(AbstractMatrix{S}, A), uplo, pivot = pivot, tol = tol)
    pivot && throw(ArgumentError("pivot only supported for Float32, Float64, Complex{Float32} and Complex{Float64}"))
    S != T && return cholfact!(convert(AbstractMatrix{S}, A), uplo)
    return cholfact!(copy(A), uplo)
end

Here (unless I'm mistaken of course :) only a staged function would help. So maybe better to use them for these cases?

EDIT -- I mean, maybe better to handle all cases equally ie using a stagedfunction

@andreasnoack
Copy link
Member

@skariel I think the cholfact example you gave would go into at least two new definitions when using the Pivot type instead of a keyword. Therefore I think than @simonster's suggestion is still enough, but I might be wrong.

In #9452, the preference seemed to be for the type version, so I think we should use Pivot instead of the pivot instance.

Would Pivot{true} and Pivot{false} be a better option here?

@simonster
Copy link
Member

@skariel Adopting @andreasnoack's suggestion of Pivot{false} and Pivot{true}, I would handle that as:

cholfact!{T}(A::StridedMatrix{T}, uplo::Symbol=:U, ::Type{Pivot{false}}) = # non-pivoting alg goes here
cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U, ::Type{Pivot{false}}) = # pivoting alg goes here
cholfact!{T}(A::StridedMatrix{T}, uplo::Symbol=:U, ::Type{Pivot{false}}) =
    throw(ArgumentError("pivot only supported for Float32, Float64, Complex{Float32} and Complex{Float64}"))
cholfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U, pivot::Union(Type{Pivot{false}},Type{Pivot{true}})=Pivot{false}) =
    cholfact!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), uplo, pivot)

# probably also want this one...
cholfact!{T}(A::StridedMatrix{T}, uplo::Symbol=:U) = cholfact!(A, uplo, Pivot{false})

which gives the same error message for types for which pivoting is unsupported for both cholfact and cholfact!.

[edit: fixed method signature for cholfact to use a union type]

@skariel
Copy link
Contributor Author

skariel commented Jan 6, 2015

@simonster, I don't think your version is equivalent to the original one. For e.g. if using T==Type{BigInt} and pivot==false then the original version would dispatch cholfact!(convert(AbstractMatrix{S}, A), uplo) while if T==Type{Float32} it would dipatch cholfact!(convert(AbstractMatrix{S}, A), uplo, tol = tol). Watch the tolerance parameter, there is no distinction betwen the two cases in your version...

(in the T==Type{BigInt} case s==Type{BigFloat})

@simonster
Copy link
Member

@skariel You can define:

cholfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U, pivot::Union(Type{Pivot{false}},Type{Pivot{true}})=Pivot{false}; tol=0.0) =
    cholfact!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), uplo, pivot; tol=tol)
cholfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U, pivot::Union(Type{Pivot{false}},Type{Pivot{true}})=Pivot{false}) =
    cholfact!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), uplo, pivot)

Alternatively, everything could accept a tol parameter and cholfact! for non-BlasFloats could warn/error if it's non-zero.

@skariel
Copy link
Contributor Author

skariel commented Jan 7, 2015

@simonster the behavior is still not the same:

consider the original version given
pivot==false,
T==Type{Rational{Int64}} and
tol>0
it dispatches cholfact!(convert(AbstractMatrix{S}, A), uplo, pivot = pivot, tol = tol)
where S==Type{Float32}.

In your version, since Rational{Int64} is not a subtype of BlasFloat this would dispatch without specifying tol. This would cause a hard to catch bug.

This function has a complex dispatch on types which is hard (or impossible) to express with current multiple-dispatch semantics, this is likly why it was implemented in a type-unstable way before stagedfunctions were available

@simonster
Copy link
Member

You're right. But I believe you can still get the current behavior by dispatch:

cholfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U, pivot::Union(Type{Pivot{false}},Type{Pivot{true}})=Pivot{false}; tol=0.0) =
    _cholfact!(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), uplo, pivot, tol)
_cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo, pivot, tol) = cholfact!(A, uplo, pivot; tol=tol)
_cholfact!(A, uplo, pivot, tol) = cholfact!(A, uplo, pivot)

Obviously it would be good to have test coverage for this to make sure it actually behaves properly. And we may still want to do something besides silently drop non-zero tol parameters for non-BlasFloats.

I'm not sure using a staged function would actually make this much simpler. @andreasnoack would know more about the history of the current implementation.

@skariel
Copy link
Contributor Author

skariel commented Jan 7, 2015

Yes, that seems to work! I guess I'm still getting used to multiple dispatch and just didn't think of breaking a single function into three levels of dispatch. But it makes sense, it's equivalent to if equals else on types and that's all we have here. So I'll use this of course

@andreasnoack
Copy link
Member

@skariel Right now when pivot=false, the tol argument is just ignored at a later point. I don't think it is a huge problem, but having a way of throwing an error when the user it trying to set tol for normal Cholesky could be a possibility.

The reason the dispatch is so convolved right now is partly that the tolerance only makes sense for pivoted Cholesky and partly that non pivoted Cholesky has a generic fall back and pivoted Cholesky is LAPACK types only.

There also used to be a separate cholpfact method for the pivoted version, but we decided to clean up the namespace a bit and only later realized that the type stability issue was too costly.

@ViralBShah ViralBShah added the linear algebra Linear algebra label Jan 8, 2015
@andreasnoack
Copy link
Member

@skariel Bump.

@skariel
Copy link
Contributor Author

skariel commented Jan 14, 2015

I'm here, I'll have time tomorrow, will try to push then

@skariel
Copy link
Contributor Author

skariel commented Jan 20, 2015

@andreasnoack I think this is it, all factorizations should be stable, all tests are passing. Is this solution good enough? If so ill start updating the metadata (documentation, news, etc.)

@skariel
Copy link
Contributor Author

skariel commented Jan 20, 2015

@andreasnoack sorry false alarm... I found a few things that still need doing. I'll try to finish it tomorrow

@skariel
Copy link
Contributor Author

skariel commented Jan 21, 2015

@andreasnoack it's good now, can you please check the interface is what we want, If so I'll just update the docs and that should be it

Thanks!

@skariel
Copy link
Contributor Author

skariel commented Jan 26, 2015

I think using Val{T} is good because it is a general solution. I just checked and having pivot a keyword argument seems no problem for type stability see attached image. So I suggest just changing Pivot{T} to Val{T} and merge. One open question: how should the union Val{ture} and Val{false} be named?

screenshot from 2015-01-26 08 00 26

@skariel skariel closed this Jan 26, 2015
@skariel
Copy link
Contributor Author

skariel commented Jan 26, 2015

(closed by mistake)

@skariel skariel reopened this Jan 26, 2015
@simonster
Copy link
Member

@skariel That test case has a default argument (which is equivalent to an ordinary positional argument for type inference purposes), not a keyword argument, and the call to @code_typed is within the function rather than outside of it. This is not type-stable:

julia> f(; n::Union(Float64,Int)=1) = 1+n;

julia> g(x) = f(n=x);

julia> @code_typed g(1)
1-element Array{Any,1}:
 :($(Expr(:lambda, Any[:x], Any[Any[],Any[Any[:x,Int64,0]],Any[]], :(begin  # none, line 1:
        return (top(kwcall))(call,1,:n,x::Int64,f,(top(ccall))(:jl_alloc_array_1d,$(Expr(:call1, :(top(apply_type)), :Array, Any, 1)),$(Expr(:call1, :(top(tuple)), :Any, :Int)),Array{Any,1},0,2,0)::Array{Any,1})::Union(Float64,Int64)
    end::Union(Float64,Int64)))))

Also, I tried building this branch and I see:

julia> Base.return_types(qrfact, (Matrix{Float64},))
1-element Array{Any,1}:
 Union(Base.LinAlg.QRCompactWY{Float64,Array{Float64,2}},Base.LinAlg.QRPivoted{Float64,Array{Float64,2}})

julia> Base.return_types(cholfact, (Matrix{Float64},))
1-element Array{Any,1}:
 Union(Base.LinAlg.Cholesky{Float64,S<:AbstractArray{T,2},UpLo},Base.LinAlg.CholeskyPivoted{Float64,S<:AbstractArray{T,2}})

_cholfact(copy_oftype(A, promote_type(typeof(chol(one(T))),Float32)), pivot, uplo, tol=tol)

_cholfact{T<:BlasFloat}(A::StridedMatrix{T}, pivot::PivotUnion, uplo::Symbol=:U; tol=0.0) =
cholfact!(A, uplo, pivot = pivot, tol = tol)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems this needs to be split into two identical methods for the different Pivot arguments to avoid a method ambiguity warning:

_cholfact{T<:BlasFloat}(A::StridedMatrix{T}, pivot::Type{Pivot{true}}, uplo::Symbol=:U; tol=0.0) =
    cholfact!(A, uplo, pivot = pivot, tol = tol)
_cholfact{T<:BlasFloat}(A::StridedMatrix{T}, pivot::Type{Pivot{false}}, uplo::Symbol=:U; tol=0.0) =
    cholfact!(A, uplo, pivot = pivot, tol = tol)

@skariel
Copy link
Contributor Author

skariel commented Jan 26, 2015

yep you're right. I'll change Pivot{T} to Val{T}, and the pivot keyword argument to a positional one. Also I'll drop the PivotUnion in favor of an in-place union this way we get one thing less to name

@skariel
Copy link
Contributor Author

skariel commented Jan 27, 2015

more type instability:

function (\)(A::StridedMatrix, B::StridedVecOrMat)
    m, n = size(A)
    if m == n
        if istril(A)
            return istriu(A) ? \(Diagonal(A),B) : \(LowerTriangular(A),B)
        end
        istriu(A) && return \(UpperTriangular(A),B)
        return \(lufact(A),B)
    end
    return qrfact(A,eltype(A)<:BlasFloat?Val{true}:Val{false})\B
end

@skariel
Copy link
Contributor Author

skariel commented Jan 27, 2015

I mean, all this can be fixed, why even allow this? is there a single instance in the standard lib where type instability is actually needed? I'd like to see that

@skariel
Copy link
Contributor Author

skariel commented Jan 27, 2015

Tests seem to fail because of a bug, I've filed a new issue #9938

@skariel
Copy link
Contributor Author

skariel commented Jan 27, 2015

@simonster this is ready. I'm just waiting for #9938 to squash

@StefanKarpinski
Copy link
Member

@skariel: I mean, all this can be fixed, why even allow this? is there a single instance in the standard lib where type instability is actually needed? I'd like to see that

Why is Julia a dynamic language rather than a static one? Kind of a deep question to sneak into a comment on a pull request. Largely because dynamic languages are easier to use. A static language assigns types to all expressions based on the structure of your program. This means that there are rules that are part of the language – pretty complicated ones, inevitably – for how types get assigned to expressions. If you write a program that doesn't obey those rules, then it won't compile, let alone run. This means that programmers have to learn those rules before they can write programs that will run. In dynamic languages, on the other hand, there are no such rules, which is a massive reduction in the learning curve before you can use the language. You can write a Julia program that is wildly type unstable and it will work. If you wonder why it's not as fast as you would like, then there's more to learn, but at least you can get a program working in the first place.

@StefanKarpinski
Copy link
Member

I also believe that \ is precisely one of those cases where type instability is intentional.

@andreasnoack
Copy link
Member

@skariel I thought for a long time that one of branches in S==T ? typeA : typeB were eliminated when possible which is part of the reason why there are so many type instabilities in the linalg code. The other reason is that, in some cases, we have decided to be type unstable, e.g. when we removed the cholpfact function in favor of keyword in cholfact we were fully aware of the type stability issue. I think we underestimated the performance impact of this so I this it is good that some if this is reverted by your fixes here.

There are a couple of problematic cases:

  • The factorize function has been mentioned a couple of times as an example of a function which cannot be type stable.
  • In cholmod.jl there will be a function that can read files in the MatrixMarket format. It dynamically determines if the elements are real or complex, the integer size needed and if the matrix is symmetric/Hermitian.
  • The eigen functions for non-symmetric problems. Here we could make it type stable, but it would require us to always return complex results which we might don't want to do.

However, I don't see a reason why Matrix{T}\VecOrMat{S} should be type unstable. Though, there is a type stability problem for SubArray rhs, which is something we'll have to consider quite generally in the linalg code.

I think it is important to be specific about for which types there is a type instability. Many Julia function have a very broad signature and therefore it will probably be type unstable for some weird input type that doesn't really matter.

Previously retured types depended on the boolean value of a keyword argument `pivot`. Now the api relies on dispatch on different tyeps either `pivot=Val{true}` or `pivot=Val{false}`. A few more type instabilities were mitigated by using a `copy_oftype` method.
@skariel
Copy link
Contributor Author

skariel commented Jan 27, 2015

Yeah I already had a long discussion of this in the forums. Sorry for all the noise :)

@StefanKarpinski I enjoy daily the interactive workflow that Julia enables
@andreasnoack Dynamically loading files is a stong point in favor of type-instability. It could be done type-stable at the price of a more convoluted solution say using some ifs etc. which might be OK for library developpers

Anyway this is done and ready to merge (no new warnings, test passing, etc.)

@skariel skariel closed this Jan 27, 2015
@skariel skariel reopened this Jan 27, 2015
@skariel
Copy link
Contributor Author

skariel commented Jan 27, 2015

(closed by mistake)

@timholy
Copy link
Member

timholy commented Jan 27, 2015

Though, there is a type stability problem for SubArray rhs, which is something we'll have to consider quite generally in the linalg code.

This sounds important; more details, please?

@andreasnoack
Copy link
Member

@timholy

julia> @code_warntype eye(3)\sub(ones(3), 1:3)
Variables:
  A::Array{Float64,2}
  B::SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1}
  #s1994::(Int64,Int64)
  #s1992::(Int64,Int64)
  m::Int64
  #s1991::(Int64,Int64)
  n::Int64
  #s1993::Int64
  _var2::Int64
  _var3::Array{Float64,1}
  _var0::Int64
  _var1::Int64
  _var4::Int64
  _var5::Int64
  _var6::Int64
  _var7::Int64

Body:
  begin  # linalg/dense.jl, line 377:
      _var0 = (top(arraysize))(A::Array{Float64,2},1)::Int64
      _var1 = (top(arraysize))(A::Array{Float64,2},2)::Int64
      #s1993 = 1
      _var4 = _var0::Int64
      _var5 = (top(box))(Int64,(top(add_int))(1,1))
      m = _var4::Int64
      #s1993 = _var5::Int64
      _var6 = _var1::Int64
      _var7 = (top(box))(Int64,(top(add_int))(2,1))
      n = _var6::Int64
      #s1993 = _var7::Int64 # line 378:
      unless m::Int64 === n::Int64::Bool goto 3 # line 379:
      unless istril(A::Array{Float64,2})::Bool goto 1 # line 380:
      unless istriu(A::Array{Float64,2})::Bool goto 0
      _var2 = (top(arraysize))(A::Array{Float64,2},1)::Int64
      _var3 = getindex(A::Array{Float64,2},diagind(_var2::Int64,(top(arraysize))(A::Array{Float64,2},2)::Int64,0)::StepRange{Int64,Int64})::Array{Float64,1}
      return $(Expr(:new, Base.LinAlg.Diagonal{Float64}, :(_var3::Array{Float64,1}))) \ B::SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1}::Union(Array{Float64,1},SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1})
      0: 
      ((top(getfield))((top(getfield))(Base,:LinAlg),:chksquare)::Any)(A::Array{Float64,2})::Int64
      return $(Expr(:new, Base.LinAlg.LowerTriangular{Float64,Array{Float64,2}}, :(A::Array{Float64,2}))) \ B::SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1}::Union(Array{Float64,1},SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1})
      1:  # line 382:
      unless istriu(A::Array{Float64,2})::Bool goto 2
      ((top(getfield))((top(getfield))(Base,:LinAlg),:chksquare)::Any)(A::Array{Float64,2})::Int64
      return $(Expr(:new, Base.LinAlg.UpperTriangular{Float64,Array{Float64,2}}, :(A::Array{Float64,2}))) \ B::SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1}::Union(Array{Float64,1},SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1})
      2:  # line 383:
      return __lufact#181__(true,A::Array{Float64,2})::Base.LinAlg.LU{Float64,Array{Float64,2}} \ B::SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1}::Union(Array{Float64,1},SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1})
      3:  # line 385:
      return (top(kwcall))(call,1,:pivot,Float64 <: BlasFloat::Bool,qrfact,(top(ccall))(:jl_alloc_array_1d,$(Expr(:call1, :(top(apply_type)), :Array, Any, 1)),$(Expr(:call1, :(top(tuple)), :Any, :Int)),Array{Any,1},0,2,0)::Array{Any,1},A::Array{Float64,2})::Union(Base.LinAlg.QRPivoted{Float64,Array{Float64,2}},Base.LinAlg.QRCompactWY{Float64,Array{Float64,2}}) \ B::SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1}::Union(Array{Float64,1},SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1})
  end::Union(Array{Float64,1},SubArray{Float64,1,Array{Float64,1},(UnitRange{Int64},),1})

@timholy
Copy link
Member

timholy commented Jan 27, 2015

Thanks, I understand better now. Not really an issue with SubArrays themselves (which is what I was worried about), but a genuine challenge in structuring the linear algebra code. For example, in this method the problem is it might return the SubArray A (if m or n is zero) or the Array C.
You could presumably always return C?

function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0)
cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U, pivot::Union(Type{Val{false}}, Type{Val{true}})=Val{false}; tol=0.0) =
_cholfact!(A, pivot, uplo, tol=tol)
function _cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, ::Type{Val{false}}, uplo::Symbol=:U; tol=0.0)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why are the underscore versions necessary? Couldn't the underscore versions just be used as the ordinary cholfact! methods instead of having cholfact! calling _cholfact!?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This would require changing the API as dispatching on Val{false} and Val{true} would require pivot to not be a default argument and hence change places with uplo. It would still require 3 functions, and since _cholfact are not exported anyway I preferred to not changed the API

What do you think?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah. I hadn't thought about that. Thanks for the explanation. Let's stick with your solution. I'll merge.

andreasnoack added a commit that referenced this pull request Jan 28, 2015
@andreasnoack andreasnoack merged commit f1973ac into JuliaLang:master Jan 28, 2015
@andreasnoack
Copy link
Member

@timholy No. It has nothing to do with SubArrays. It is about the linear algebra methods. I dont' think they have been widely used with SubArrays yet so the type problems haven't been revealed. We should begin to include @inferred tests in the linalg test suite.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants