Merge pull request #75 from matrixfunctions/multilincomb_refactor

More than two matrices in lincomb

Author: jarlebring

src/GraphMatFun.jl

@@ -25,6 +25,7 @@ include("generators/horner.jl")
 include("generators/newton_schulz.jl")
 include("generators/exp.jl")
 include("generators/sastre.jl")
+include("generators/bigraph.jl")
 
 # Error bounds
 include("error_bounds.jl")

src/code_gen/gen_c_code.jl

@@ -107,6 +107,12 @@ function get_blas_type(::LangC_OpenBLAS, T::Type{Complex{Float64}})
     return ("openblas_complex_double", "z")
 end
 
+
+function preprocess_codegen(graph, lang::LangC)
+    (g,_)=graph_bigraph(graph)
+    return g # Lang C only supports bigraphs
+end
+
 function function_definition(lang::LangC, graph, T, funname, precomputed_nodes)
     (blas_type, blas_prefix) = get_blas_type(lang, T)
     code = init_code(lang)

src/code_gen/gen_code.jl

@@ -1,8 +1,6 @@
 include("gen_code_mem.jl")
 include("gen_code_snippets.jl")
 
-include("multilincomb.jl");
-
 include("gen_c_code.jl")
 include("gen_julia_code.jl")
 include("gen_matlab_code.jl")

src/code_gen/gen_julia_code.jl

@@ -4,7 +4,6 @@ export LangJulia
 struct LangJulia
     overwrite_input::Any # Overwrite input
     inline::Any
-    dot_fusing::Any # Allow dot fusion
     axpby_header::Any
     alloc_function
     only_overwrite
@@ -14,18 +13,14 @@ end
 default_alloc_function(k)="similar(A,T)"
 
 """
-    LangJulia(overwrite_input=true,inline=true,dot_fusing=true,axpby_header=:auto,alloc_function,only_overwrite=false)
+    LangJulia(overwrite_input=true,inline=true,axpby_header=:auto,alloc_function,only_overwrite=false)
 
-Code generation in julia language, with optional overwriting of input, inlining
-the function and optional usage of dot fusion. The `axpby_header` specifies if axpby function calls should be included in the beginning of the file. The parameter `alloc_function` is a function of three parameters `alloc_function(k)` where `k` is the memory slot (default is `alloc_function(k)=similar(A,T)`). The `only_overwrite` specifies if `f` should be created if the overwrite funtion `f!` contains the actual code.
+Code generation in julia language, with optional overwriting of input. The `axpby_header` specifies if axpby function calls should be included in the beginning of the file. The parameter `alloc_function` is a function of three parameters `alloc_function(k)` where `k` is the memory slot (default is `alloc_function(k)=similar(A,T)`). The `only_overwrite` specifies if `f` should be created if the overwrite funtion `f!` contains the actual code.
 """
-LangJulia() = LangJulia(true, true, true, :auto, default_alloc_function,false,"ValueOne","matfun_axpby!")
+LangJulia() = LangJulia(true, true,  :auto, default_alloc_function,false,"ValueOne","matfun_axpby!")
 LangJulia(overwrite_input) = LangJulia(overwrite_input, true, true, :auto, default_alloc_function,false,"ValueOne","matfun_axpby!")
-function LangJulia(overwrite_input, inline)
-    return LangJulia(overwrite_input, inline, true, :auto, default_alloc_function,false,"ValueOne","matfun_axpby!")
-end
-function LangJulia(overwrite_input, inline, dot_fusing; value_one_name="ValueOne",axpby_name="matfun_axpby!")
-    return LangJulia(overwrite_input, inline, dot_fusing, :auto, default_alloc_function,false,value_one_name,axpby_name)
+function LangJulia(overwrite_input, inline ; value_one_name="ValueOne",axpby_name="matfun_axpby!")
+    return LangJulia(overwrite_input, inline,  :auto, default_alloc_function,false,value_one_name,axpby_name)
 end
 
 # Language specific operations.
@@ -60,11 +55,7 @@ function assign_coeff_basic(lang::LangJulia, v, i)
 end
 
 function preprocess_codegen(graph, lang::LangJulia)
-    if (lang.dot_fusing)
-        return MultiLincombCompgraph(graph) # Merge many lincombs for dot fusion
-    else
-        return graph
-    end
+    return graph # Merge many lincombs for dot fusion
 end
 
 # Code generation.
@@ -124,7 +115,7 @@ function function_definition(
     code = init_code(lang)
     axpby_header = lang.axpby_header
     if (lang.axpby_header == :auto)
-        axpby_header = !lang.dot_fusing
+        axpby_header = false
     end
 
     push_code!(code, "using LinearAlgebra", ind_lvl = 0)
@@ -297,7 +288,7 @@ end
 function execute_operation!(
     lang::LangJulia,
     T,
-    graph::MultiLincombCompgraph,
+    graph,
     node,
     dealloc_list,
     mem,
@@ -399,10 +390,6 @@ function execute_operation!(
     end
 end
 
-function execute_operation!(lang::LangJulia, T, graph, node, dealloc_list, mem)
-    return execute_operation_basic!(lang, T, graph, node, dealloc_list, mem)
-end
-
 # The general base case. Separated for dispatch.
 function execute_operation_basic!(
     lang::LangJulia,

src/code_gen/gen_matlab_code.jl

@@ -61,21 +61,26 @@ end
 function execute_operation!(lang::LangMatlab, T, graph, node, dealloc_list, mem)
     op = graph.operations[node]
     parent1 = graph.parents[node][1]
-    parent2 = graph.parents[node][2]
 
     code = init_code(lang)
     push_comment!(code, "Computing $node with operation: $op")
     # Needs to be updated
     if op == :mult
+        parent2 = graph.parents[node][2]
         push_code!(code, "$node = $parent1 * $parent2;")
     elseif op == :ldiv
+        parent2 = graph.parents[node][2]
         push_code!(code, "$node = $parent1 \\ $parent2;")
     elseif op == :lincomb
-        (coeff1, coeff1_code) = assign_coeff(lang, graph.coeffs[node][1], 1)
-        push_code!(code, "$coeff1_code;")
-        (coeff2, coeff2_code) = assign_coeff(lang, graph.coeffs[node][2], 2)
-        push_code!(code, "$coeff2_code;")
-        push_code!(code, "$node = $coeff1*$parent1 + $coeff2*$parent2;")
+
+        coeff_names=Vector();
+        for (i,coeff) = enumerate(graph.coeffs[node])
+            (coeff, coeff_code) = assign_coeff(lang, graph.coeffs[node][i], i)
+            push_code!(code, "$coeff_code;")
+            push!(coeff_names,coeff)
+        end
+        sum_code=join((coeff_names.*"*") .* string.(graph.parents[node])," + ")
+        push_code!(code, "$node = $sum_code;")
     end
     return (code, "$node")
 end