diff --git a/src/code_gen/gen_c_code.jl b/src/code_gen/gen_c_code.jl
index fdb7e10..a1410e9 100644
--- a/src/code_gen/gen_c_code.jl
+++ b/src/code_gen/gen_c_code.jl
@@ -4,7 +4,7 @@ export LangC_MKL, LangC_OpenBLAS
 """
      LangC_MKL(gen_main::Bool)
 
-Code generation in C using MKL implementation of BLAS.
+Code generation in C using the oneAPI MKL implementation of BLAS.
 """
 struct LangC_MKL
     gen_main::Any
@@ -14,7 +14,7 @@ end
 """
      LangC_OpenBLAS(gen_main::Bool)
 
-Code generation in C using OpenBLAS implementation of BLAS.
+Code generation in C using the OpenBLAS implementation of BLAS.
 """
 struct LangC_OpenBLAS
     gen_main::Any
@@ -32,56 +32,58 @@ end
 
 # Language specific operations.
 comment(::LangC, s) = "/* $s */"
-
 slotname(::LangC, i) = "memslots[$(i-1)]"
 
 # Variable declaration, initialization, and reference.
-# In C, complex types are structures and are passed by reference.
-function assign_coeff(::LangC_MKL, val::T, i) where {T<:Complex}
-    assignment_string = [
-        "coeff$i.real = " * string(real(val)) * ";",
-        "coeff$i.imag = " * string(imag(val)) * ";",
-    ]
-    variable_string = "&coeff$i"
-    return (variable_string, assignment_string)
-end
-function assign_coeff(::LangC_OpenBLAS, val::T, i) where {T<:Complex}
-    assignment_string =
-        ["coeff$i = " * string(real(val)) * " + " * string(imag(val)) * "*I;"]
-    variable_string = "&coeff$i"
-    return (variable_string, assignment_string)
-end
-function assign_coeff(::LangC, val::T, i) where {T<:Real}
-    assignment_string = ["coeff$i = $val;"]
-    variable_string = "coeff$i"
-    return (variable_string, assignment_string)
+# In MKL, complex types are structures and are passed by reference.
+function initialization_string(::LangC_MKL, real_part, imag_part)
+    return "{.real = " * string(real_part) * ", " *
+            ".imag = " * string(imag_part) * "}"
+end
+function initialization_string(::LangC_OpenBLAS, real_part, imag_part)
+    return string(real_part) * " + " * string(imag_part) * "*I"
+end
+function initialization_string(lang::LangC, val::T) where {T<:Complex}
+    return initialization_string(lang, real(val), imag(val))
+end
+function initialization_string(lang::LangC, val::T) where {T<:Real}
+    return "$val"
+end
+
+function declare_var(lang::LangC, val, id, type)
+    return "$type $id = " * initialization_string(lang, val) * ";"
+end
+function declare_coeff(lang::LangC, val::T, id) where T
+    # Use real type regardless of T.
+    actual_val = isreal(val) ? real(val) : val
+    (blas_t, blas_prefix) = get_blas_type(lang, typeof(actual_val))
+    variable_string = "coeff_$id"
+    dec_init_string = declare_var(lang, actual_val, variable_string, blas_t)
+    return (val, variable_string, dec_init_string)
+end
+
+function assignment_string(::LangC, var, real_part)
+    return "$var = $(real_part)"
+end
+function assignment_string(::LangC_MKL, var, real_part, imag_part)
+    return "$(string(var)).real = $(real_part); " *
+           "$(string(var)).imag = $(imag_part);"
+end
+function assignment_string(lang::LangC_OpenBLAS, var, real_part, imag_part)
+    return "$(string(var)) = " *
+            initialization_string(lang, real_part, imag_part) * ";"
 end
 
 # Constant declaration and reference.
-function declare_constant(::LangC_MKL, val::T, id, type) where {T<:Complex}
-    return "const $type $id = {.real = " *
-           string(real(val)) *
-           ", " *
-           ".imag = " *
-           string(imag(val)) *
-           "};"
-end
-function declare_constant(::LangC_OpenBLAS, val::T, id, type) where {T<:Complex}
-    return "const $type $id = " *
-           string(real(val)) *
-           " + " *
-           string(imag(val)) *
-           "*I;"
-end
-function declare_constant(::LangC, val::T, id, type) where {T<:Real}
-    return "const $type $id = $val;"
-end
-reference_constant(::LangC, T, id) = (T <: Complex) ? "&$id" : "$id"
+function declare_const(lang::LangC, val, id, type)
+    return "const " * declare_var(lang::LangC, val, id, type)
+end
+reference_value(::LangC, T, id) = (T <: Complex) ? "&$id" : "$id"
 
 # Code generation.
 # Return MKL/OpenBLAS-specific includes.
 function get_blas_includes(::LangC_MKL)
-    return "#include<mkl/mkl.h>"
+    return "#include<mkl.h>"
 end
 function get_blas_includes(::LangC_OpenBLAS)
     return "#include<cblas.h>\n#include<lapacke.h>"
@@ -101,16 +103,80 @@ function get_blas_type(::LangC_MKL, T::Type{Complex{Float64}})
     return ("MKL_Complex16", "z")
 end
 function get_blas_type(::LangC_OpenBLAS, T::Type{Complex{Float32}})
-    return ("openblas_complex_float", "c")
+    return ("lapack_complex_float", "c")
 end
 function get_blas_type(::LangC_OpenBLAS, T::Type{Complex{Float64}})
-    return ("openblas_complex_double", "z")
+    return ("lapack_complex_double", "z")
 end
 
+# Add functions to work with MKL complex types.
+function add_auxiliary_functions(code, ::LangC, T)
+    return # In general, nothing is needed.
+end
+function create_auxiliary_functions(complex_type, real_type)
+    return "// Compute c <- c + a * b.
+void fma_$(complex_type)_$(real_type)($(complex_type) *c,
+                    const $(real_type) *a,
+                    const $(complex_type) *b) {
+    c->real += *a * b->real;
+    c->imag += *a * b->imag;
+}
+
+// Compute c <- a * b.
+void prod_$(complex_type)_$(real_type)($(complex_type) *c,
+                    const $(real_type) *a,
+                    const $(complex_type) *b) {
+    c->real = *a * b->real;
+    c->imag = *a * b->imag;
+}
+
+// Compute c <- c + a.
+void acc_$(complex_type)_$(real_type)($(complex_type) *c,
+                    const $(real_type) *a) {
+        c->real += *a;
+}
+
+// Compute c <- c + a * b.
+void fma_$(complex_type)($(complex_type) *c,
+                    const $(complex_type) *a,
+                    const $(complex_type) *b) {
+    $(real_type) real_part = a->real * b->real - a->imag * b->imag;
+    $(real_type) imag_part = a->real * b->imag + a->imag * b->real;
+    c->real += real_part;
+    c->imag += imag_part;
+}
+
+// Compute c <- a * b.
+void prod_$(complex_type)($(complex_type) *c,
+                    const $(complex_type) *a,
+                    const $(complex_type) *b) {
+    $(real_type) real_part = a->real * b->real - a->imag * b->imag;
+    $(real_type) imag_part = a->real * b->imag + a->imag * b->real;
+    c->real = real_part;
+    c->imag = imag_part;
+}
+
+// Compute c <- c + a.
+void acc_$(complex_type)($(complex_type) *c,
+                    const $(complex_type) *a) {
+    c->real += a->real;
+    c->imag += a->imag;
+}"
+end
+function add_auxiliary_functions(code, ::LangC_MKL, T::Type{Complex{Float32}})
+    push_code!(
+        code,
+        create_auxiliary_functions("MKL_Complex8", "float"),
+        ind_lvl = 0
+    )
+end
 
-function preprocess_codegen(graph, lang::LangC)
-    (g,_)=graph_bigraph(graph)
-    return g # Lang C only supports bigraphs
+function add_auxiliary_functions(code, ::LangC_MKL, T::Type{Complex{Float64}})
+    push_code!(
+        code,
+        create_auxiliary_functions("MKL_Complex16", "double"),
+        ind_lvl = 0
+    )
 end
 
 function function_definition(lang::LangC, graph, T, funname, precomputed_nodes)
@@ -121,6 +187,11 @@ function function_definition(lang::LangC, graph, T, funname, precomputed_nodes)
     push_code!(code, "#include<stdlib.h>", ind_lvl = 0)
     push_code!(code, "#include<string.h>", ind_lvl = 0)
     push_code!(code, "")
+    # Add auxiliary functions to work with MKL complex types in lincomb.
+    if :lincomb ∈ values(graph.operations)
+        add_auxiliary_functions(code, lang, T)
+    end
+    push_code!(code, "")
     push_comment!(code, "Code for matrix function evaluation.", ind_lvl = 0)
     input_variables =
         join(map(x -> "$blas_type *" * string(x), precomputed_nodes), ", ")
@@ -145,7 +216,7 @@ function init_mem(lang::LangC, max_nof_nodes, precomputed_nodes)
 end
 
 function function_init(lang::LangC, T, mem, graph, precomputed_nodes)
-    (blas_type, blas_prefix) = get_blas_type(lang, T)
+    (blas_t, blas_prefix) = get_blas_type(lang, T)
     code = init_code(lang)
     max_nodes = size(mem.slots, 1)
 
@@ -162,26 +233,14 @@ function function_init(lang::LangC, T, mem, graph, precomputed_nodes)
     push_code!(code, "")
 
     push_comment!(code, "Declarations and initializations.")
-    # Allocate coefficients for linear combinations, if needed.
-    if :lincomb in graph_ops
-        push_code!(code, "$blas_type coeff1, coeff2;")
-    end
     # Declare constants ZERO and ONE, if needed.
     if :mult in graph_ops
-        push_code!(
-            code,
-            declare_constant(lang, convert(T, 0.0), "ZERO", blas_type),
-        )
-    end
-    if :lincomb in graph_ops || :mult in graph_ops
-        push_code!(
-            code,
-            declare_constant(lang, convert(T, 1.0), "ONE", blas_type),
-        )
+        push_code!(code, declare_const(lang, convert(T, 0.0), "ZERO", blas_t))
+        push_code!(code, declare_const(lang, convert(T, 1.0), "ONE", blas_t))
     end
     # Array ipiv for pivots of GEPP (needed only if graph has linear systems).
     if :ldiv in graph_ops
-        push_code!(code, "lapack_int *ipiv = malloc(n*sizeof(*ipiv));")
+        push_code!(code, "lapack_int *ipiv = malloc(n * sizeof(*ipiv));")
     end
     push_code!(code, "size_t j;")
     push_code!(code, "")
@@ -189,12 +248,12 @@ function function_init(lang::LangC, T, mem, graph, precomputed_nodes)
     push_comment!(code, "Memory management.")
     push_code!(
         code,
-        "$blas_type *master_mem = malloc(n*n*max_memalloc" *
-        "*sizeof(*master_mem));",
+        "$blas_t *master_mem = malloc(n * n * max_memalloc" *
+        " * sizeof(*master_mem));",
     )
     push_code!(
         code,
-        "$blas_type *memslots[max_memslots]; " *
+        "$blas_t *memslots[max_memslots]; " *
         comment(lang, "As many slots as nodes"),
     )
 
@@ -211,26 +270,18 @@ function function_init(lang::LangC, T, mem, graph, precomputed_nodes)
             # Otherwise make a copy
             push_code!(
                 code,
-                "memcpy($Ak_slot_name,$n,n*n*sizeof(*master_mem));",
+                "memcpy($Ak_slot_name, $n, n * n * sizeof(*master_mem));",
             )
         end
     end
 
     # Initialize pointers.
     push_comment!(code, "The other slots are pointers to allocated memory.")
-    push_code!(code, "for (j=$(num_preallocated_slots-1); j<max_memalloc; j++)")
-    push_code!(code, "memslots[j+1] = master_mem+j*n*n;", ind_lvl = 2)
-
-    # Store identity explicitly only if graph has a linear combination of I.
-    if has_identity_lincomb(graph)
-        push_comment!(code, "Graph has linear combination of identities.")
-        push_comment!(code, "The matrix I is explicitly allocated.")
-        alloc_slot!(mem, num_precomputed_nodes + 1, :I)
-        nodemem = get_slot_name(mem, :I)
-        push_code!(code, "memset($nodemem, 0, n*n*sizeof(*master_mem));")
-        push_code!(code, "for(j=0; j<n*n; j+=n+1)")
-        push_code!(code, "*($nodemem+j) = ONE;", ind_lvl = 2)
-    end
+    push_code!(
+        code,
+        "for (j = $(num_preallocated_slots-1); j < max_memalloc; j++)"
+    )
+    push_code!(code, "memslots[j+1] = master_mem + j * n * n;", ind_lvl = 2)
 
     return code
 end
@@ -241,7 +292,7 @@ function function_end(lang::LangC, graph, mem)
     retval = get_slot_name(mem, retval_node)
     push_code!(code, "")
     push_comment!(code, "Prepare output.")
-    push_code!(code, "memcpy(output, $retval, n*n*sizeof(*output));")
+    push_code!(code, "memcpy(output, $retval, n * n * sizeof(*output));")
     if :ldiv in values(graph.operations)
         push_code!(code, "free(ipiv);")
     end
@@ -252,41 +303,81 @@ function function_end(lang::LangC, graph, mem)
     return code
 end
 
+function add_lincomb_body(code, lang::LangC, T, nodemem,
+                          coeff_names, coeff_values, parent_mems)
+    rhs = join((coeff_names .* " * ") .* ("*(" .* parent_mems .* " + i)"),
+            " + \n            ")
+    for_body = "*(" * nodemem * " + i) = " * (isempty(rhs) ? "0" : rhs) * ";"
+    push_code!(code, for_body; ind_lvl = 2)
+end
+
+prod_function_name(T::Type{Complex{Float32}}, coeff_real) =
+    coeff_real ? "prod_MKL_Complex8_float" : "prod_MKL_Complex8"
+prod_function_name(T::Type{Complex{Float64}}, coeff_real) =
+    coeff_real ? "prod_MKL_Complex16_double" : "prod_MKL_Complex16"
+fma_function_name(T::Type{Complex{Float32}}, coeff_real) =
+    coeff_real ? "fma_MKL_Complex8_float" : "fma_MKL_Complex8"
+fma_function_name(T::Type{Complex{Float64}}, coeff_real) =
+    coeff_real ? "fma_MKL_Complex16_double" : "fma_MKL_Complex16"
+function add_lincomb_body(code, lang::LangC_MKL, T::Type{Complex{S}}, nodemem,
+                          coeff_names, coeff_values,
+                          parent_mems) where S <: Real
+    for (it, (coeff, coeff_value, parent)) in
+        enumerate(zip(coeff_names, coeff_values, parent_mems))
+        operation = it == 1 ?
+                    prod_function_name(T, isreal(coeff_value)) :
+                    fma_function_name(T, isreal(coeff_value))
+        statement = operation *
+                "(" * nodemem * " + i, " *
+                reference_value(lang, T, coeff) * ", " *
+                parent * " + i);"
+        push_code!(code, statement; ind_lvl = 2)
+    end
+end
+
+function add_lincomb_identity_body(code, lang::LangC, T,
+                                   nodemem, coeff_id, coeff_value)
+    push_code!(
+        code,
+        "*(" * nodemem * " + i) += " * coeff_id * ";",
+        ind_lvl = 2
+    )
+end
+
+acc_function_name(T::Type{Complex{Float32}}, coeff_real) =
+    coeff_real ? "acc_MKL_Complex8_float" : "acc_MKL_Complex8"
+acc_function_name(T::Type{Complex{Float64}}, coeff_real) =
+    coeff_real ? "acc_MKL_Complex16_double" : "acc_MKL_Complex16"
+function add_lincomb_identity_body(code, lang::LangC_MKL, T::Type{Complex{S}},
+                                   nodemem, coeff_id, coeff_value) where S <: Real
+    statement = acc_function_name(T, isreal(coeff_value)) * "(" * nodemem * " + i, " *
+        reference_value(lang, T, coeff_id) * ");"
+    push_code!(code, statement, ind_lvl = 2)
+end
+
 function execute_operation!(lang::LangC, T, graph, node, dealloc_list, mem)
     (blas_type, blas_prefix) = get_blas_type(lang, T)
 
     op = graph.operations[node]
-    parent1 = graph.parents[node][1]
-    parent2 = graph.parents[node][2]
+
+    code = init_code(lang)
+    push_code!(code, "")
+    push_comment!(code, "Computing $node with operation: $op")
 
     # Keep deallocation list (used for smart memory management).
     dealloc_list = deepcopy(dealloc_list)
     setdiff!(dealloc_list, keys(mem.special_names))
 
-    # Check if parents have a memory slot.
-    if has_identity_lincomb(graph) # Linear combination of I.
-        p1_is_identity = p2_is_identity = false
-    else
-        p1_is_identity = parent1 == :I
-        p2_is_identity = parent2 == :I
-        setdiff!(dealloc_list, [:I])
-    end
+    rzero = reference_value(lang, T, "ZERO")
+    rone = reference_value(lang, T, "ONE")
 
-    # Get memory slots of parents, if explicitly stored.
-    if !p1_is_identity
+    if op == :mult
+        # The graph has been compressed, thus we can assume that neither parent1
+        # is the identity.
+        parent1 = graph.parents[node][1]
+        parent2 = graph.parents[node][2]
         parent1mem = get_slot_name(mem, parent1)
-    end
-    if !p2_is_identity
         parent2mem = get_slot_name(mem, parent2)
-    end
-
-    rzero = reference_constant(lang, T, "ZERO")
-    rone = reference_constant(lang, T, "ONE")
-
-    code = init_code(lang)
-    push_code!(code, "")
-    push_comment!(code, "Computing $node with operation: $op")
-    if op == :mult
         (nodemem_i, nodemem) = get_free_slot(mem)
         alloc_slot!(mem, nodemem_i, node)
         push_code!(
@@ -298,6 +389,11 @@ function execute_operation!(lang::LangC, T, graph, node, dealloc_list, mem)
             "            $rzero, $nodemem, n);",
         )
     elseif op == :ldiv
+        # Left parent cannot be the identity.
+        parent1 = graph.parents[node][1]
+        parent1mem = get_slot_name(mem, parent1)
+        # Parent 2 can be the identity, but this is dealt with below.
+        parent2 = graph.parents[node][2]
         # Find a memory slot for the LU factors.
         # As ?detrs computes the LU decomposition of parent1 in-place, we need
         # to make a copy of $parent1mem, unless parent1 is on the deallocation
@@ -315,7 +411,7 @@ function execute_operation!(lang::LangC, T, graph, node, dealloc_list, mem)
             alloc_slot!(mem, lhsmem_i, :dummy)
             push_code!(
                 code,
-                "memcpy($lhsmem, $parent1mem, " * "n*n*sizeof(*master_mem));",
+                "memcpy($lhsmem, $parent1mem, " * "n * n * sizeof(*master_mem));",
             )
         end
 
@@ -356,12 +452,13 @@ function execute_operation!(lang::LangC, T, graph, node, dealloc_list, mem)
                 setdiff!(dealloc_list, [parent2])
                 set_slot_number!(mem, get_slot_number(mem, parent2), node)
             else
+                parent2mem = get_slot_name(mem, parent2)
                 (nodemem_i, nodemem) = get_free_slot(mem)
                 alloc_slot!(mem, nodemem_i, node)
                 push_code!(
                     code,
                     "memcpy($nodemem, $parent2mem, " *
-                    "n*n*sizeof(*master_mem));",
+                    "n * n*sizeof(*master_mem));",
                 )
             end
 
@@ -382,140 +479,96 @@ function execute_operation!(lang::LangC, T, graph, node, dealloc_list, mem)
         end
 
     elseif op == :lincomb
-        coeff_vars = Vector{String}(undef, 2)
 
-        # Coefficients of graph should have type T.
-        (coeff1, coeff1_code) = assign_coeff(lang, graph.coeffs[node][1], 1)
-        for s in coeff1_code
-            push_code!(code, s)
-        end
-        (coeff2, coeff2_code) = assign_coeff(lang, graph.coeffs[node][2], 2)
-        for s in coeff2_code
-            push_code!(code, s)
-        end
+        setdiff!(dealloc_list, [:I])
 
-        # Reuse parent in deallocation list for output (saves memcopy).
-        if ((parent1 in dealloc_list) || (parent2 in dealloc_list))
-            recycle_parent = (parent1 in dealloc_list) ? parent1 : parent2
-            push_comment!(code, "Smart lincomb recycle $recycle_parent")
-
-            # Perform operation.
-            nodemem = get_slot_name(mem, recycle_parent)
-            if (recycle_parent == parent1)
-                if p2_is_identity
-                    push_code!(
-                        code,
-                        "cblas_$blas_prefix" *
-                        "scal(n*n, $coeff1, " *
-                        "$nodemem, 1);",
-                    )
-                    push_code!(
-                        code,
-                        "cblas_$blas_prefix" *
-                        "axpby(n, " *
-                        "$coeff2, &ONE, 0,\n" *
-                        "             $rone, $nodemem, n+1);",
-                    )
-                else
-                    push_code!(
-                        code,
-                        "cblas_$blas_prefix" *
-                        "axpby(n*n, " *
-                        "$coeff2, $parent2mem, 1,\n" *
-                        "             $coeff1, $nodemem, 1);",
-                    )
-                end
+        fused_sum = (join("x*" .* string.(graph.parents[node]), " + "))
+        push_comment!(code, "Computing $node = $fused_sum")
+
+        # Set coefficients.
+        coeff_names = Vector()
+        coeff_values = Vector()
+        parent_mems = Vector()
+        id_coefficient = 0
+        counter = 1
+        for (i, v) in enumerate(graph.coeffs[node])
+            n = graph.parents[node][i]
+            if v == 0
+                continue
+            end
+            if (n == :I)
+                # Coefficient of identities.
+                id_coefficient += v
             else
-                if p1_is_identity
-                    push_code!(
-                        code,
-                        "cblas_$blas_prefix" *
-                        "scal(n*n, $coeff2, " *
-                        "$nodemem, 1);",
-                    )
-                    push_code!(
-                        code,
-                        "cblas_$blas_prefix" *
-                        "axpby(n, " *
-                        "$coeff1, &ONE, 0,\n" *
-                        "             $rone, $nodemem, n+1);",
-                    )
-                else
-                    push_code!(
-                        code,
-                        "cblas_$blas_prefix" *
-                        "axpby(n*n, " *
-                        "$coeff1, $parent1mem, 1,\n" *
-                        "             $coeff2, $nodemem, 1);",
-                    )
-                end
+                # Coefficient of other nodes.
+                (coeff_value, coeff_i, coeff_i_code) = declare_coeff(lang, v,
+                                                 "$node" * "_" * "$counter")
+                counter += 1
+                push_code!(code, coeff_i_code)
+                push!(coeff_names, coeff_i)
+                push!(coeff_values, coeff_value)
+                push!(parent_mems, get_slot_name(mem, n))
+            end
+        end
+
+        # Determine a target mem: nodemem
+        nodemem = nothing
+        recycle_parent = nothing
+        for n in graph.parents[node]
+            if (n in dealloc_list) # If it's about to be deallocated
+                nodemem = get_slot_name(mem, n)
+                recycle_parent = n
+                break
             end
+        end
+        if (isnothing(recycle_parent))
+            (nodemem_i, nodemem) = get_free_slot(mem)
+            alloc_slot!(mem, nodemem_i, node)
+        else
+            push_comment!(code, "Smart lincomb recycle $recycle_parent.")
+            # Move parent to be recycled to the beginning of the vector.
+            index = findfirst(x -> occursin(nodemem, x), parent_mems)
+
+            splice!(parent_mems, index)
+            pushfirst!(parent_mems, nodemem)
+
+            coeff_name = coeff_names[index]
+            splice!(coeff_names, index)
+            pushfirst!(coeff_names, coeff_name)
+
+            coeff_value = coeff_values[index]
+            splice!(coeff_values, index)
+            pushfirst!(coeff_values, coeff_value)
+        end
+
+        # Write the linear combination.
+        if isempty(coeff_names)
+            push_code!(code, "memset($nodemem, 0, n * n * sizeof(*master_mem));")
+        else
+            push_code!(code, "for (size_t i = 0; i < n * n; i++) {")
+            add_lincomb_body(code, lang, T, nodemem,
+                             coeff_names, coeff_value, parent_mems)
+            push_code!(code, "}")
+        end
+
+        if id_coefficient != 0
+            (coeff_value, coeff_id, coeff_id_code) = declare_coeff(lang,
+                    id_coefficient, "$node" * "_0")
+            push_code!(code, coeff_id_code)
+            push_code!(code, "for (size_t i = 0; i < n * n; i += n + 1) {")
+            add_lincomb_identity_body(code, lang, T, nodemem, coeff_id, coeff_value)
+            push_code!(code, "}")
+        end
 
-            # Remove output from deallocation list.
+        if !(isnothing(recycle_parent))
+            # Avoid deallocate recycled parent
             setdiff!(dealloc_list, [recycle_parent])
 
-            # Update CodeMem.
+            # Set the memory pointer
             j = get_slot_number(mem, recycle_parent)
             set_slot_number!(mem, j, node)
-        else
-            # Allocate new slot for result.
-            (nodemem_i, nodemem) = get_free_slot(mem)
-            alloc_slot!(mem, nodemem_i, node)
-
-            # Compute linear combination.
-            if p1_is_identity
-                push_code!(
-                    code,
-                    "memcpy($nodemem, $parent2mem, " *
-                    "n*n*sizeof(*master_mem));",
-                )
-                push_code!(
-                    code,
-                    "cblas_$blas_prefix" *
-                    "scal(n*n, $coeff2, " *
-                    "$nodemem, 1);",
-                )
-                push_code!(
-                    code,
-                    "cblas_$blas_prefix" *
-                    "axpby(n, " *
-                    "$coeff1, &ONE, 0,\n" *
-                    "             $rone, $nodemem, n+1);",
-                )
-            elseif p2_is_identity
-                push_code!(
-                    code,
-                    "memcpy($nodemem, $parent1mem, " *
-                    "n*n*sizeof(*master_mem));",
-                )
-                push_code!(
-                    code,
-                    "cblas_$blas_prefix" *
-                    "scal(n*n, $coeff1, " *
-                    "$nodemem, 1);",
-                )
-                push_code!(
-                    code,
-                    "cblas_$blas_prefix" *
-                    "axpby(n*n, " *
-                    "$coeff2, &ONE, 0,\n" *
-                    "             $rone, $nodemem, n+1);",
-                )
-            else
-                push_code!(
-                    code,
-                    "memcpy($nodemem, $parent2mem, " *
-                    "n*n*sizeof(*master_mem));",
-                )
-                push_code!(
-                    code,
-                    "cblas_$blas_prefix" *
-                    "axpby(n*n, " *
-                    "$coeff1, $parent1mem, 1,\n" *
-                    "             $coeff2, $nodemem, 1);",
-                )
-            end
         end
+
     else
         error("Unknown operation")
     end
@@ -531,36 +584,22 @@ function execute_operation!(lang::LangC, T, graph, node, dealloc_list, mem)
     return (code, nodemem)
 end
 
-function scalar_to_string(::LangC, z)
-    return "$z"
-end
-function scalar_to_string(::LangC_OpenBLAS, z::T) where {T<:Complex}
-    return "$(real(z)) + $(imag(z))*I"
+random_string() = "rand() / (1.0 * RAND_MAX)";
+function generate_random_entry(lang::LangC, ::Type{T}) where {T<:Real}
+    return assignment_string(lang, "A[i]", random_string())
 end
-function scalar_to_string(::LangC_MKL, z::T) where {T<:Complex}
-    return "{$(real(z)), $(imag(z))}"
-end
-function print_indented_matrix(lang::LangC, code, A; ind_lvl = 1)
-    (m, n) = size(A)
-    for j = 1:n
-        column_string = ""
-        for i = 1:m
-            column_string *=
-                scalar_to_string(lang, A[i, j]) * (i != m ? ", " : "")
-        end
-        column_string *= (j != n ? "," : "")
-        push_code!(code, column_string, ind_lvl = ind_lvl)
-    end
+function generate_random_entry(lang::LangC, ::Type{T}) where {T<:Complex}
+    return assignment_string(lang, "A[i]", random_string(), random_string())
 end
 
 function compilation_string(::LangC_OpenBLAS, fname)
-    return "gcc -o main_compiled $fname -labials -llapacke"
+    return "gcc -o main_compiled $fname -lblas -llapacke"
 end
 function compilation_string(::LangC_MKL, fname)
     return "gcc -o main_compiled $fname -lmkl_rt"
 end
 
-function gen_main(lang::LangC, T, fname, funname; A = 10::Union{Integer,Matrix})
+function gen_main(lang::LangC, T, fname, funname, graph)
     code = init_code(lang)
     if (lang.gen_main)
         (blas_type, blas_prefix) = get_blas_type(lang, T)
@@ -586,25 +625,17 @@ function gen_main(lang::LangC, T, fname, funname; A = 10::Union{Integer,Matrix})
         push_code!(code, "size_t i;")
 
         # Generate matrix.
-        if isa(A, Matrix)
-            n = LinearAlgebra.checksquare(A)
-            push_code!(code, "size_t n = $n;")
-            push_code!(code, "blas_type A[$(n*n)] = {")
-            print_indented_matrix(lang, code, A, ind_lvl = 2)
-            push_code!(code, "};")
-        else # A is an integer
-            n = A
-            push_code!(code, "size_t n = $n;")
-            push_code!(code, "srand(0);")
-            push_code!(code, "blas_type *A = malloc(n*n*sizeof(*A));")
-            push_code!(code, "for(i=0; i<n*n; i++){")
-            push_code!(code, "A[i] = rand() / (1.0*RAND_MAX);", ind_lvl = 2)
-            push_code!(code, "}")
-        end
+        n = 4
+        push_code!(code, "size_t n = $n;")
+        push_code!(code, "srand(0);")
+        push_code!(code, "blas_type *A = malloc(n * n * sizeof(*A));")
+        push_code!(code, "for(i = 0; i < n * n; i++){")
+        push_code!(code, generate_random_entry(lang, eltype(graph)), ind_lvl = 2)
+        push_code!(code, "}")
 
         # Call polynomial evaluation function.
-        push_code!(code, "blas_type *B = malloc(n*n*sizeof(*A));")
-        push_code!(code, "$blas_prefix$funname(A,n,B);")
+        push_code!(code, "blas_type *B = malloc(n * n * sizeof(*A));")
+        push_code!(code, "$blas_prefix$funname(A, n, B);")
         push_code!(code, "return 0;")
         push_code!(code, "}", ind_lvl = 0)
     end
diff --git a/src/code_gen/gen_code.jl b/src/code_gen/gen_code.jl
index ecc49ad..d0a6ebf 100644
--- a/src/code_gen/gen_code.jl
+++ b/src/code_gen/gen_code.jl
@@ -19,7 +19,7 @@ export gen_code
 #                    dealloc_list,   mem)
 
 # Fallback. Default to no main function
-function gen_main(lang, T, fname, funname)
+function gen_main(lang, T, fname, funname, graph)
     return init_code(lang)
 end
 
@@ -137,7 +137,7 @@ function _gen_code(fname, graph, lang, priohelp, funname, precomputed_nodes)
     println(file, to_string(function_end(lang, graph, mem)))
 
     # Generate main function, if necessary.
-    exec_code = gen_main(lang, T, fname, funname)
+    exec_code = gen_main(lang, T, fname, funname, graph)
     println(file, to_string(exec_code))
 
     if (fname isa String)
diff --git a/src/code_gen/gen_julia_code.jl b/src/code_gen/gen_julia_code.jl
index 50ab7de..70a7887 100644
--- a/src/code_gen/gen_julia_code.jl
+++ b/src/code_gen/gen_julia_code.jl
@@ -200,7 +200,7 @@ function function_init(lang::LangJulia, T, mem, graph, precomputed_nodes)
         push_code!(code,"$thisslotname = "*lang.alloc_function(i));
     end
 
-    push_comment!(code,"Assign precomputed nodes memslots ");
+    push_comment!(code,"Assign precomputed nodes memslots");
 
     precomp_nodes_string = join(repeat("A",jj),",");
     for (i, n) in enumerate(precomputed_nodes)
diff --git a/src/compress_graph.jl b/src/compress_graph.jl
index 78562e4..93c9f3e 100644
--- a/src/compress_graph.jl
+++ b/src/compress_graph.jl
@@ -3,6 +3,7 @@
 
 export compress_graph_dangling!
 export compress_graph_zero_coeff!
+export has_lincomb_with_identity
 export has_identity_lincomb
 export has_trivial_nodes
 export compress_graph_trivial!
@@ -110,24 +111,19 @@ end
 
 # Check if node is trivial because of an identity on the left (:mult or :ldiv).
 function is_trivial_left(graph, node)
-    op = graph.operations[node]
-    if (op == :lincomb)
-        return false;
-    else
-        return graph.parents[node][1] == :I &&
-            (op == :mult || op == :ldiv)
-    end
+    return graph.operations[node] ∈ [:mult, :ldiv] &&
+        graph.parents[node][1] == :I
 end
 
 # Check if node is trivial because of an identity on the right (:mult).
 function is_trivial_right(graph, node)
+    return graph.operations[node] == :mult &&
+        graph.parents[node][2] == :I
+end
 
-    op = graph.operations[node]
-    if (op == :lincomb)
-        return false;
-    else
-        return (graph.parents[node][2] == :I && op == :mult)
-    end
+# Check if node is trivial.
+function is_trivial(graph, node)
+    return is_trivial_left(graph, node) || is_trivial_right(graph, node)
 end
 
 # Return a parent that can replace a node.
@@ -159,6 +155,18 @@ function replace_node!(graph, node, replacement, cref)
     return delete_crefs!(cref, node)
 end
 
+"""
+    has_lincomb_with_identity(graph) -> Bool
+
+Checks whether the graph has a linear combination that includes the identity
+matrix.
+"""
+function has_lincomb_with_identity(graph)
+    lincombs = [key for (key, value) in graph.operations if value == :lincomb]
+    lincomb_parents = [graph.parents[key] for key in lincombs]
+    return any(x -> :I in x, lincomb_parents)
+end
+
 """
     has_identity_lincomb(graph) -> Bool
 
@@ -166,15 +174,9 @@ Checks whether the graph has a node which is a linear combination of identity
 matrices.
 """
 function has_identity_lincomb(graph)
-    has_identity_lincomb = false
-    for (key, parents) in graph.parents
-        if graph.operations[key] == :lincomb &&
-           all(parents .== :I)
-            has_identity_lincomb = true
-            break
-        end
-    end
-    return has_identity_lincomb
+    lincombs = [key for (key, value) in graph.operations if value == :lincomb]
+    lincomb_parents = [graph.parents[key] for key in lincombs]
+    return any(x -> all(x .== :I), lincomb_parents)
 end
 
 """
@@ -184,14 +186,7 @@ Checks whether the graph has trivial nodes, that is, multiplications by the
 identity or linear systems whose coefficient is the identity matrix.
 """
 function has_trivial_nodes(graph)
-    has_trivial_nodes = false
-    for (key, parents) in graph.parents
-        if is_trivial_left(graph, key) || is_trivial_right(graph, key)
-            has_trivial_nodes = true
-            break
-        end
-    end
-    return has_trivial_nodes
+    return any(is_trivial(graph, key) for (key, value) in graph.parents)
 end
 
 """
@@ -362,7 +357,8 @@ end
 """
     compress_graph_passthrough!(graph,cref=[];verbose=false);
 
-Identifies lincombs lincomb of length one with coeff equal to one. The node has no effect. It redirect appropriately.
+Identifies lincombs of length one with coeff equal to one.
+The node has no effect and is redirected appropriately.
 """
 
 function compress_graph_passthrough!(graph,
diff --git a/test/code_gen_test.jl b/test/code_gen_test.jl
index 8ff8085..fa2d52f 100644
--- a/test/code_gen_test.jl
+++ b/test/code_gen_test.jl
@@ -11,8 +11,11 @@ using LinearAlgebra, StaticArrays
     # Test julia code generation
     for i = 1:length(a)
         (graph, crefs) = graph_ps_degopt([3 4 2 -1 2 a[i]])
-        add_lincomb!(graph,:Q,[2.0],graph.outputs) # Check
-        add_ldiv!(graph, :R0, :B4, :Q)
+        add_lincomb!(graph, :Q, [2.0], graph.outputs)# Check
+        add_ldiv!(graph, :R2, :B4, :Q)
+        add_lincomb!(graph, :Is, [2., 3., 5.], [:I, :I, :I])
+        add_lincomb!(graph, :R1, [2., 2.], [:Is, :R2])
+        add_ldiv!(graph, :R0, :R1, :I)
         clear_outputs!(graph)
         add_output!(graph, :R0)
 
@@ -70,9 +73,13 @@ using LinearAlgebra, StaticArrays
 
     for i = 1:4 #Not high-precision test for Matlab and C
         (graph, crefs) = graph_ps_degopt([3 4 2 a[i] 1 0])
-        add_ldiv!(graph, :R0, :B4, graph.outputs[1])
+        add_ldiv!(graph, :R1, :B4, graph.outputs[1])
+        add_ldiv!(graph, :R2, :R1, :I)
+        add_ldiv!(graph, :R3, :B4, :A)
+        add_lincomb!(graph,:R4, [1., 2., 3.], [:I, :I, :I])
+        add_lincomb!(graph,:R5, [2., 3., 4., 5., 6.], [:I, :R1, :R2, :R2, :R4])
         clear_outputs!(graph)
-        add_output!(graph, :R0)
+        add_output!(graph, :R5)
 
         # Test Matlab code generation (not execution)
         fname = tempname() * ".m"
@@ -80,20 +87,16 @@ using LinearAlgebra, StaticArrays
         rm(fname)
 
         # Test C code generation
-        fname = tempname() * ".c"
-        gen_code(fname, graph, lang = LangC_MKL())
-        rm(fname)
-
-        fname = tempname() * ".c"
-        gen_code(fname, graph, lang = LangC_MKL(true))
-        rm(fname)
-
-        fname = tempname() * ".c"
-        gen_code(fname, graph, lang = LangC_OpenBLAS())
-        rm(fname)
+        for gen_main in (true, false)
+            for overwrite_input in (true, false)
+                fname = tempname() * ".c"
+                gen_code(fname, graph, lang = LangC_MKL(gen_main, overwrite_input))
+                rm(fname)
 
-        fname = tempname() * ".c"
-        gen_code(fname, graph, lang = LangC_OpenBLAS(true))
-        rm(fname)
+                fname = tempname() * ".c"
+                gen_code(fname, graph, lang = LangC_OpenBLAS(gen_main, overwrite_input))
+                rm(fname)
+            end
+        end
     end
 end
diff --git a/test/compress.jl b/test/compress.jl
index d16575e..e326003 100644
--- a/test/compress.jl
+++ b/test/compress.jl
@@ -34,13 +34,25 @@ using LinearAlgebra
     # Check that we actually removed something
     @test size(get_all_cref(graph2), 1) < size(get_all_cref(graph1), 1)
 
+    ## Test function to detect if graph has linear combination with identities
     ## Test function to detect if graph has linear combination of identities
     graph = Compgraph()
     add_lincomb!(graph, :P1, 1.0, :A, 1.0, :I)
     add_lincomb!(graph, :P2, 1.0, :I, 1.0, :I)
     add_mult!(graph, :P0, :P1, :P2)
+    @test has_lincomb_with_identity(graph) == true
     @test has_identity_lincomb(graph) == true
 
+    graph = Compgraph()
+    add_lincomb!(graph, :P, 1.0, :A, 1.0, :I)
+    @test has_lincomb_with_identity(graph) == true
+    @test has_identity_lincomb(graph) == false
+
+    graph = Compgraph()
+    add_mult!(graph, :P, :A, :A)
+    @test has_lincomb_with_identity(graph) == false
+    @test has_identity_lincomb(graph) == false
+
     ## Test trivial node detection and removal
     graph = Compgraph()
     add_mult!(graph, :AI, :A, :I)
@@ -58,10 +70,15 @@ using LinearAlgebra
     @test has_trivial_nodes(graph) == true
 
     graph = Compgraph()
-    add_mult!(graph, :I2, :I, :I) # This node is trivial
-    add_mult!(graph, :AI, :A, :I2) # This node is trivial
+    add_lincomb!(graph, :A2, [2], [:A])
+    add_output!(graph, :A2)
+    @test has_trivial_nodes(graph) == false
+
+    graph = Compgraph()
+    add_mult!(graph, :I2, :I, :I)   # This node is trivial
+    add_mult!(graph, :AI, :A, :I2)  # This node is trivial
     add_mult!(graph, :IA, :I2, :AI) # This node is trivial
-    add_ldiv!(graph, :P0, :I, :IA) # This node is trivial
+    add_ldiv!(graph, :P0, :I, :IA)  # This node is trivial
     add_output!(graph, :P0)
     graph1 = graph
     graph2 = deepcopy(graph)