-
A small project to do a pure port of ARPACK Symmetric solver to Julia - Part 5 - Understanding Ref
Here, we pick up with debugging the call to
dsconvwithin the ARPACK library itself from the last post. This is when the repo was at this commitWhen we run the tests against the Arpack wrapped function directly, it errors.
using Pkg; Pkg.test("ArpackInJulia"; test_args=["arpackjll"])In particular, the following code always returns “-1”
import Arpack_jll, LinearAlgebra function arpack_dsconv(n::Int, ritz::Vector{Float64}, bounds::Vector{Float64}, tol::Float64) nconv::LinearAlgebra.BlasInt = -1 ccall((:dsconv_, Arpack_jll.libarpack), Cvoid, (Ref{LinearAlgebra.BlasInt}, Ptr{Float64}, Ptr{Float64}, Ref{Float64}, Ref{LinearAlgebra.BlasInt}), n, ritz, bounds, tol, nconv) return nconv endUgh. This is tough to debug. It seems like nconv always isn’t getting updated. One of my common refrains is that problems occur at boundaries between code. Interfaces are very hard to get right because slight differences in models can manifest as bugs or just outright incorrect answers. This is exactly what we are seeing here.
Now, I could just jump to the solution, which is blindingly obvious in retrospect. On the other hand, there’s value in walking through the process I used. (Or so I tell myself.)
Given that the code above doesn’t work, what were my working hypotheses about the causes?
- So maybe I don’t understand how the fortran calling stuff works and I shouldn’t have expected this to get updated.
- Maybe Julia is doing something I don’t understand?
- There is some odd interaction for this particular piece of code?
- The Fortran / Arpack code doesn’t do what I think it does.
Well, the last one is the the most important as everything follows from my assumption the fortran code works. So that seems like a good place to start. (In truth, I did futz around with some of the others, but then quickly settled on this one.)
Step 1. Do something you know or strongly suspect will work.
Arpack is designed to work with C and C++. There are just fewer places for bugs to hide there. So I wanted to make sure I could call this from C and get the response I expected.
This gave me the following C program.
#include <inttypes.h> typedef int64_t blasint; void dsconv_(blasint* n, double* ritz, double* bounds, double* tol, blasint *nconv); #include <printf.h> int main(int argc, char** argv) { double ritz[] = {1.0, 1.0, 1.0}; double bounds[] = {0.0, 0.0, 0.0}; double tol = 1e-8; blasint n = 3; blasint nconv = 4; dsconv_(&n, ritz, bounds, &tol, &nconv); printf("nconv = %i\n", (int)nconv); return 0; }Which I saved into a file:
test_dsconv.c. This should print 3 if I’ve understood this correctly. (And the other Julia code is correct.)I wanted to compile this against the Julia libraries. Here’s how I got it working. Get the path to libarpack from Julia
@show Arpack_jll.libarpack_path /Users/dgleich/.julia/artifacts/096eefeb84e5b1a463aefac5a561d6109f6e9467/lib/libarpack.2.0.0.dylib(Ugh, security breach… yes, my local username is dgleich. Entirely surprising!) After much fussing with library paths, etc. here’s how to get this all compiled.
gcc test_dsconv.c \ /Users/dgleich/.julia/artifacts/096eefeb84e5b1a463aefac5a561d6109f6e9467/lib/libarpack.2.0.0.dylib \ -Wl,-rpath /Users/dgleich/.julia/artifacts/096eefeb84e5b1a463aefac5a561d6109f6e9467/lib/ \ -Wl,-rpath /Applications/Julia-1.6.app/Contents/Resources/julia/lib/juliaLet’s pick this apart.
gcc test_dsconv.c # above is the call to compile this piece of code /Users/dgleich/.julia/artifacts/096eefeb84e5b1a463aefac5a561d6109f6e9467/lib/libarpack.2.0.0.dylib # above is telling it we want to use this library. -Wl,-rpath /Users/dgleich/.julia/artifacts/096eefeb84e5b1a463aefac5a561d6109f6e9467/lib/ # above is telling gcc to past an option to the linker command. Specifically # rpath means where to get runtime libraries ("runtime path") # this can also be set with LD_LIBRARY_PATH, but that's... controversialLet’s try without the last link to Julia’s libraries to see what happens!
dgleich@circulant test_dsconv % gcc test_dsconv.c \ /Users/dgleich/.julia/artifacts/096eefeb84e5b1a463aefac5a561d6109f6e9467/lib/libarpack.2.0.0.dylib \ -Wl,-rpath /Users/dgleich/.julia/artifacts/096eefeb84e5b1a463aefac5a561d6109f6e9467/lib/ \ dgleich@circulant test_dsconv % ./a.out dyld: Library not loaded: @rpath/libopenblas64_.dylib Referenced from: /Users/dgleich/.julia/artifacts/096eefeb84e5b1a463aefac5a561d6109f6e9467/lib/libarpack.2.0.0.dylib Reason: image not found zsh: abort ./a.outWhich is because it can’t find libopenblas, referenced by libarpack. So we need it tell it about all the other common Julia libraries. That’s what the final
-Wl,-rpath /Applications/Julia-1.6.app/Contents/Resources/julia/lib/julia # tell it about Julia's libraries.And now it works great.
dgleich@circulant test_dsconv % gcc test_dsconv.c \ /Users/dgleich/.julia/artifacts/096eefeb84e5b1a463aefac5a561d6109f6e9467/lib/libarpack.2.0.0.dylib \ -Wl,-rpath /Users/dgleich/.julia/artifacts/096eefeb84e5b1a463aefac5a561d6109f6e9467/lib/ \ -Wl,-rpath /Applications/Julia-1.6.app/Contents/Resources/julia/lib/julia dgleich@circulant test_dsconv % ./a.out nconv = 3Okay, so my understanding of the Arpack function was correct. Ugh. That means I need to learn something about Julia.
Step 2: Better understanding Julia.
Okay, let’s go back to Julia because what we are doing ought to work. Here’s where @code_native and some hints about assembly code come in handy.
julia> @code_native arpack_dsconv(3, ones(3), zeros(3), 1e-8) .section __TEXT,__text,regular,pure_instructions ; ┌ @ none:1 within `arpack_dsconv' subq $40, %rsp ; │ @ none:4 within `arpack_dsconv' ; │┌ @ essentials.jl:396 within `cconvert' ; ││┌ @ refpointer.jl:104 within `convert' ; │││┌ @ refvalue.jl:8 within `RefValue' movq %rdi, 32(%rsp) vmovsd %xmm0, 16(%rsp) movq $-1, (%rsp) ; │└└└ ; │┌ @ pointer.jl:65 within `unsafe_convert' movq (%rsi), %rsi movq (%rdx), %rdx leaq 32(%rsp), %rdi leaq 16(%rsp), %rcx movq %rsp, %r8 movabsq $dsconv_, %rax ; │└ callq *%rax ; │ @ none:11 within `arpack_dsconv' movq $-1, %rax addq $40, %rsp retq nopw %cs:(%rax,%rax) ; └If you look at this code long enough and google around about how
amd64assembly works (orx86_64if you prefer, although it’s really amd64…) then the return value of the function is in registerrax. This codemovq $-1, %raxsets the return register to -1 right before returning (
retq)If you want to learn more about calling conventions and register conventions on
amd64, see this powerpoint from cs217 at PrincetonOkay, so what’s happening is that we are just always returning
-1independently of the function call. This means that the Julia compiler is treating the value ofnconvas a constant and optimizing the code.Another debugging tenant is try things you think will likely work. Here’s something I strongly suspect will work. And just for fun, I tried the following code.
import Arpack_jll, LinearAlgebra function arpack_dsconv(n::Int, ritz::Vector{Float64}, bounds::Vector{Float64}, tol::Float64) nconv::Vector{LinearAlgebra.BlasInt} = [-1] ccall((:dsconv_, Arpack_jll.libarpack), Cvoid, (Ref{LinearAlgebra.BlasInt}, Ptr{Float64}, Ptr{Float64}, Ref{Float64}, Ref{LinearAlgebra.BlasInt}), n, ritz, bounds, tol, nconv) return nconv[1] endWhy do I think this will work? Because I’m passing an array and then explicitly looking at an entry. This is what I want to do, just without the array allocation that creating the vector causes. (Julia’s memory allocater is good, but not that good.)
julia> arpack_dsconv(3, ones(3), zeros(3), 1e-8) 1-element Vector{Int64}: 3Hmm… so I need to learn more about Ref to understand why the other call didn’t work.
This line stuck out to me.
In Julia, Ref objects are dereferenced (loaded or stored) with [].
Anyway, at this point, I realized my mental model was wrong and got to …
The fix
Now, here’s the fix. The issue was the
ccallwill auto wrap an intrinsic type (i.e. int, float, etc.) with a Ref if you tell it the ref type. But this doesn’t refer to the original value, it makes a new copy of it and refers to that. So here’s working code.import Arpack_jll, LinearAlgebra function arpack_dsconv(n::Int, ritz::Vector{Float64}, bounds::Vector{Float64}, tol::Float64) nconv = Ref{LinearAlgebra.BlasInt}(-1) ccall((:dsconv_, Arpack_jll.libarpack), Cvoid, (Ref{LinearAlgebra.BlasInt}, Ptr{Float64}, Ptr{Float64}, Ref{Float64}, Ref{LinearAlgebra.BlasInt}), n, ritz, bounds, tol, nconv) return nconv[] endThe error in my mental model
As I said above, bugs exist on boundaries. Here’s the origin of this bug and my incorrect mental model of the code.
My mental model based on C/C++
nconv::Int = -1 ccall((:func, libfunc), Cvoid, (Ref{Int}), nconv) # <-> func(&nconv) in CBut that’s not right! Here’s what the code actually does.
nconv::Int = -1 ccall((:func, libfunc), Cvoid, (Ref{Int}), nconv) # <-> int newint = nconv; func(&int(newint)) in CWhich of course explains why nconv isn’t updated. It wasn’t even passed in!
Suggestion
- Add a note to
Refdocumentation on this - Add a note to
ccalldocumentation on this
Is it worth it for one person with an incorrect mental model? I think yes because the
ccallarea is likely to be a source of similar bugs and it’s needed to document the assumptions! This is already partially done in this paragraphWhen passed as a
ccallargument (either as aPtrorReftype), aRefobject will be converted to a native pointer to the data it references. For mostT, or when converted to aPtr{Cvoid}, this is a pointer to the object data. WhenTis anisbitstype, this value may be safely mutated, otherwise mutation is strictly undefined behavior. -
A small project to do a pure port of ARPACK Symmetric solver to Julia - Part 4 - Checking against libarpack
See Overview for a high-level picture, also see ArpackInJulia.jl on Github
One of the things I wanted to do was check against the true ARPACK functions. The idea is that we can swap out true ARPACK routines for our routines easily. So I wanted to be able to call the true ARPACK routines individually myself.
This is remarkably easy and how Arpack.jl works with the functions
aupdandeupdin ARPACK to implement eigs. So this isn’t surprisingly. It’s still fun to see that ease of use translate into new codes.The idea is that
Arpack.jlusesArpack_jllto access the compiled FORTRAN code as a dynamic library. All the binary compiling, etc. is handled byArpack_jllandArpack.jlcan just use that. So we are going toimport Arpack_jllso we can get access to
@show Arpack_jll.libarpackwhich gives the path to the ARPACK library in a system independent fashion. (At some point, this will work on the Apple M1 too, but not yet.) [Probably will be about the same time as I finish posting these notes though and getting all of the routines ported.]
Once we have that, we can simply setup a ccall to get what we need! See the blas.jl functions for many examples of how to do this. Also see Arpack.jl for even more examples that are specific to ARPACK.
import Arpack_jll, LinearAlgebra function arpack_dsconv(n::Int, ritz::Vector{Float64}, bounds::Vector{Float64}, tol::Float64) nconv::Int = 0 ccall((:dsconv_, Arpack_jll.libarpack), Cvoid, (Ref{LinearAlgebra.BlasInt}, Ptr{Float64}, Ptr{Float64}, Ref{Float64}, Ref{LinearAlgebra.BlasInt}), n, ritz, bounds, tol, nconv) return nconv endOf course, I only want to depend on
Arpack_jllin the test routines. (Actually, I’d really only like to depend on it in a subset of the test routines. But I don’t see how to do that right now and it isn’t relevant.)Put another way, I want to enable the tests to run without
Arpack_jllworking on a particular platform. And where we don’t run the tests that need it.To control which tests are run, we can check
ARGSas described in this pull request https://github.com/JuliaLang/Pkg.jl/pull/1226 where this functionality was implemented in Julia.if "arpackjll" in ARGS include("arpackjll.jl") # get the helpful routines to call stuff in "arpackjll" @testset "arpackjll" begin # tests that check against libarpack directly end endSo I stuck the call to dsconv into
arpackjll.jlimport Arpack_jll, LinearAlgebra function arpack_dsconv(n::Int, ritz::Vector{Float64}, bounds::Vector{Float64}, tol::Float64) nconv::LinearAlgebra.BlasInt = -1 ccall((:dsconv_, Arpack_jll.libarpack), Cvoid, (Ref{LinearAlgebra.BlasInt}, Ptr{Float64}, Ptr{Float64}, Ref{Float64}, Ref{LinearAlgebra.BlasInt}), n, ritz, bounds, tol, nconv) return nconv endThis does require me to add
Arpack_jllandLinearAlgebrato the Test package itself. We can do this by activating the test package and then adding them.# change to the test subdirectory ] activate . add LinearAlgebra add Arpack_jllAnd our full piece of new test code
if "arpackjll" in ARGS include("arpackjll.jl") # get the helpful routines to call stuff in "arpackjll" @testset "arpackjll" begin @testset "dsconv" begin soln = arpack_dsconv(10, ones(10), zeros(10), 1e-8) [@test](https://micro.blog/test) ArpackInJulia.dsconv(10, ones(10), zeros(10), 1e-8)==soln soln = arpack_dsconv(10, zeros(10), ones(10), 1e-8) [@test](https://micro.blog/test) ArpackInJulia.dsconv(10, zeros(10), ones(10), 1e-8)==soln end end endThis flips ones and zeros, which should show converged in one of the cases.
And what’s weird, it fails. I always get back arpack_dsconv==-1. Hmm… will debug in the next post!
-
A funny thing in @JuliaLanguage. I often illustrate that finite differences in floating point are inaccurate with small changes. But this Julia code works and it’s fun to understand why!
xx=range(0,1,length=100) F(x) = 5.0.*x .+ 1/3 h = eps(1.0) z = (F(xx.+h)-F(xx))./h
-
A small project to do a pure port of ARPACK Symmetric solver to Julia - Part 3 - First routine
See Overview for a high-level picture, also see ArpackInJulia.jl on Github
We are going to start our ARPACK conversion with about the simplest routine:
dsconv. This routine for a double-precision, symmetric eigenvalue problem just checks if the eigenvalue estimates have converged at all.c----------------------------------------------------------------------- c\BeginDoc c c\Name: dsconv c c\Description: c Convergence testing for the symmetric Arnoldi eigenvalue routine. c c\Usage: c call dsconv c ( N, RITZ, BOUNDS, TOL, NCONV ) c c\Arguments c N Integer. (INPUT) c Number of Ritz values to check for convergence. c c RITZ Double precision array of length N. (INPUT) c The Ritz values to be checked for convergence. c c BOUNDS Double precision array of length N. (INPUT) c Ritz estimates associated with the Ritz values in RITZ. c c TOL Double precision scalar. (INPUT) c Desired relative accuracy for a Ritz value to be considered c "converged". c c NCONV Integer scalar. (OUTPUT) c Number of "converged" Ritz values. c c\EndDoc c c----------------------------------------------------------------------- c c\BeginLib c c\Routines called: c arscnd ARPACK utility routine for timing. c dlamch LAPACK routine that determines machine constants. c c\Author c Danny Sorensen Phuong Vu c Richard Lehoucq CRPC / Rice University c Dept. of Computational & Houston, Texas c Applied Mathematics c Rice University c Houston, Texas c c\SCCS Information: @(#) c FILE: sconv.F SID: 2.4 DATE OF SID: 4/19/96 RELEASE: 2 c c\Remarks c 1. Starting with version 2.4, this routine no longer uses the c Parlett strategy using the gap conditions. c c\EndLib c c----------------------------------------------------------------------- c subroutine dsconv (n, ritz, bounds, tol, nconv) c c %----------------------------------------------------% c | Include files for debugging and timing information | c %----------------------------------------------------% c include 'debug.h' include 'stat.h' c c %------------------% c | Scalar Arguments | c %------------------% c integer n, nconv Double precision & tol c c %-----------------% c | Array Arguments | c %-----------------% c Double precision & ritz(n), bounds(n) c c %---------------% c | Local Scalars | c %---------------% c integer i Double precision & temp, eps23 c c %-------------------% c | External routines | c %-------------------% c Double precision & dlamch external dlamch c %---------------------% c | Intrinsic Functions | c %---------------------% c intrinsic abs c c %-----------------------% c | Executable Statements | c %-----------------------% c call arscnd (t0) c eps23 = dlamch('Epsilon-Machine') eps23 = eps23**(2.0D+0 / 3.0D+0) c nconv = 0 do 10 i = 1, n c c %-----------------------------------------------------% c | The i-th Ritz value is considered "converged" | c | when: bounds(i) .le. TOL*max(eps23, abs(ritz(i))) | c %-----------------------------------------------------% c temp = max( eps23, abs(ritz(i)) ) if ( bounds(i) .le. tol*temp ) then nconv = nconv + 1 end if c 10 continue c call arscnd (t1) tsconv = tsconv + (t1 - t0) c return c c %---------------% c | End of dsconv | c %---------------% c endIn Julia, this routine is basically a two-liner (could likely do in one) that computes:
eps23 = eps(1.0)^(2/3) return sum( (ritzi, boundsi)->boundsi <= max(eps23, abs(ritzi))*tol ? 1 : 0, zip(ritz, bounds))(Full disclosure, I didn’t test this, so there could be a subtle error, but I hope the idea is clear that this COULD be done very easily.)
On the other hand, the idea was to port this as directly as possible, so let’s write out the loops, etc.
function dsconv( n::Int, ritz::Vector{Float64}, bounds::Vector{Float64}, tol::Float64 ) eps23::Float64 = eps(1.0)^(2.0/3.0) nconv::Int = 0 if n > length(ritz) throw(ArgumentError("range 1:n=$n out of bounds for ritz, of length $(length(ritz))")) end if n > length(bounds) throw(ArgumentError("range 1:n=$n out of bounds for bounds, of length $(length(bounds))")) end @inbounds for i=1:n tmp = max( eps23, abs(ritz[i]) ) if ( bounds[i] <= tol*tmp ) nconv += 1 end end return nconv endHere, we have borrowed some of the out of bounds reporting from blas.jl in order to try and be somewhat standard.
Of course, the issue is that the ARPACK code also includes some cool timing information
c c %----------------------------------------------------% c | Include files for debugging and timing information | c %----------------------------------------------------% c include 'debug.h' include 'stat.h' call arscnd (t0) c ... many lines ... call arscnd (t1) tsconv = tsconv + (t1 - t0)The function
arscndis the ARPACK second function, which gives the current time in seconds. So this simply tracks the seconds used in this function.To store the information, this code uses the Fortran common block, which looks like a set of global variables, and update them based on the time used in this routine. In the port, I’d like to have this functionality as well. However while the the FORTRAN common block could be replaced with a set of global variables, I wanted to use a more modern construct enabled by Julia’s macros that will give thread-safety and also compiler-optimized code removal when not needed.
Here’s a better design. Take in a keyword argument `timing’ that is a julia structure.
Base.@kwdef struct ArpackStats tconv::Float64 = 0.0 # the time sp # ... there will be more! endWe use Base.@kwdef to auto-initialize all the times to zero so we don’t have to write a construction and can just call ArpackStats(). (Check out help on Base.@kwdef for more!)
Then we can simply use the following construct to update time:
function dsconv( n::Int, ritz::Vector{Float64}, bounds::Vector{Float64}, tol::Float64; stats=nothing ) t0 = timing != nothing ? time() : 0.0 ## NEW eps23::Float64 = eps(1.0)^(2.0/3.0) nconv::Int = 0 if n > length(ritz) throw(ArgumentError("range 1:n=$n out of bounds for ritz, of length $(length(ritz))")) end if n > length(bounds) throw(ArgumentError("range 1:n=$n out of bounds for bounds, of length $(length(bounds))")) end @inbounds for i=1:n tmp = max( eps23, abs(ritz[i]) ) if ( bounds[i] <= tol*tmp ) nconv += 1 end end if stats != nothing; stats.tconv += time() - t0; end ## NEW return nconv endOf course, these if statements get a bit tedious. So let’s use Julia’s macros to handle this more seemlessly to get what I want to write.
""" Convergence testing for the symmetric Arnoldi eigenvalue routine. In Julia, this can be replaced with sum( ) Julia port of ARPACK dsconv routine. Usage: call dsconv ( N, RITZ, BOUNDS, TOL; debug, stats ) Arguments N Integer. (INPUT) Number of Ritz values to check for convergence. RITZ Double precision array of length N. (INPUT) The Ritz values to be checked for convergence. BOUNDS Double precision array of length N. (INPUT) Ritz estimates associated with the Ritz values in RITZ. TOL Double precision scalar. (INPUT) Desired relative accuracy for a Ritz value to be considered "converged". Optional values (they default to nothing) debug (INPUT/OUTPUT) nothing or an ArpackStats struct stats (INPUT/OUTPUT) nothing or an ArpackStats struct These are used to collect debugging and performance stats and are mostly internal. Return value NCONV Integer scalar. (RETURN) Number of "converged" Ritz values. """ function dsconv( n::Int, ritz::Vector{Float64}, bounds::Vector{Float64}, tol::Float64; stats=nothing, debug=nothing ) #= c\Author c Danny Sorensen Phuong Vu c Richard Lehoucq CRPC / Rice University c Dept. of Computational & Houston, Texas c Applied Mathematics c Rice University c Houston, Texas Julia port by David F. Gleich, Purdue University =# #= c | Executable Statements | =# t0 = @jl_arpack_time if n > length(ritz) throw(ArgumentError("range 1:n=$n out of bounds for ritz, of length $(length(ritz))")) end if n > length(bounds) throw(ArgumentError("range 1:n=$n out of bounds for bounds, of length $(length(bounds))")) end eps23::Float64 = eps(1.0)^(2.0/3.0) nconv::Int = 0 @inbounds for i=1:n #= c | The i-th Ritz value is considered "converged" | c | when: bounds(i) .le. TOL*max(eps23, abs(ritz(i))) | =# tmp = max( eps23, abs(ritz[i]) ) if ( bounds[i] <= tol*tmp ) nconv += 1 end end # call arscnd (t1) # tsconv = tsconv + (t1 - t0) @jl_update_time(tconv, t0) return nconv endHere, we can use the following macros to accomplish this.
macro jl_arpack_time() return esc(:( stats == nothing ? zero(ArpackTime) : time() )) end macro jl_update_time(field, t0) return esc( quote if stats != nothing stats.$field += time() - $t0 end end ) endAnd thus, we have our actually strategy. Note that some of the little code snippets here were taken from development versions and some of the variable names may have changed. (e.g.
tconvvs.tsconv) because there is no reason to track symmetric vs. nonsymmetric time separately. (I mean, I guess we could, but that seems unlike Julia language stuff.) Could be revisited in the future!At this point, we also put a bunch of stuff in Github to keep the project coming. See the ArpackInJulia.jl Github repo.
-
A small project to do a pure port of ARPACK Symmetric solver to Julia - Part 2 - Overview
See Overview for a high-level picture.
Here is a quick call graph of the various ARPACK functions for the symmetric eigenvalue/eigenvector problem
dsaupdd= double precisions= symmetricaupd= Arnoldi Update (since we’ll call this in a loop to update the internal Arnoldi information.)
Note to the pedantic. Yes, I know the Arnoldi of a symmetric matrix is usually called the Lanczos process. So I should be saying Lanczos everywhere here. Except, maybe I really shouldn’t be saying either and we could use the phrase: iteratively orthogonalized power basis instead of Arnoldi or three-term recurrent power basis instead of Lanczos? Either way, since the package is called ARPACK and I’m writing about that, I’ll stick with Arnoldi here.
Back to the call graph.
dsaupd calls dsaup2 calls dgetv0 calls printing, reverse communication dsaitr calls printing, blas, lapack dsapps calls printing, blas, lapack dsconv dseigt dsgets dsortr ivout - write out an integer vector arscnd - internal wrapper avoiding LAPACK second function dvout - write out a double vectorthen we also need to handle getting the eigenvector information out of the computation.
dseupd calls dsesrt dsortr - same as aboveThis suggests a plan.
Port the low-level/internal routines first. Then port the high-level routines.
Wait what you say? Why not do it the opposite way? Well, I started doing it that way and realized it’s just harder ;). So we are going to switch and do the easy stuff first!
Just as a quick overview of what the code actually does.
dsaupd - high level interface dsaup2 - medium level interface, manages the overall process dgetv0 - gets the starting vector, or returns control to user to get it. dsaitr - haven't looked yet... dsapps - haven't looked yet... dsconv - check what has converged (easiest routine!) dseigt - haven't looked yet... dsgets - get shifts based on eigenvalue estimates dsortr - sorts a vector and/or applies the sort to an aux vector tooAlso, let me note that the ARPACK documentation is stellar. Stellar. So a few clicks and I could easily figure out what those other functions are doing.
-
A small project to do a pure port of ARPACK Symmetric solver to Julia - Part 1 - Goals
I’m working on a small “me-time” project to port the ARPACK symmetric solver to Julia.
My goal is to port the double-precision ARPACK for symmetric matrices in Julia. Including all ARPACK stuff. So this should give “exactly” what ARPACK does but be a pure Julia implementation. (Where exactly is … it should be executing roughly the same sequence of floating point operations and can differ on levels that would be expected for different compilers compiling the same code.)
- not a goal to “Julia-ize” the package; I want to keep as close to the FORTRAN as possible so that I might be able to replace calls to Julia’s Arpack.saupd / Arpack.seupd (which call the Fortran library) with this code.
- small internal function changes are okay, e.g. ARPACK has various debugging and timing stuff that would need to be done differently in Julia.
- small simplifications, e.g. if a function computes a single Int, we can rewrite that to return the Int rather than writing it into an array like in FORTRAN.
- Why? Why not implement my own ImplicitRestart/Eigensolver/Etc.? Simple: I trust ARPACK. Also, I want to understand exactly what the symmetric ARPACK solver is doing.
- Why not use a Fortran->Julia compiler? Well, I could. But I could also do this and learn intimate details of how it works to help out in teaching :)
- I want to document some of the cool stuff in ARPACK!
For context, I spent about a semester studying these solvers with Gene Golub to see if there were obvious ideas to make computing the Fiedler vectors faster, but beyond the high level algorithmics and lots of experiments, nothing came of it. This project is different in that the end goal is easy – port the code, with a focus on explaining in detail what each step does. (hah, how about that name-dropping?)
Since it’s good to have a driving use-case, the idea is to make the following MatrixNetworks.jl code work without the ARPACK dependency.
""" `_symeigs_smallest_arpack` --- Compute eigenvalues and vectors using direct calls to the ARPACK wrappers to get type-stability. This function works for symmetric matrices. This function works on Float32 and Float64-valued sparse matrices. It returns the smallest set of eigenvalues and vectors. It only works on matrices with more than 21 rows and columns. (Use a dense eigensolver for smaller problems.) Functions --------- - `(evals,evecs) = _symeigs_smallest_arpack(A::SparseMatrixCSC{V,Int}, nev::Int,tol::V,maxiter::Int, v0::Vector{V})` Inputs ------ - `A`: the sparse matrix, must be symmetric - `nev`: the number of eigenvectors requested - `tol`: the relative tolerance of the eigenvalue computation - `maxiter`: the maximum number of restarts - `v0`: the initial vector for the Lanczos process Example ------- This is an internal function. """ function _symeigs_smallest_arpack( A::SparseMatrixCSC{V,Int},nev::Int,tol::V,maxiter::Int, v0::Vector{V}) where V n::Int = checksquare(A) # get the size @assert n >= 21 # setup options mode = 1 sym = true iscmplx = false bmat = String("I") ncv = min(max(2*nev,20),n-1) whichstr = String("SA") ritzvec = true sigma = 0. #TOL = Array{V}(undef,1) #TOL[1] = tol TOL = Ref(tol) lworkl = ncv*(ncv + 8) v = Array{V}(undef, n, ncv) workd = Array{V}(undef, 3*n) workl = Array{V}(undef, lworkl) resid = Array{V}(undef, n) resid[:] = v0[:] info = zeros(BlasInt, 1) info[1] = 1 iparam = zeros(BlasInt, 11) ipntr = zeros(BlasInt, 11) ido = zeros(BlasInt, 1) iparam[1] = BlasInt(1) iparam[3] = BlasInt(maxiter) iparam[7] = BlasInt(mode) # this is a helpful indexing vector zernm1 = 0:(n-1) while true # This is the reverse communication step that ARPACK does # we need to extract the desired vector and multiply it by A # unless the code says to stop Arpack.saupd( ido, bmat, n, whichstr, nev, TOL, resid, ncv, v, n, iparam, ipntr, workd, workl, lworkl, info) load_idx = ipntr[1] .+ zernm1 store_idx = ipntr[2] .+ zernm1 x = workd[load_idx] if ido[1] == 1 workd[store_idx] = A*x elseif ido[1] == 99 break else error("unexpected ARPACK behavior") end end # Once the process terminates, we need to extract the # eigenvectors. # calls to eupd howmny = String("A") select = Array{BlasInt}(undef, ncv) d = Array{V}(undef, nev) sigmar = ones(V,1)*sigma ldv = n Arpack.seupd(ritzvec, howmny, select, d, v, ldv, sigmar, bmat, n, whichstr, nev, TOL, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info) if info[1] != 0 error("unexpected ARPACK exception") end # Now we want to return them in sorted order (smallest first) p = sortperm(d) d = d[p] vectors = v[1:n,p] return (d,vectors,iparam[5]) end -
Anyone else find the documentation in old Fortran codes is amazingly beautiful? e.g. look at the documentation in ARPACK. github.com/opencolla… Wonderfully annotated code.
-
Julia (@JuliaComputing) and VSCode on Native M1.
Install M1 Native VSCode.
Install Xcode & Command Line tools.
Clone Julia master.
make
./julia
ENV[“JULIA_GR_PROVIDER”] = “GR”
]add Plots
using Plots
plot(randn(10))
Point VSCode at new julia and 🍾🤩🥳
-
My thoughts on getting a COVID vaccine today.
My thoughts on getting a COVID vaccine today. I’m of mixed feelings about this because, while vaccine!, strictly speaking, the eligibility of University faculty (as I understand them) are unclear. I’m also not particularly old, nor high-risk.
Meijer (and Kroger, and others) are vaccinating university faculty as educators. (Which I definitely am!) This fact has been posted on Twitter by others, and announced on a university mailing list. While I’d like to think of myself as a as-scrupulous-as-reasonable rule-follower, this is a murky area.
I never lied to get this. Also, remember – any type of medical/healthcare employer was vaccinating nearly everyone when vaccines were really scarce, even those capable of work at home or not doing patient facing care. Here’s part of my rationalization.
-The US has a somewhat inefficient allocation of vaccine supply to where it is needed. I booked this appointment with <= 24 not-sleeping hours notice. I cannot believe this keeps anyone who ought to be eligible from getting a vaccine. (The pharmacy still had same day availability when I got my vaccine.)
-I have the capacity/flexibility/ability to exchange my time/energy to get a vaccine when/where others can’t. Should I not use this to help address the inefficient allocation above?
-I have a responsibility to help keep the daycare safe that we participate in, I struggled with how I would feel if I were to get COVID at this point (knowing that I passed up a vaccine opportunity) and pass it onto anyone at the daycare.
I’m posting this note for a few reasons.
We need to remind people that vaccines help everyone! There is some social signaling and sociology Even if you don’t need it, it will help others by curtailing possible spread, and others may need some encouragement to get it.
By trying to understand the appointment landscape, I helped a number of people get the vaccine weeks earlier than they were originally scheduled.
If you, for any reason, want to get the vaccine. Realize that there is capacity in the system in some places.
And use this information to help others too.
If people want to know my thoughts on the current state of getting a vaccine in Indiana, I’m happy to share more from looking at appointment availability for an elementary teacher I work with. My sense is that appointments are reasonably highly available with <= 24 hours notice if you are willing to travel <= 100 miles.
-
Me. Thinking as I’m looking at my shelf. See a book with obscured title… “Lag…” which I read as “Lagrange.” Me. What book about Lagrange again? The rest of the title “and the brain.” Of course, the real title is “Language and the brain” which is … oh so so so appropriate.
-
arXiv vanity works amazingly well with some small warts.
Running our 80-page survey paper on flow-based network cluster improvement through the processor arXiv vanity shows some impressive results. www.arxiv-vanity.com/papers/20… This is not perfect, but pretty amazing considering some of the tables and embedded “aside” paragraphs translate.
The biggest bummer for me is that our macros aren’t translated at all. e.g. www.arxiv-vanity.com/papers/20…
-
Visualizing COVID Cases and Tests - Purdue vs. Vermont
As a fun little exercise with my students, we sought to improve the Purdue COVID tracking dashboard. Our quick take was that this display focused too heavily on tests instead of the more relevant positive cases.
After a few iterations, we wanted a design that had the same bits of information, but presented in a way that made them easier to understand. For instance, we ditched the two-axis design with tests and cases mixed and instead scaled the testing volume down to get a ‘positive rate’ line that acts as a means of evaluating the case counts. (After all, that’s what you’d use the test count data for!)
Here’s the design.
The solid vertical bars are the daily case counts. The thick-black line is the running 7-day average. The internal color of this line gives the positivity rate where 1% is green and 9% is red. The dashed colored lines show the positivity rate given implied by running 7 day count case count. (There are strong day-of-week effects, so the raw data isn’t as meaningful, but 7 day trends are reliable.)
Here, you can see that testing is basically constant through Nov 10 because the constant rate of positivity bars. You can also quickly see that cases are growing recently. (That’s also true everywhere! Maybe we need another plot that compares these to successively larger regions!)
For fun, I tried the same plot for the state of Vermont (I have a relative at UVM!). Vermont has ~623k people. (12x Purdue’s 50k).
In Vermont, we can see a few testing spikes in the Vermont data, and also the recent uptick in cases. (The last few days are actually quite a bit worse in Vermont, but I wanted the same data region.)
I think the Vermont data shows the visualization isn’t quite as successful in general as we might have hoped because testing there seems less consistent than Purdue.
Caveats
- These were meant to illustrate a visualization (hah, but the title implies a comparison). A good visualization facilities comparison, in which case we might ask more refined questions.
- Critically evaluate the data yourself before drawing any conclusions from these pictures (e.g. maybe I used the wrong field from the covid tracking project data for Vermont!) There may be a typo in my Purdue case counts, which we transcribed by hand.
Here’s the code
## Download vermont using CSV data = CSV.read(download("https://covidtracking.com/data/download/vermont-history.csv")) ## align dates ds = data[6:6+101,:] ## ds.totalTestsPeopleViralIncrease ## ds.positiveIncrease using Plots gr() function covidplot(ntests, npositive) A = [ntests npositive] avgs = [] ratios = [] tvol = [] for i = 1:size(A,1) si = max(1,i-6) ei = min(size(A,1),i) push!(avgs,mean(A[si:ei,2])) push!(ratios,sum(A[si:ei,2])/sum(A[si:ei,1])) push!(tvol, sum(A[si:ei,1])/length(si:ei)) end cmap = cgrad([:limegreen, :goldenrod1, :firebrick1], [0.2, 0.5, 0.8]) xx=-size(A,1)+1:0 Plots.bar(xx,A[:,2],alpha=0.4, label="daily cases", linewidth=0, legend=:topleft, #fill_z = 100*A[:,2]./A[:,1], #color=cmap, clims=(1,9), legendfontvalign=:top) for p in [1,3,5,7,9] plot!(xx, (p/100)*tvol, linestyle=:dash, line_z=p, label="", color=cmap, clims=(1,9)) annotate!(1, (p/100)*tvol[end], Plots.text("$p%", 10, get(cmap,p/10), :left)) end plot!() #plot!(1:size(A,1),avgs, color=:black, linewidth=3.5, #label="") plot!(xx,avgs, color=:black, linewidth=3.5, label="avg. cases over 7 days", legendfontvalign=:top) plot!(xx,avgs, line_z = 100*ratios, color=cmap, linewidth=1.5, label="", clims=(1,9), colorbar=false, foreground_color_legend=nothing, background_color_legend=nothing, legendtitle="Colors show positivity rate", legendtitlefontsize=10, legendtitlefontvalign=:bottom, legendtitlefonthalign=:right, size=(500,300), legendfontvalign=:top) plot!(xlabel="days since November 10, 2020") plot!(ylabel="PCR positive SARS-nCOV-2 tests") end covidplot(reverse(ds.totalTestsPeopleViralIncrease), reverse(ds.positiveIncrease)) ylims!(0,121) plot!(legend=:topright) -
Project and Forget Metric Optimization - Explained
In a new video, I do a deep dive with Rishi Sonthalia and Anna Gilbert about their paper on optimization with metric constraints – and their associated code together with my former student Nate Veldt who worked on similar problems with me. Their stuff is super cool and we go into why it’s awesome!
The video - <youtu.be/OBC4FmUui…>
Their paper - <arxiv.org/abs/2005….>
Their code - <github.com/rsonthal/…>
More context – a month or so ago (who can keep track of time – maybe this was September, actually) I sat down with Rishi Sonthalia, my former student Nate Veldt, and Anna Gilbert to understand Rishi and Anna’s latest paper on “Project and Forget Metric Optimization” <arxiv.org/abs/2005….>. For fun, we decided to record this so that I’d be able to look at this again when I wanted to keep working on this. But it turned out so great we figured we’d share it! (After spending too much time editing it, of course…!)
I care about this because it gives a faster way to solve a quadratic or linear program based on a relaxation of the correlation clustering problem. I think it also gives a great place to build future HPC-style work on graphs where there is a problem that takes considerable time and computing effort. The particular bottleneck in the current code is computing an all-pairs-shortest path computation, which is somewhere where there has been a tremendous amount of work that could be used to make this go even faster!
-
Fixing the NYTimes Horrible Exit Poll page.
I have so many problems with this display of data in the attached screenshot from the nytimes https://www.nytimes.com/interactive/2020/11/03/us/elections/exit-polls-president.html. It’s so easy to misread this in very subtle ways. E.g. this makes it seem like most trump voters don’t trust their counts accurately? But that’s not what the data show at all.
The second one (without the pictures) is an improvement I managed to make in excel. This wouldn’t pass muster (yet) if the PhD students I work with tried to include it in a paper, but it’s close enough and I did it in Excel. Mainly to demonstrate it’s possible. (In retrospect, Inkscape would have been easier…)
This shows almost no difference among the groups of voters! Which is, of course, because we are all Americans and largely trust our voting and make up our minds earlier.
-
Designing for Julia is hard
As a longtime Matlab user, someone who enjoys Python syntax, and loves crafting C codes for maximum efficiency, I often struggle with thinking about code interfaces in Julia.
Why? Because the interface possibilities in Julia are extremely powerful. This power can be difficult to deploy, however. Simply put, other languages have conventions or requirements that are essentially constraints on what interfaces ought to look like. These constraints restrict the design space, once they are understood, and simplify design. In contrast, Julia offers more possibilities.
Here are a few relevant thoughts.
Matlab interfaces tended to be holistic functional interfaces that ‘do everything’. Although one might think these are hard to do, packing all the relevant information into a single function call is a useful design constraint that forces or suggest certain choices.
Why isn’t C harder than Julia? The reason that C isn’t harder than Julia is because C interfaces essentially need to fix a small set of types to be reasonable.
Aside C++ templated interfaces are essentially equivalent to many of the design challenges in Julia. Designing these right and Julia interfaces present similar challenges.
Python promotes the class as a general structure to organize interfaces.
Consider what linear programming might look like in these different languages.
Matlab
optval, optarg = linear_program(c,A,b) # solve min c'*x s.t. Ax = b, x >= 0 optval, optarg, dualvars = linear_program(c,A,b) # and get dual variables too optval, optarg, dualvars = linear_program(c,A,b,'method','ip','tol',1e-3)C99 for single line comments
p = alloc_problem(c,A,b) s = alloc_solution(p) solve_lp(p, &s) # s has fields, optval, optarg, dualvars o = interior_point_options() o.tol = 1e-3 solve_lp_interior_point(p, &s, o) // s has fields, optval, optarg, dualvars free(p); free(s); // very important in CPython
p = opt.LinearProgram(c,A,b) M = opt.algs.LinearProgramDefault() soln = M.solve(p) # soln has optval, optarg # get optval, optarg, dualvars from soln M2 = opt.algs.LinearProgramInteriorPoint({'tol':1e-3)} soln = M2.solve(p)Notice that the Matlab interface is probably the easiest and simplest, if far more limited.
There are a number of common challenges that these interfaces have to deal with.
- How to specify the problem
- How to specify a solver
- How to parametrize the solver
- How to get solution pieces.
Matlab’s do everything solver makes this interface easy to write. It also enforces a design funnel that would channel alternatives into this uber-function. (Matlab loves the uber-function, it’s a great and timeless design!) The other two interfaces are largely similar, but for Python, where should the solve function be inside the class for the problem (P.solve(M) or M.solve(P)). (Probably M.solve(P) as then the solver can extract the relevant information it needs from the problem, whereas otherwise the problem has to fix an interface to a method.)
In C, it’s just a function-fest everywhere, and don’t forget to free.
Julia with Convex.jl
x = Variable() problem = minimize(c'*x, [A*x == b, x >= 0]) solve!(problem, SCS.Optimizer) problem.optval, xJulia with JuMP
model = Model(SCS.Optimizer) @variable(model, 0 <= x) @objective(model, Min, c'*x) @constraint(model, A*x == b) optimize!(model)Both of these interfaces exploit Julia’s powerful syntax transformations for all sorts of clever features that are largely impossible in C or Matlab. But back to Matlab’s advantages, I don’t want a super-general optimization tool, I just want to solve a simple LP.
Maybe I’m writing a package where I expect people to play around with various LP solvers, or I want to give them that option. How would I do this? This is where Julia design paralysis sets in. It’s definitely possible, but there are many possible ways to treat these issues and few conventions that would simplify or suggest the design.
So what’s my recommendation or advice? The simplest thing is to treat Julia like C and fix a small set of types that will work with the code. Then try and grow this into more general interfaces over time. This is essentially a restatement of Knuth’s “premature optimization dictum” https://en.wikiquote.org/wiki/Donald_Knuth
The real problem is that programmers have spent far too much time worrying about efficiency in the wrong places and at the wrong times; premature optimization is the root of all evil (or at least most of it) in programming.
Don’t be afraid to hardcode Float64, Int, etc. and don’t feel obliged to be fully type general. (Type general functions are very hard to get right and are best done after a few tries with simpler hardcoded type functions.) This will have the side effect of making fewer calls to the Julia compiler when you call functions.
Avoid trying to use clever Julia tricks in initial designs. Keep it simple and straightforward like C.
Don’t provide flexibility where it isn’t (yet) useful.
Don’t be afraid to make horrible breaking changes in future code revisions. The Julia devs drove some of us nuts with constant breaking syntax changes through 0.2 - 0.7 as they evolved and revisited decisions. Designing without reality is impossible.
Back to Knuth though, with a variant with yet more relevant context.
Programmers waste enormous amounts of time thinking about, or worrying about, the speed of noncritical parts of their programs, and these attempts at efficiency actually have a strong negative impact when debugging and maintenance are considered. We should forget about small efficiencies, say about 97% of the time: premature optimization is the root of all evil. Yet we should not pass up our opportunities in that critical 3%.
The goal of these initial design ideas is to find that 3%. That can only be done once an initial framework is established.
And too close, let’s go back to Knuth.
Science is knowledge which we understand so well that we can teach it to a computer; and if we don’t fully understand something, it is an art to deal with it.
This exemplifies good Julia design, and why it is so hard. Good Julia design flows from a deep understanding of the topic and an appreciation of where similarities and differences lie that is often a contribution in itself.
subscribe via RSS