From 74ca269e72007152aec1ea5925e047b341483519 Mon Sep 17 00:00:00 2001
From: Jukka Aho <jukka.aho@kapsi.fi>
Date: Wed, 2 Dec 2015 08:36:17 +0200
Subject: [PATCH] little tweaking of assembly

---
 src/assembly.jl     | 25 +++++++++++++++++--------
 src/directsolver.jl |  8 ++++----
 2 files changed, 21 insertions(+), 12 deletions(-)

diff --git a/src/assembly.jl b/src/assembly.jl
index 39061b6..b1dfb92 100644
--- a/src/assembly.jl
+++ b/src/assembly.jl
@@ -47,9 +47,10 @@ function reduce(assembly::Assembly, boundary_dofs_::Vector{Int})
     # empty assembly to release memory for factorization
     empty!(assembly.stiffness_matrix)
     empty!(assembly.force_vector)
+    gc()
 
     if dim < 100000
-        # no need to do any reduction of matrix size at all
+        # no need to do any reduction of matrix size at all, just \ it.
         return CAssembly([], all_dofs, Matrix{Float64}(), K, f, spzeros(0, 0), spzeros(0,1))
     end
 
@@ -64,9 +65,10 @@ function reduce(assembly::Assembly, boundary_dofs_::Vector{Int})
     fb = f[boundary_dofs]
 
     F = cholfact(K[interior_dofs, interior_dofs])
-    K = spzeros(0, 0)
+    K = 0
+    gc()
+
 
-#=
     if dim < 100000
         # for small problems we don't need to care about memory usage
         Kd = Kib' * (F \ Kib)
@@ -78,31 +80,38 @@ function reduce(assembly::Assembly, boundary_dofs_::Vector{Int})
         for bi in 1:nb
             mod(bi, p) == 0 && info("Reduction: ", round(Int, bi/nb*100), " % done")
             C = full(F \ Kib[:, bi])
-            for bj in 1:nb
+            for bj in bi:nb
                 d = Kib[:, bj]
-                Kd[bj,bi] = dot(C[rowvals(d)], nonzeros(d))
+                @inbounds Kd[bj,bi] = dot(C[rowvals(d)], nonzeros(d))
             end
         end
+        Kd += tril(Kd, -1)'
     end
     Kc = spzeros(dim, dim)
     Kc[boundary_dofs, boundary_dofs] = Kbb - Kd
-=#
 
+#=  # this is slightly faster but uses more memory
     chunks = round(Int, dim/3000)
     info("Reduction is done in $chunks chunks.")
     nb = length(boundary_dofs)
     kk = round(Int, collect(linspace(0, nb, chunks+1)))
     sl = [kk[j]+1:kk[j+1] for j=1:length(kk)-1]
-    Kc = spzeros(dim, dim)
+    Kd = zeros(Float64, nb, nb)
+    #Kd = SharedArray(Float64, nb, nb)
     for (k,sli) in enumerate(sl)
         b1 = boundary_dofs[sli]
         Sc = F \ Kib[:,sli]
         for slj in sl
             b2 = boundary_dofs[slj]
-            Kc[b2,b1] = Kbb[slj,sli] - Kib[:,slj]'*Sc
+            #Kc[b2,b1] = Kbb[slj,sli] - Kib[:,slj]'*Sc
+            Kd[slj, sli] = Kib[:,slj]'*Sc
         end
         info("Reduction: ", round(k/chunks*100, 0), " % done")
     end
+    
+    Kc = spzeros(dim, dim)
+    Kc[boundary_dofs, boundary_dofs] = Kbb - Kd
+=#
 
     fc = spzeros(dim, 1)
     fc[boundary_dofs] = fb - Kib' * (F \ fi)
diff --git a/src/directsolver.jl b/src/directsolver.jl
index e32838f..4148dfb 100644
--- a/src/directsolver.jl
+++ b/src/directsolver.jl
@@ -163,16 +163,16 @@ function call(solver::DirectSolver, ::Type{Val{:noreduce}}, time::Number=0.0)
         tic(timing, "solution of system")
         # solve increment for linearized problem
         nz = unique(rowvals(A))  # take only non-zero rows
-        sol = zeros(b)
-        sol[nz] = A[nz,nz] \ full(b[nz])
+        sol = zeros(length(b))
+        sol[nz] = full(A[nz,nz]) \ full(b[nz])
         #info("solution vector before reconstruction")
         #info(full(sol)')
         la = sol[dim+1:end]
         #for assembly in assemblies
         #    reconstruct!(assembly, sol)
         #end
-        la = vec(full(la))
-        sol = vec(full(sol))
+#       la = vec(full(la))
+#       sol = vec(full(sol))
         #info("la = ", la')
         #info("sol = ", sol')