Sparse Linear Algebra (native Kratos implementation)

WORK IN PROGRESS

The Kratos prvides a basic implementation of sparse linear algebra. The implementation is designed to ease FEM operations and to provide a similar interface both in the SMP and MPI case.

Construction and assembly of CSR Matrices (system matrices)

The library provides basic capabilities for the construction of CSR matrices by FEM assembly. A classical challenge for the construction of CSR matrices is that the graph of the matrix needs to be provided at once. The kratos implementation provides a "graph" object designed to allow the sparse insertion of items. CSR matrices can be constructing by providing such a graph.

SMP case

Graph objects

in the SMP case, two "graph-type" objects are implemented with slightly different properties:

• SparseGraph which allows constructing a graph for which neither the number of rows nor the number of columns is known
• SparseContiguousRowGraph which allows constructing a graph for which the number of rows is fixed.

The SparseContiguousRowGraph is threadsafe, while the SparseGraph is not (although diffeent graphs can be constructed in parallel and then merged)

The graph is done by a loop of the type

    SparseGraph<> Agraph; //alternaively here we could use the SparseContiguousRowGraph but we would need to specify the number of rows in the constructo
    for(const auto& c : connectivities)
        Agraph.AddEntries(c); //assuming that c=[1,2], this would fill in the graph all the connectivities (1,1),(1,2),(2,1),(2,2) 
    Agraph.Finalize(); //the graph is considered complete and should not any longer by touched after this call

The graphs offer 3 type of "AddMethods"

• AddEntry(I,J) which adds the single entry I,J
• AddEntries(list_of_ids) which adds all the entries in the positions list_of_ids(i),list_of_ids(j) for any i and j within the size of list_of_ids
• AddEntries(row_ids,column_ids) similar to the previous but for rectangular matrices

A method Has(I,J) allows querying if a given I,J is present in the graph

The method Size() returns the number of rows in the graph.

Graphs can be iterated to obtain the "I","J" positions of the nonzeros in c++ by a code of the type

    for(auto it=rAgraph.begin(); it!=rAgraph.end(); ++it)
    {
        const auto I = it.GetRowIndex();
        for(auto J : *it )
        {
            //do something
        }
    }

A "csr-type" representation of the graphs can be obtained by the function ExportCSRArrays

A very compact "single vector" representation of the graph into a single array can be obtained by the function ExportSingleVectorRepresentation. Such representation can also be used for a fast reconstruction of the graph or for serialization purposes

CSR matrix A native implementation of a CSR (Compressed Sparse Row) matrix, optimized to ease finite element assembly is also available. CsrMatrix objects are **threadsafe" and are designed to be constructed from a Graph, and to be built by fem assemble. An example of construction is as follows:

    CsrMatrix<double> A(Agraph); //construct from an existing graph

    A.BeginAssemble(); //this call need to be performed at the beginning of the assembly phase
    for(const auto& equation_ids : connectivities) //this loop CAN happen in parallel
    {   
        A.Assemble(local_matrix,equation_ids ); //this assembles a SQUARE local_matrix into the positions contained in equation_ids 
    }
    A.FinalizeAssemble(); //this call need to be performed at the end of the assembly phase

The class implements the 3 entry points:

• Assemble(value,I,J) assembling the single scalar "value" into I,J (I,J needs to be in the graph)
• Assemble(square_matrix,positions) "standard" FEM assembly of a square_matrix
• Assemble(rectangual_matrix,row_ids,col_ids) assembling of a rectangular matrix

Useful methods are:

• size1() number of rows
• size2() number of columns
• nnz() number of nonzeros
• index1_data() gives access to the underlying vector of row pointers
• index2_data() gives access to the underlying vector of column indices
• value_data() gives access to the underlying vector of values
• SpMV(x,y) implements y += A*x
• SpMV(alpha,x,beta,y) implements y = alphay + betaA*x
• TransposeSpMV(x,y) implements y += A^T*x
• TransposeSpMV(alpha,x,beta,y) implements y = alphay + betaA^T*x
• ToMap() returns a dictionary of type { (I,J) : value }

MPI case

Construction and assembly of Vectors (system vectors)

While the construction of Vectors ready for FEM assembly provides no difficulty in the SMP case, the effective use in a MPI context requires storing the communication patterns. The exact