Generate code with precomputed nodes

[jarlebring]

Some implementations of the matrix exponential does a decision of which computation should be done based on A^2, A^4 (and A^8 ?). This does not fit perfectly with the way our generated code, since the computation of these powers are integrated into the generated code.

Proposed approach: Generated code can be generated without the computation of specific nodes, which are assumed to be precomputed.
For example

gen_code(graph,precomputed_nodes=[:A2,:A4])

would generate computation where the operations to get :A2 and :A4 are treated differently. Maybe effectively they are also inputs? cf discussion in #9

More details to be specified.

[jarlebring]

Question for @mfasi: In the generated c-code, all memory is allocated with one malloc-call. In the julia code they are generated as a vector of matrices. Would there be a disadvantage in the C-code if the matrices are stored as an array of pointers to each matrix rather than an array with all matrices? cf

    (double*) memslots[MAXSLOTS];
    for (i=0;i<MAXSLOTS;i++){
       memslots[i]=malloc(n*n*sizeof(**memslots));
    }

I think "array of arrays" could make the precomputation feature easier since one could just assign the memslot pointers to the precomputed matrices. This would be difficult to achieve (without an additional copy) with if we have just one big array.

In fact, atm the C-code uses one more memory slot in comparison to the julia code because the julia code assigns memslot[2]=A and later recycles it (i.e. overwrites A), e.g., for m6_SID_2_22. I think this would also be fixable with an array of arrays.

[mfasi]

Would there be a disadvantage in the C-code if the matrices are stored as an array of pointers to each matrix rather than an array with all matrices?

In general, a single malloc is faster than many small ones if the system has free RAM to spare, which is generally the case currently. Things would probably be different on a system with high memory pressure, where trying to allocate a larger chunk of memory is likely to be slower and could in principle fail. But as far as I can tell the difference should be minimal, thus it's probably a good idea to switch to an array of pointers if that simplifies the implementation and reduces the memory footprint.

[jarlebring]

Okej. If the single vs multiple malloc-calls is the issue, we can still use arrays of arrays so we can get the best of both worlds.
Something like:

main_memory=malloc(n*n*sizeof(*main_memory)*memsolts)
for (i=1;i<(MAXSLOTS+1);i++){
       memslots[i]=main_memory+n*n*sizeof(*main_memory)*i;
}

and set the zero memslots to A:

memslots[0]=A

(Memslot counter maybe has to be done a bit more carefully see Julia generated code.)

[jarlebring]

The c-generator now has the above memslots functionality implemented in 21cc3fe and 6b03446. It generates code starting like this

void ddummy(double *A, const size_t n, double *output) {
    size_t max_memslots = 7;

    /* Initializations. */
    double coeff1, coeff2;
    const double ZERO = 0.0;
    const double ONE = 1.0;
    double *master_mem = malloc(n*n*max_memslots*sizeof(*master_mem));
    int j;
    double *memslots[8]; // One more than max_memslots
    memslots[0] = A; // First slot is A (and can be recycled)
    /* The other slots are pointers to allocated memory */
    for (j=0; j<max_memslots; j++) memslots[j+1] = master_mem+j*n*n;
/*     Computation order: Bb2 Ba4_2 Bb4_2 Ba2 B2 Ba4_3 Bb4_3 Ba3_2 Ba3 Ba5_2 Ba5_3 Bb5_2 Bb5_3 T2k2 T2k3 Bb3_2 Bb3 B3 Ba4 Bb4 B4 Ba5_4 Ba5 T2k4 Bb5_4 T2k5 Bb5 B5 T2k7 */

    /* Computing Bb2 with operation: lincomb */
    coeff1 = 0.0;
    coeff2 = 1.0;
    memcpy(memslots[2-1], A, n*n*sizeof(*master_mem));
    cblas_dscal(n*n, coeff2, memslots[2-1], 1);
    cblas_daxpby(n, coeff1, &ONE, 0,

I had to change const double *A to double *A in the arg to avoid a warning. I tried it on some small examples. Seems to work okay. Could you have a brief look @mfasi ?

[jarlebring]

Actually, this needs a careful second look since I'm uncertain if what I did interferes with this code

function langc_get_slot_name(mem,key)
    return key==:A ? "A" : get_slot_name(mem,key)

[mfasi]

The change makes sense to me. As long as A is not referenced after being deallocated, which should be the case, I don't think that there is a problem with the current version of langc_get_slot_name(). I do believe that the code can be simplified to

langc_get_slot_name(mem,key)=get_slot_name(mem,key)

in the current setting.

One minor issue is that now the code generation fails for the following degenerate graphs:

graph1=CompGraph()
add_output!(graph1,:I)

graph2=CompGraph()
add_output!(graph2,:A)

For LangJulia, the generation also fails for the following

graph=CompGraph()
add_lincomb!(graph,:I2,2.1,:I,2.2,:I)
add_output!(graph,:I2)

All degnerate cases are handled correctly for LangMatlab.

[jarlebring]

Concerning minor issue: I see the code generation functionality as a tool to help construct efficient matrix function algorithm, with some expectation on the programmer to do manual post processing. For that reason, I think we don't have to think so much about corner cases such which do not correspond to efficient graphs. I don't see a situation where one would want to start by forming a lincomb of two identites. The code generation is already difficult code.

[mfasi]

Concerning minor issue: I see the code generation functionality as a tool to help construct efficient matrix function algorithm, with some expectation on the programmer to do manual post processing. For that reason, I think we don't have to think so much about corner cases such which do not correspond to efficient graphs. I don't see a situation where one would want to start by forming a lincomb of two identites. The code generation is already difficult code.

Good point. I added a check in code_gen that returns an error message for those trivial cases. I think it's the least involved solution.

[jarlebring]

More thoughts on how precomputed nodes can be handled.

In a generated code we can add several parameters

function dummy(A,A2,A4)

In the gen_code function we add a kwarg which is precomputed_nodes=[:A]. In the two main loops in the code generation
https://github.com/jarlebring/GraphMatFun.jl/blob/1a32bf51933996ed2333113f589210dee6bd5cc8/src/code_gen/gen_code.jl#L101-L110
we can add a check and basically not add any operations to compute the precomputed nodes:

if (node in precomputed_nodes) 
    (exec_code, ..) =execute_..