BigFloat matrix in partialschur() does not work #94

sverek · 2019-01-11T17:52:13Z

MWE:

using ArnoldiMethod, LinearAlgebra, SparseArrays
function mwe()
A = spdiagm(
           -1 => -1.0*ones(9),
            0 =>  2.0*ones(10), 
            1 => -1.0*ones(9)
       );
B = spdiagm(
           -1 => -big(1.0)*ones(9),
            0 => big(2.0)*ones(10), 
            1 => -big(1.0)*ones(9)
       );
display(A)
decomp, history = partialschur(A, nev=10, tol=1e-6, which=SR());
display(decomp)

display(B)
decomp, history = partialschur(B, nev=10, tol=1e-6, which=SR());
display(decomp)
end

Output:


julia> mwe()
10×10 SparseMatrixCSC{Float64,Int64} with 28 stored entries:
  [1 ,  1]  =  2.0
  [2 ,  1]  =  -1.0
  [1 ,  2]  =  -1.0
  [2 ,  2]  =  2.0
  [3 ,  2]  =  -1.0
  [2 ,  3]  =  -1.0
  [3 ,  3]  =  2.0
  [4 ,  3]  =  -1.0
  [3 ,  4]  =  -1.0
  ⋮
  [8 ,  7]  =  -1.0
  [7 ,  8]  =  -1.0
  [8 ,  8]  =  2.0
  [9 ,  8]  =  -1.0
  [8 ,  9]  =  -1.0
  [9 ,  9]  =  2.0
  [10,  9]  =  -1.0
  [9 , 10]  =  -1.0
  [10, 10]  =  2.0
10×10 SparseMatrixCSC{BigFloat,Int64} with 28 stored entries:
  [1 ,  1]  =  2.0
  [2 ,  1]  =  -1.0
  [1 ,  2]  =  -1.0
  [2 ,  2]  =  2.0
  [3 ,  2]  =  -1.0
  [2 ,  3]  =  -1.0
  [3 ,  3]  =  2.0
  [4 ,  3]  =  -1.0
  [3 ,  4]  =  -1.0
  ⋮
  [8 ,  7]  =  -1.0
  [7 ,  8]  =  -1.0
  [8 ,  8]  =  2.0
  [9 ,  8]  =  -1.0
  [8 ,  9]  =  -1.0
  [9 ,  9]  =  2.0
  [10,  9]  =  -1.0
  [9 , 10]  =  -1.0
  [10, 10]  =  2.0
PartialSchur decomposition (Float64) of dimension 10
eigenvalues:
10-element Array{Complex{Float64},1}:
 0.08101405277100668 + 0.0im
 0.31749293433763753 + 0.0im
  0.6902785321094296 + 0.0im
   1.169169973996227 + 0.0im
   1.715370323453429 + 0.0im
   2.284629676546568 + 0.0im
   2.830830026003772 + 0.0im
  3.3097214678905678 + 0.0im
   3.682507065662358 + 0.0im
  3.9189859472289927 + 0.0im
ERROR: MethodError: no method matching gemv!(::Char, ::BigFloat, ::SubArray{BigFloat,2,Array{BigFloat,2},Tuple{Base.Slice{Base.OneTo{Int64}},UnitRange{Int64}},true}, ::SubArray{BigFloat,1,Array{BigFloat,2},Tuple{UnitRange{Int64},Int64},true}, ::BigFloat, ::SubArray{BigFloat,1,Array{BigFloat,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true})
Closest candidates are:
  gemv!(::AbstractChar, ::Float64, ::Union{AbstractArray{Float64,1}, AbstractArray{Float64,2}}, ::AbstractArray{Float64,1}, ::Float64, ::AbstractArray{Float64,1}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/blas.jl:565
  gemv!(::AbstractChar, ::Float32, ::Union{AbstractArray{Float32,1}, AbstractArray{Float32,2}}, ::AbstractArray{Float32,1}, ::Float32, ::AbstractArray{Float32,1}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/blas.jl:565
  gemv!(::AbstractChar, ::Complex{Float64}, ::Union{AbstractArray{Complex{Float64},1}, AbstractArray{Complex{Float64},2}}, ::AbstractArray{Complex{Float64},1}, ::Complex{Float64}, ::AbstractArray{Complex{Float64},1}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/blas.jl:565

The text was updated successfully, but these errors were encountered:

haampie · 2019-01-13T09:47:17Z

Thanks, I will look into that. There should be a fallback for this in Julia base when the mul!api is done

andreasnoack · 2019-01-14T08:22:17Z

Ref: JuliaLang/julia#29634 and
JuliaLang/LinearAlgebra.jl#473. Unfortunately, it seems that the discussion has dead locked.

jarlebring · 2019-01-19T09:10:10Z

Workaround if you really need it and don't care about performance so much: Write your own fallback.

import LinearAlgebra.BLAS.gemv!
function gemv!(C,alpha,A,X,beta,Y)
     if (C=='N')
        Y[:]=alpha*A*X+beta*Y
    elseif (C=='T')
        Y[:]=alpha*transpose(A)*X+beta*Y
    elseif (C=='C')
        Y[:]=alpha*A'*X+beta*Y
    else 
        error("Unknown transpose character")
    end
end

Output:

julia> mwe()
10×10 SparseMatrixCSC{Float64,Int64} with 28 stored entries:
  [1 ,  1]  =  2.0
  [2 ,  1]  =  -1.0
  [1 ,  2]  =  -1.0
  [2 ,  2]  =  2.0
  [3 ,  2]  =  -1.0
  [2 ,  3]  =  -1.0
  [3 ,  3]  =  2.0
  [4 ,  3]  =  -1.0
  [3 ,  4]  =  -1.0
  [4 ,  4]  =  2.0
  [5 ,  4]  =  -1.0
  [4 ,  5]  =  -1.0
  ⋮
  [7 ,  6]  =  -1.0
  [6 ,  7]  =  -1.0
  [7 ,  7]  =  2.0
  [8 ,  7]  =  -1.0
  [7 ,  8]  =  -1.0
  [8 ,  8]  =  2.0
  [9 ,  8]  =  -1.0
  [8 ,  9]  =  -1.0
  [9 ,  9]  =  2.0
  [10,  9]  =  -1.0
  [9 , 10]  =  -1.0
  [10, 10]  =  2.0
PartialSchur decomposition (Float64) of dimension 10
eigenvalues:
10-element Array{Complex{Float64},1}:
 0.08101405277100521 + 0.0im
  0.3174929343376377 + 0.0im
  0.6902785321094296 + 0.0im
  1.1691699739962276 + 0.0im
  1.7153703234534308 + 0.0im
   2.284629676546571 + 0.0im
   2.830830026003773 + 0.0im
  3.3097214678905718 + 0.0im
   3.682507065662364 + 0.0im
   3.918985947228997 + 0.0im
10×10 SparseMatrixCSC{BigFloat,Int64} with 28 stored entries:
  [1 ,  1]  =  2.0
  [2 ,  1]  =  -1.0
  [1 ,  2]  =  -1.0
  [2 ,  2]  =  2.0
  [3 ,  2]  =  -1.0
  [2 ,  3]  =  -1.0
  [3 ,  3]  =  2.0
  [4 ,  3]  =  -1.0
  [3 ,  4]  =  -1.0
  [4 ,  4]  =  2.0
  [5 ,  4]  =  -1.0
  [4 ,  5]  =  -1.0
  ⋮
  [7 ,  6]  =  -1.0
  [6 ,  7]  =  -1.0
  [7 ,  7]  =  2.0
  [8 ,  7]  =  -1.0
  [7 ,  8]  =  -1.0
  [8 ,  8]  =  2.0
  [9 ,  8]  =  -1.0
  [8 ,  9]  =  -1.0
  [9 ,  9]  =  2.0
  [10,  9]  =  -1.0
  [9 , 10]  =  -1.0
  [10, 10]  =  2.0
PartialSchur decomposition (BigFloat) of dimension 10
eigenvalues:
10-element Array{Complex{BigFloat},1}:
 8.101405277100522021926388586734460187509030315567608991100988200961938340741451e-02 + 0.0im
 3.174929343376376622763767021612645649734150031589242027146997646366977861967341e-01 + 0.0im
 6.902785321094298718861498550674128936324176013261431447874087909882098201012722e-01 + 0.0im
     1.169169973996227148941451701540753592951990179070926375147553335243120438110192 + 0.0im
     1.715370323453429719112414662767260662417897277748031343163834426881989020798096 + 0.0im
     2.284629676546570280887585337232739337582102722251968656836165573118010979201559 + 0.0im
     2.830830026003772851058548298459246407048009820929073624852446664756879561889998 + 0.0im
     3.309721467890570128113850144932587106367582398673856855212591209011790179898892 + 0.0im
     3.682507065662362337723623297838735435026584996841075797285300235363302213803357 + 0.0im
     3.918985947228994779780736114132655398124909696844323910088990117990380616592766 + 0.0im
julia>

sverek · 2019-01-19T09:54:34Z

Tack! :)

jarlebring · 2019-01-20T20:28:09Z

Any thoughts on adding tutorials in the online documentation? I think this workaround would be suitable there rather than doing some (dirty?) fix in the code.

sverek · 2019-01-20T21:20:35Z

I really appreciate this example, very illuminating.

fgerick · 2019-10-31T10:59:49Z

What needs to be done to get this running for arbitrary number types? I'm very interested to have this working (I don't dare try to compile PETSc and SLEPc in 128 bit versions and somehow interface to it...).

haampie · 2019-10-31T13:00:26Z

JuliaLang/julia#29634 landed in Julia in the meantime, so what has to be changed is

Remove
https://github.com/haampie/ArnoldiMethod.jl/blob/da85958072ec54c64c7dca8362eb42e572252f62/src/expansion.jl#L3

And replace https://github.com/haampie/ArnoldiMethod.jl/blob/da85958072ec54c64c7dca8362eb42e572252f62/src/expansion.jl#L35 by the corresponding mul! function

fgerick · 2019-11-11T10:53:14Z

Is the mul! change going to be in 1.3 when it is released? Anyways, for the moment I have just defined the gemv! for generic number types, as suggested by @jarlebring . There is also a problem in the lu function in schursort.jl related to StaticArrays of non-isbitstype numbers:

using ArnoldiMethod, LinearAlgebra
import LinearAlgebra.BLAS.gemv!
function gemv!(C,alpha,A,X,beta,Y)
    if (C=='N')
       Y.=alpha*A*X+beta*Y
    elseif (C=='T')
       Y.=alpha*transpose(A)*X+beta*Y
    elseif (C=='C')
       Y.=alpha*A'*X+beta*Y
    else
       error("Unknown transpose character")
    end
end

a=big.(randn(10,10))
partialschur(a)

leads to:

ERROR: setindex!() with non-isbitstype eltype is not supported by StaticArrays. Consider using SizedArray.
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] setindex! at /Users/gerickf/.julia/packages/StaticArrays/DBECI/src/MArray.jl:109 [inlined]
 [3] macro expansion at /Users/gerickf/.julia/packages/StaticArrays/DBECI/src/indexing.jl:65 [inlined]
 [4] _setindex!_scalar at /Users/gerickf/.julia/packages/StaticArrays/DBECI/src/indexing.jl:46 [inlined]
 [5] setindex! at /Users/gerickf/.julia/packages/StaticArrays/DBECI/src/indexing.jl:42 [inlined]
 [6] lu(::StaticArrays.SArray{Tuple{2,2},BigFloat,2,4}, ::Type{ArnoldiMethod.CompletePivoting}) at /Users/gerickf/.julia/packages/ArnoldiMethod/5fDBS/src/schursort.jl:98
 [7] sylv at /Users/gerickf/.julia/packages/ArnoldiMethod/5fDBS/src/schursort.jl:184 [inlined]
 [8] swap12!(::Array{BigFloat,2}, ::Int64, ::Array{BigFloat,2}) at /Users/gerickf/.julia/packages/ArnoldiMethod/5fDBS/src/schursort.jl:408
 [9] swap! at /Users/gerickf/.julia/packages/ArnoldiMethod/5fDBS/src/schursort.jl:479 [inlined]
 [10] rotate_right!(::Array{BigFloat,2}, ::Int64, ::Int64, ::Array{BigFloat,2}) at /Users/gerickf/.julia/packages/ArnoldiMethod/5fDBS/src/schursort.jl:27
 [11] partition_schur_three_way!(::Array{BigFloat,2}, ::Array{BigFloat,2}, ::Array{Int64,1}) at /Users/gerickf/.julia/packages/ArnoldiMethod/5fDBS/src/run.jl:348
 [12] _partialschur(::Array{BigFloat,2}, ::Type{BigFloat}, ::Int64, ::Int64, ::Int64, ::BigFloat, ::Int64, ::LM) at /Users/gerickf/.julia/packages/ArnoldiMethod/5fDBS/src/run.jl:254
 [13] #partialschur#1 at /Users/gerickf/.julia/packages/ArnoldiMethod/5fDBS/src/run.jl:103 [inlined]
 [14] partialschur(::Array{BigFloat,2}) at /Users/gerickf/.julia/packages/ArnoldiMethod/5fDBS/src/run.jl:101
 [15] top-level scope at none:0

Does changing to SizedArray come with a big caveat?

fgerick mentioned this issue Feb 17, 2021

BigFloat compatibility #105

Merged

haampie closed this as completed in #105 Feb 18, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

BigFloat matrix in partialschur() does not work #94

BigFloat matrix in partialschur() does not work #94

sverek commented Jan 11, 2019

haampie commented Jan 13, 2019

andreasnoack commented Jan 14, 2019

jarlebring commented Jan 19, 2019 •

edited

Loading

sverek commented Jan 19, 2019

jarlebring commented Jan 20, 2019

sverek commented Jan 20, 2019

fgerick commented Oct 31, 2019

haampie commented Oct 31, 2019

fgerick commented Nov 11, 2019

BigFloat matrix in partialschur() does not work #94

BigFloat matrix in partialschur() does not work #94

Comments

sverek commented Jan 11, 2019

haampie commented Jan 13, 2019

andreasnoack commented Jan 14, 2019

jarlebring commented Jan 19, 2019 • edited Loading

sverek commented Jan 19, 2019

jarlebring commented Jan 20, 2019

sverek commented Jan 20, 2019

fgerick commented Oct 31, 2019

haampie commented Oct 31, 2019

fgerick commented Nov 11, 2019

jarlebring commented Jan 19, 2019 •

edited

Loading