SparseAlgebra #12

FelixBangerter · 2018-07-19T14:43:21Z

…nclude MatrixModule.jl

coveralls · 2018-07-19T14:58:26Z

Coverage decreased (-16.1%) to 69.308% when pulling e6a0d41 on PicricAcid96:pull-request/e6a0d41a into 69bde0d on vonDonnerstein:master.

vonDonnerstein

Hi Felix, Thank you for your pull request. I've just had a quick read through your changes (see my comments). Looks very promising! I've mostly commented on more general aspects. Since you know best how to write your code, I hope I wasn't too intrusive about any implementational details. Only if you like, I could have a go on those as well. Please add a few unit or integration tests, so that I or anyone else working with your code doesn't unintentionally break your work. Thanks!

vonDonnerstein · 2018-07-19T14:54:22Z

src/SubModules/MatrixModule.jl

@@ -60,7 +66,7 @@ function writePixmap(filename::String,matrix::Matrix)
            ',' => "#ffeebb",
            '.' => "#ffeeee",
            ' ' => "#ffffff",
-	    '0' => "#bbbbbb", # zero by matrix format (blocking)
+		    '0' => "#bbbbbb", # zero by matrix format (blocking)


intentional?

I've added a comment to http://schurkus.com/2015/12/05/quantumlab-design-guidelines/ to clarify indentation:
I'm assuming one tabstop corresponds to 8 spaces (github default)

Wasn't intentional I changed it back.

vonDonnerstein · 2018-07-19T14:54:29Z

src/SubModules/MatrixModule.jl

+	rowptr::Array{Int,1}
+	rowpattern::Array{Tuple{Int64,Int64},1}
+	colpattern::Array{Tuple{Int64,Int64},1}
+    BCSRSpM(val,col,rowptr,rowpattern,colpattern) = new(val,col,rowptr,rowpattern,colpattern)


isn't this just the default constructor? In that case you could leave it out.

vonDonnerstein · 2018-07-19T14:56:27Z

src/SubModules/MatrixModule.jl

+
+
+"""
+    BCSRSpM(val,col,rowptr,rowpattern,colpattern)


The CSC matrix format is called SparseMatrixCSC in julia, so I would think SparseMatrixBCSR would be a consistent name for your type. What do you think?

I agree, I changed it.

vonDonnerstein · 2018-07-19T15:02:59Z

src/SubModules/MatrixModule.jl

+end
+
+"""
+    symBCSRSpM(val,col,rowptr,pattern)


This I would probably call SparseMatrixBCSRSymmetric or something in that vein. Similar to the one above. As I understand it, this is a type on its own and not a Symmetric-view on the BCSR type, is it?

I agree, I changed it.
Yes SparseMatrixBCSRSymmetric is a new type for symmetric matrices, it exploits the fact that for symmetric matrices only the lower triangular blocks have to be stored and the corresponding upper triangular part can just be constructed by mirroring along the diagonal.
Not sure if it's worth it to save the diagonal elements as LowerTriangulars, as they need more memory than normal arrays:

julia> Base.summarysize(randn(100,100)) 80000 julia> Base.summarysize(LowerTriangular(randn(100,100))) 80008

Your thoughts on that?

As discussed here, the point of LowerTriangular seems not to be memory reduction but instead runtime improvements in BLAS. Since BLAS requires a certain memory layout, this makes it necessary to save the full matrix (as you noticed).

I have just done a few timings on my machine. It seems, that the use of LowerTriangular does save one third to half of the matrix multiplication time (depending on the size of the matrix blocks, 100x100 - 1000x1000)*. But since your matrix is really symmetric and not triangular, this can't be sufficient since you have to do two multiplications if you choose LowerTriangulars.

Another approach would be to use Symmetric views on the diagonal matrix blocks. But as Symmetric(matrix) seems not to be exploited and the underlying data is still the full matrix, these timings don't really improve.
julia> r = Matrix(Symmetric(rand(1000,1000)));
julia> s = Symmetric(r);
julia> rr = zeros(r); @time for i in 1:100 rr += r*r end
24.662503 seconds (400 allocations: 1.490 GiB, 0.25% gc time)
julia> rs = zeros(r); @time for i in 1:100 rs += r*s end
24.902586 seconds (1.89 k allocations: 2.235 GiB, 0.34% gc time)
julia> sr = zeros(r); @time for i in 1:100 sr += s*r end
27.563856 seconds (7.70 k allocations: 1.491 GiB, 0.24% gc time)
julia> ss = zeros(r); @time for i in 1:100 ss += s*s end
26.048422 seconds (1.58 k allocations: 2.980 GiB, 0.58% gc time)

* results for LowerTriangular
julia> l = LowerTriangular(r);
julia> ll = zeros(r); @time for i in 1:100 ll += l*l end
13.332010 seconds (700 allocations: 2.235 GiB, 0.50% gc time)
julia> rl = zeros(r); @time for i in 1:100 rl += r*l end
12.494506 seconds (2.68 k allocations: 1.490 GiB, 0.46% gc time)
julia> lr = zeros(r); @time for i in 1:100 lr += l*r end
12.227783 seconds (2.02 k allocations: 1.490 GiB, 0.45% gc time)

Thank you for the explanation, I'll leave that feature in for now. While matrix multiplications with LowerTriangulars are faster I am not convinced that that is much of a benefit as only blocks on the main diagonal get stored as LowerTriangular matrices.
I'll think about a way of using the fact that the corresponding dense matrix is symmetric in the future.

vonDonnerstein · 2018-07-19T15:06:36Z

src/SubModules/MatrixModule.jl

+type BCSRSpM
+	val::Array{Array{Float64,2},1}
+	col::Array{Int,1}
+	rowptr::Array{Int,1}


As this is a very specific nomenclature, I would suggest adding a reference to the documentation, where the details of the implementation (with these names) can be found, or - if no such publication exists - explain the basic idea with a few minimal examples in the documentation. You could have a look, e.g., at the documentation for Symmetric to get an idea.

Added a small example in the documentation, also did that for other BCSRSpM constructors (Note that BCSRSpM is now called SparseMatrixBCSR).

vonDonnerstein · 2018-07-19T15:51:15Z

src/SubModules/MatrixModule.jl

+    computeDifferenceRowptr(Integer,rowptr)
+returns the running indices to run over all blocks in a row; b-a = # of blocks in a row
+"""
+function computeDifferenceRowptr(i::Int,rowptr)


for this and many other functions of this file, I'm still unclear about how "internal" they are meant to be. While I don't get the idea of this function at all, maybe that's not a problem, if it isn't ever to be called by any user of the library. If that is your intent, maybe you could group these "helper" functions in the beginning or end of the file (e.g., putting "# HELPERS" before them)?

I reordered the code and made sections SPARSE MATRIX CONSTRUCTORS, HELPERS, BASE EXTENSIONS and MATRIX OPERATIONS for better readability. computeDifferenceRowptr is not supposed to be called by the average user just a helper function.

vonDonnerstein · 2018-07-19T15:52:27Z

src/SubModules/MatrixModule.jl

+    multiplySpMV(SpM,Vector)
+returns the product Vector2 = SpM*Vector
+"""
+function multiplySpMV(SpM::BCSRSpM,vec::Array{Float64,1})


isn't this an extension to Base.*?

Yes, I added the extension in the code as *(SpM::SparseMatrixBCSR,vec::Array{Float64,1}) = multiplySpMV(SpM,vec).

Why not just rename this function?

function Base.*(SpM::SparseMatrixBCSR, vec::Vector{Float64}) if SpM.colpattern[... end

Maybe I'm just missing the advantage of having it as a separate function. I've noticed that this style of naming a separate function does appear in Base in some places, but never understood what it's good for, so I would be eager to learn.

Honestly I didn't know that this was an option, I see no point of having it as a separate function I just thought that was the way to do it. I changed it. You don't actually need Base.* as a function name just * is enough if Base.* is imported.

vonDonnerstein · 2018-07-19T15:55:42Z

src/SubModules/MatrixModule.jl

+	return M
+end
+
+#TODO: macht aehnliches wie computeBlock -> vereinigen


I'll wait with my comments on those functions for which you still have TODOs until you are ready.

vonDonnerstein · 2018-07-19T16:03:18Z

src/SubModules/MatrixModule.jl

+	return SpMT
+end
+
+#@Henry: laesst sich so nicht in MatrixModule aufrufen weil ShellModule benutzt


I'm not sure if benchmarking should be a part of QuantumLab itself. Maybe, if you want to, we could have a separate benchmark.jl script in QL/test to separate it from the code functionalities. Or you could keep that script separate (not commiting it to QuantumLab and just run it when you are ready and report the results in your lab report.

I'll just leave it out, it's only meant for my report anyway.

vonDonnerstein · 2018-07-19T16:05:13Z

src/SubModules/ShellModule.jl

+    computePattern(basis)
+returns the matrix blocking pattern for a given basis
+"""
+function computePattern(basis::Array{ShellModule.Shell,1})


basis is a bit misleading for the argument name, because it isn't a Basis (e.g. GaussianBasis). How about shells::Vector{Shell}?

Good point I changed all instances of basis to shells.

vonDonnerstein · 2018-07-21T08:28:38Z

src/SubModules/MatrixModule.jl

+returns a matrix M in BCSR format with the given blocking row- and colpattern
+"""
+function BCSRSpM(M::Array{Float64,2},rowpattern::Array{Tuple{Int64,Int64},1},colpattern::Array{Tuple{Int64,Int64},1})
+	(checkPattern(M,rowpattern,1) && checkPattern(M,colpattern,2)) || error("Hey idiot your pattern doesn't match the matrix ᕦ(ò_óˇ)ᕤ")


Oh, I must have missed this on my first read. While I, too, enjoy humor in error messages from time to time, I very much prefer if they stay polite and clear. When the user uses the library in an incorrect way, the fact that it doesn't work as expected is already quite frustrating. Being insulted by the code doesn't really help with that ;).
I've added a corresponding bullet point to the coding guidelines.

Good point I wouldn't want to receive such an error message either tbh. I changed the error messages to be more informative and polite but I left the ASCII smileys in.

vonDonnerstein · 2018-07-21T11:52:28Z

Felix, after you have made your changes, could you update your pull request? See the last point of Detailed Description for how to do that. No hurry, just take the time that you need. Thanks.

…new code in runtests.jl

added sparse algebra to MatrixModule.jl, modified ShellModule.jl to i…

e6a0d41

…nclude MatrixModule.jl

vonDonnerstein requested changes Jul 19, 2018

View reviewed changes

vonDonnerstein reviewed Jul 21, 2018

View reviewed changes

FelixBangerter added 2 commits July 23, 2018 08:26

reworked MatrixModule.jl and ShellModule.jl, added test coverage for …

db41aca

…new code in runtests.jl

reworked MatrixModule.jl and ShellModule.jl, added test coverage for …

a8cc399

…new code in runtests.jl

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

SparseAlgebra #12

SparseAlgebra #12

FelixBangerter commented Jul 19, 2018

coveralls commented Jul 19, 2018 •

edited

Loading

vonDonnerstein left a comment

vonDonnerstein Jul 19, 2018

vonDonnerstein Jul 19, 2018

FelixBangerter Jul 20, 2018

vonDonnerstein Jul 19, 2018

FelixBangerter Jul 20, 2018

vonDonnerstein Jul 19, 2018

FelixBangerter Jul 20, 2018

vonDonnerstein Jul 19, 2018

FelixBangerter Jul 20, 2018

vonDonnerstein Jul 21, 2018

FelixBangerter Jul 22, 2018

vonDonnerstein Jul 19, 2018

FelixBangerter Jul 20, 2018 •

edited

Loading

vonDonnerstein Jul 19, 2018

FelixBangerter Jul 20, 2018

vonDonnerstein Jul 19, 2018

FelixBangerter Jul 20, 2018

vonDonnerstein Jul 21, 2018

FelixBangerter Jul 22, 2018

vonDonnerstein Jul 19, 2018

vonDonnerstein Jul 19, 2018

FelixBangerter Jul 20, 2018

vonDonnerstein Jul 19, 2018

FelixBangerter Jul 20, 2018

vonDonnerstein Jul 21, 2018

FelixBangerter Jul 22, 2018

vonDonnerstein commented Jul 21, 2018



		"""
		BCSRSpM(val,col,rowptr,rowpattern,colpattern)

SparseAlgebra #12

Are you sure you want to change the base?

SparseAlgebra #12

Conversation

FelixBangerter commented Jul 19, 2018

coveralls commented Jul 19, 2018 • edited Loading

vonDonnerstein left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

FelixBangerter Jul 20, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

vonDonnerstein commented Jul 21, 2018

coveralls commented Jul 19, 2018 •

edited

Loading

FelixBangerter Jul 20, 2018 •

edited

Loading