Functionality extension guidelines#

This document provides information about the different ways to implement user-specific routines, solvers, or preconditioners to the rocALUTION library package. Additional features can be added in multiple ways. Additional solver and preconditioner functionality that uses the existing backend functionality performs well on accelerator devices without the need for expert GPU programming knowledge. Also, those not interested in using accelerators are not required to perform HIP and GPU-related programming tasks to add additional functionality.

In the following sections, different levels of functionality enhancements are illustrated. These examples can be used as guidelines to extend rocALUTION step by step with your own routines. Please note, that user-added routines can also be added to the main GitHub repository using pull requests.

`LocalMatrix` functionality extension#

This section demonstrates how to extend the LocalMatrix class with an additional routine. The routine supports both Host and Accelerator backend. Furthermore, the routine requires the matrix to be in CSR format. Here are the steps to extend the LocalMatrix functionality:

1. API enhancement#

To make the new routine available through the API, modify the LocalMatrix class. The corresponding header file local_matrix.hpp is located in src/base/. The new routines can be added as public member functions as shown below:

...
void ConvertTo(unsigned int matrix_format, int blockdim);

void MyNewFunctionality(void);

virtual void Apply(const LocalVector<ValueType>& in, LocalVector<ValueType>* out) const;
virtual void ApplyAdd(const LocalVector<ValueType>& in,
...

For the implementation of the new API function, it is important to know the location of the availability of this functionality. To add support for any backend and matrix format, format conversions are required if MyNewFunctionality() is only supported for CSR metrices. This is subject to the API function implementation:

template <typename ValueType>
void LocalMatrix<ValueType>::MyNewFunctionality(void)
{
    // Debug logging
    log_debug(this, "LocalMatrix::MyNewFunctionality()");

#ifdef DEBUG_MODE
    // If we are in debug mode, perform an additional matrix sanity check
    this->Check();
#endif

    // If no non-zero entries, do nothing
    if(this->GetNnz() > 0)
    {
        // As we want to implement our function only for CSR, we first need to convert
        // the matrix to CSR format
        unsigned int format = this->GetFormat();
        int blockdim = this->GetBlockDimension();
        this->ConvertToCSR();

        // Call the corresponding base matrix implementation
        bool err = this->matrix_->MyNewFunctionality();

        // Check its return type
        if((err == false) && (this->is_host_() == true))
        {
            // If our matrix is on the host, the function call failed.
            LOG_INFO("Computation of LocalMatrix::MyNewFunctionality() failed");
            this->Info();
            FATAL_ERROR(__FILE__, __LINE__);
        }

        // Run backup algorithm on host, in case the accelerator version failed
        if(err == false)
        {
            // Move matrix to host
            bool is_accel = this->is_accel_();
            this->MoveToHost();

            // Try again
            if(this->matrix_->MyNewFunctionality() == false)
            {
                LOG_INFO("Computation of LocalMatrix::MyNewFunctionality() failed");
                this->Info();
                FATAL_ERROR(__FILE__, __LINE__);
            }

            // On a successful host call, move the data back to the accelerator
            // if initial data was on the accelerator
            if(is_accel == true)
            {
                // Print a warning, that the algorithm was performed on the host
                // even though the initial data was on the device
                LOG_VERBOSE_INFO(2, "*** warning: LocalMatrix::MyNewFunctionality() was performed on the host");

                this->MoveToAccelerator();
            }
        }

        // Convert the matrix back to CSR format
        if(format != CSR)
        {
            // Print a warning, that the algorithm was performed in CSR format
            // even though the initial matrix format was different
            LOG_VERBOSE_INFO(2, "*** warning: LocalMatrix::MyNewFunctionality() was performed in CSR format");

            this->ConvertTo(format, blockdim);
        }
    }

#ifdef DEBUG_MODE
    // Perform additional sanity check in debug mode, because this is a non-const function
    this->Check();
#endif
}

Similarly, you can implement host-only functions. In this case, initial data explicitly needs to be moved to the host backend using the API implementation.

The next step is to implement the actual functionality in the BaseMatrix class.

2. Enhancement of the `BaseMatrix` class#

To make the new routine available in the base class, first modify the BaseMatrix class. The corresponding header file base_matrix.hpp is located in src/base/. The new routines can be added as public member functions, e.g.

...
virtual bool ILU0Factorize(void);

/// Perform MyNewFunctionality algorithm
virtual bool MyNewFunctionality(void);

/// Perform LU factorization
...

We don’t implement the purely virtual MyNewFunctionality() as we don’t supply an implementation for all base classes. We decided to implement it only for CSR format and hence need to return an error flag, so that the LocalMatrix class is aware of the failure and can convert it to CSR.

template <typename ValueType>
bool MyNewFunctionality(void)
{
    return false;
}

3. Platform-specific host implementation#

To satisfy the rocALUTION host backup philosophy, there must be a host implementation available. Hence, for the new function to succeed, there must be backend implementation available. Place the host implementation in src/base/host/host_matrix_csr.cpp as we decided to make it available for CSR format.

...
virtual bool ILUTFactorize(double t, int maxrow);

virtual bool MyNewFunctionality(void);

virtual void LUAnalyse(void);
...

template <typename ValueType>
bool HostMatrixCSR<ValueType>::MyNewFunctionality(void)
{
    // Place some asserts to verify sanity of input data

    // Our algorithm works only for squared metrices
    assert(this->nrow_ == this->ncol_);
    assert(this->nnz_ > 0);

    // place the actual host based algorithm here:
    // for illustration, we scale the matrix by its inverse diagonal
    for(int i = 0; i < this->nrow_; ++i)
    {
        int row_begin = this->mat_.row_offset[i];
        int row_end   = this->mat_.row_offset[i + 1];

        bool diag_found = false;
        ValueType inv_diag;

        // Find the diagonal entry
        for(int j = row_begin; j < row_end; ++j)
        {
            if(this->mat_.col[j] == i)
            {
                diag_found = true;
                inv_diag = static_cast<ValueType>(1) / this->mat_.val[j];
            }
        }

        // Our algorithm works only with full rank
        assert(diag_found == true);

        // Scale the row
        for(int j = row_begin; j < row_end; ++j)
        {
            this->mat_.val[j] *= inv_diag;
        }
    }

    return true;
}

4. Platform-specific HIP implementation#

You can now add an additional implementation for the HIP backend using HIP programming framework. This is required to make your algorithm available on accelerators so that rocALUTION doesn’t need to switch to the host backend on function calls anymore. Add the HIP implementation src/base/hip/hip_matrix_csr.cpp in this case.

...
virtual bool ILU0Factorize(void);

virtual bool MyNewFunctionality(void);

virtual bool ICFactorize(BaseVector<ValueType>* inv_diag = NULL);
...

template <typename ValueType>
bool HIPAcceleratorMatrixCSR<ValueType>::MyNewFunctionality(void)
{
    // Place some asserts to verify sanity of input data

    // Our algorithm works only for squared metrices
    assert(this->nrow_ == this->ncol_);
    assert(this->nnz_ > 0);

    // Enqueue the HIP kernel
    hipLaunchKernelGGL((kernel_csr_mynewfunctionality),
                       dim3((this->nrow_ - 1) / this->local_backend_.HIP_block_size + 1),
                       dim3(this->local_backend_.HIP_block_size),
                       0,
                       0,
                       this->mat_.row_offset,
                       this->mat_.col,
                       this->mat_.val);

    // Check for HIP execution error before successfully returning
    CHECK_HIP_ERROR(__FILE__, __LINE__);

    return true;
}

Place the corresponding HIP kernel in src/base/hip/hip_kernels_csr.hpp.

Adding a solver#

This section demonstrates how to add a new solver to rocALUTION. Here are the steps:

Define the API for the new solver

As an example, we add a new IterativeLinearSolver. To achieve this, we use CG as a template. Thus, we first copy src/solvers/krylov/cg.hpp to src/solvers/krylov/mysolver.hpp and src/solvers/krylov.cg.cpp to src/solvers/krylov/mysolver.cpp (assuming we add a krylov subspace solvers).

Modify the cg.hpp and cg.cpp as per your requirement (e.g. change the solver name from CG to MySolver)

Implement each of the following virtual functions present in the class. Follow the implementation details given below:

MySolver(): The constructor of the new solver class.
~MySolver(): The destructor of the new solver class. It calls the Clear() function.
void Print(void) const: Prints some informations about the solver.
void Build(void): Creates all required structures of the solver, e.g. allocates memory and sets the backend of temporary objects.
void BuildMoveToAcceleratorAsync(void): Moves all solver-related objects asynchronously to the accelerator device.
void Sync(void): Synchronizes all solver related objects.
void ReBuildNumeric(void): Rebuilds the solver only numerically.
void Clear(void): Cleans up all solver-relevant structures that have been created using Build().
void SolveNonPrecond_(const VectorType& rhs, VectorType* x): Performs the solving phase Ax=y without the use of a preconditioner.
void SolvePrecond_(const VectorType& rhs, VectorType* x): Performs the solving phase Ax=y with the use of a preconditioner.
void PrintStart_(void) const: Protected function. Called when the solver starts.
void PrintEnd_(void) const: Protected function. Called when the solver ends.
void MoveToHostLocalData_(void): Protected function. Moves all local solver objects to the host.
void MoveToAcceleratorLocalData_(void): Protected function. Moves all local solver objects to the accelerator.

You can also introduce any additional solver-specific member functions.

Make the new solver visible

To make the new solver visible, add it to the src/rocalution.hpp header:

...
#include "solvers/krylov/cg.hpp"
#include "solvers/krylov/mysolver.hpp"
#include "solvers/krylov/cr.hpp"
...

Add the new solver to the CMake compilation list

The CMake compilation list is found in src/solvers/CMakeLists.txt:

...
set(SOLVERS_SOURCES
  solvers/krylov/cg.cpp
  solvers/krylov/mysolver.cpp
  solvers/krylov/fcg.cpp
...

Functionality extension guidelines

Contents