LUFactorization

#include <vtkMath.h>

#include <iostream>

namespace {
template <class TReal> TReal** create_matrix(long nrow, long ncol)
{
  typedef TReal* TRealPointer;
  TReal** m = new TRealPointer[nrow];

  TReal* block = static_cast<TReal*>(calloc(nrow * ncol, sizeof(TReal)));
  m[0] = block;
  for (int row = 1; row < nrow; ++row)
  {
    m[row] = &block[row * ncol];
  }
  return m;
}

// Free a TReal matrix allocated with create_matrix().
template <class TReal> void free_matrix(TReal** m)
{
  free(m[0]);
  delete[] m;
}

void OutputMatrix(double** a)
{
  std::cout << "[ " << a[0][0] << " " << a[0][1] << std::endl;
  std::cout << "  " << a[1][0] << " " << a[1][1] << " ]" << std::endl;
}
} // namespace

int main(int, char*[])
{
  // Create and populate the matrix.
  int n = 2;
  double** a = create_matrix<double>(n, n);
  a[0][0] = 4;
  a[0][1] = 3;
  a[1][0] = 6;
  a[1][1] = 3;
  // Here we interchange the rows.
  // a[1][0] = 4;
  // a[1][1] = 3;
  // a[0][0] = 6;
  // a[0][1] = 3;

  /*
  Note:
   [4 3; 6 3] should decompose to [1 0; 3/2 1] * [4 3; 0 -3/2]
   [6 3; 4 3] should decompose to [1 0; 2/3 1] * [6 3; 0 1]
   LUFactorLinearSystem does the following:
    [4 3; 6 3] decomposes to [1 0; 2/3 1] * [6 3; 0 1]
    [6 3; 4 3] decomposes to [1 0; 3/2 1] * [4 3; 0 -3/2]
  */

  std::cout << "a" << std::endl;
  OutputMatrix(a);

  int pivotIndices[2] = {0, 0};

  // Decompose matrix A into the LU form.
  vtkMath::LUFactorLinearSystem(a, pivotIndices, n);

  std::cout << "A decomposed into (unit lower triangular) L and U:"
            << std::endl;
  OutputMatrix(a);
  std::cout << "pivotIndices = (" << pivotIndices[0] << ", " << pivotIndices[1]
            << ")" << endl;

  /*
  The resulting matrix: [6 3; 0.666667 1] is a superposition of L and U,
   with L being a unit lower triangular matrix [[1 0; 0.666667 1]]
   and U being the upper triangular matrix [6 3; 0 1].

  In other words, the diagonal of the resulting matrix A is the diagonal of U.
  The upper right triangle of A is the upper right triangle of U.
  The lower left triangle of A is the lower left triangle of L
  (and remember, the diagonal of L is all 1's).

  To show that the resulting interpretation of the output matrix A is correct,
  we form the matrices following the description above and show that they
  multiply to the original A matrix:
    [1 0; 0.666667 1] * [6 3; 0 1] = [6 3; 4 3]
  We would need to interchange the rows to get [4 3; 6 3]
  */

  free_matrix(a);
  return EXIT_SUCCESS;
}

cmake_minimum_required(VERSION 3.12 FATAL_ERROR)

project(LUFactorization)

find_package(VTK COMPONENTS 
  CommonCore
)

if (NOT VTK_FOUND)
  message(FATAL_ERROR "LUFactorization: Unable to find the VTK build folder.")
endif()

# Prevent a "command line is too long" failure in Windows.
set(CMAKE_NINJA_FORCE_RESPONSE_FILE "ON" CACHE BOOL "Force Ninja to use response files.")
add_executable(LUFactorization MACOSX_BUNDLE LUFactorization.cxx )
  target_link_libraries(LUFactorization PRIVATE ${VTK_LIBRARIES}
)
# vtk_module_autoinit is needed
vtk_module_autoinit(
  TARGETS LUFactorization
  MODULES ${VTK_LIBRARIES}
)

cd LUFactorization/build

cmake ..

cmake -DVTK_DIR:PATH=/home/me/vtk_build ..

make

./LUFactorization