R side

I run use_r("cpp11gaussjordan-package") and added the following contents.

#' @useDynLib cpp11gaussjordan, .registration = TRUE
NULL

#' Invert (some) square matrices
#' @export
#' @param A numeric matrix
#' @return numeric matrix
#' @examples
#' A <- matrix(c(2,1,3,-1), nrow = 2, ncol = 2)
#' invert_matrix(A)
invert_matrix <- function(A) {
  invert_matrix_(A)
}

#' Solve (some) linear systems
#' @export
#' @param A numeric matrix
#' @param b numeric matrix
#' @return numeric matrix
#' @examples
#' A <- matrix(c(2,1,3,-1), nrow = 2, ncol = 2)
#' b <- matrix(c(7,4), nrow = 2, ncol = 1)
solve_system <- function(A,b) {
  solve_system_(A, b)
}

Here solve_system() is the “end-user” function, and solve_system_() is the “development” function. The “development” functions were written in C++ and cpp11 handles adding the .Call() when required to link the C++ code and be able to use the C++ functions from R.

C++ side

I replaced code.cpp contents with the following lines.

#include <cpp11.hpp>
#include <cpp11/doubles.hpp>

using namespace cpp11;

[[cpp11::register]] doubles_matrix<> invert_matrix_(doubles_matrix<> A)
{
    // Obtain (X'X)^-1 via Gauss-Jordan

    // Check dimensions

    int N = A.nrow();
    int M = A.ncol();

    if (N != M)
    {
        stop("X must be a square matrix");
    }

    // Copy the matrix

    writable::doubles_matrix<> Acopy(N, N);

    for (int i = 0; i < N; i++)
    {        
        for (int j = 0; j < N; j++)
        {
            Acopy(i, j) = A(i, j);
        }
    }

    // Create the identity matrix as a starting point for Gauss-Jordan

    writable::doubles_matrix<> Ainv(N, N);

    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < N; j++)
        {
            if (i == j)
            {
                Ainv(i, j) = 1.0;
            }
            else
            {
                Ainv(i, j) = 0.0;
            }
        }
    }

    // Overwrite Ainv by steps with the inverse of A
    // in other words, find the echelon form of A

    for (int i = 0; i < M; i++)
    {
        double a = Acopy(i, i);

        for (int j = 0; j < M; j++)
        {
            Acopy(i, j) /= a;
            Ainv(i, j) /= a;
        }

        for (int j = 0; j < M; j++)
        {
            if (i != j)
            {
                a = Acopy(j, i);

                for (int k = 0; k < M; k++)
                {
                    Acopy(j, k) -= Acopy(i, k) * a;
                    Ainv(j, k) -= Ainv(i, k) * a;
                }
            }
        }
    }

    return Ainv;
}

[[cpp11::register]] doubles_matrix<> solve_system_(doubles_matrix<> A, doubles_matrix<> b)
{
    // Use Gauss-Jordan to solve the system Ax = b

    // Check dimensions

    int N1 = A.nrow();
    int M1 = A.ncol();

    if (N1 != M1)
    {
        stop("A must be a square matrix");
    }

    int N2 = b.nrow();
    int M2 = b.ncol();

    if (N1 != N2)
    {
        stop("b must have the same number of rows as A");
    }

    if (M2 != 1)
    {
        stop("b must be a column vector");
    }

    // Obtain the inverse

    // writable::doubles_matrix<> Acopy(N1,N1);

    // for (int i = 0; i < N1; i++)
    // {        
    //     for (int j = 0; j < N1; j++)
    //     {
    //         Acopy(i, j) = A(i, j);
    //     }
    // }

    // doubles_matrix<> Ainv = invert_matrix_(Acopy);

    // I don´t need to create a copy in this case
    
    doubles_matrix<> Ainv = invert_matrix_(A);

    // Multiply Ainv by b

    writable::doubles_matrix<> x(N1, 1);

    for (int i = 0; i < N1; i++)
    {
        x(i, 0) = 0.0;

        for (int j = 0; j < N1; j++)
        {
            x(i, 0) += Ainv(i, j) * b(j, 0);
        }
    }

    return x;
}

I also needed a Makevars file to indicate the compiler (either clang or gcc) to use the vendored headers. The content is the following.

PKG_CPPFLAGS = -I../inst/include

Because of the vendoring option, I also had to remove the LinkingTo: cpp11 line from DESCRIPTION.

Putting all together

I created a vscode-install.r file in the root of the project, and I added it to .Rbuildignore. It contains the following lines.

# cpp_vendor() # run only when updating C++ headers
clean_dll()
cpp_register()
document()
install()
# load_all()

To test that the functions work, after running install() I solved a simple system.

> library(cpp11gaussjordan)
> A <- matrix(c(2,1,3,-1), nrow = 2, ncol = 2)
> b <- matrix(c(7,4), nrow = 2, ncol = 1)
> solve_system(A,b)
     [,1]
[1,]  3.8
[2,] -0.2
> solve(A) %*% b
     [,1]
[1,]  3.8
[2,] -0.2