README.md

eliminate.cpp

#ifdef __APPLE__
  extern "C" {
    #include "vecLib/clapack.h"
    #include "vecLib/cblas.h"
  }
#else
  #include <cblas.h>
  #include <lapacke.h>
#endif

/*

  Do the partial elimination in equation system by calculating Schur complement.
  
  IMPORTANT: input matrix is columnwise ordered (Fortran-like) not in C (row-wise) fashion.

  +---+---+        +---+             +----+---+        +---+
  | A | B |}m      | E |}n           |Diag| B'|}m      | E'|}n
  +---+---+  * x = +---+       =>    +----+---+  * x = +---+
  | C | D |}k      | F |}k           | 0  | D'|}k      | F'|}k
  +---+---+        +---+             +----+---+        +---+
    ^   ^                             ^   ^
    m   k                             m   k

*/

void eliminate(double *matrix, double *rhs, int n, int m) {
  int info;
  int ipiv[m];
  int k = n - m;
  char LapackNoTrans = 'N';

  double *A = matrix, 
         *B = matrix + m * n, 
         *C = matrix + m,
         *D = matrix + n * m + m, 
         *E = rhs,
         *F = rhs + m;

  // 1: LU factorize A
  dgetrf_(
      /* M    */ &m, // size
      /* N    */ &m,
      /* A    */ A,    // pointer to data
      /* LDA  */ &n,   // LDA = matrix_size
      /* IPIV */ ipiv, // pivot vector
      /* INFO */ &info);

  if (info != 0) {
    throw lapack_exception("DGETRF failed: LU decomposition not possible");
  }

  // 2: B = A^-1 * B
  dgetrs_(
      /* TRANS */ &LapackNoTrans,
      /* N     */ &m,
      /* NRHS  */ &k,
      /* A     */ A,
      /* LDA   */ &n,
      /* IPIV  */ ipiv,
      /* B     */ B,
      /* LDB   */ &n,
      /* INFO  */ &info);

  if (info != 0) {
    throw lapack_exception("DGETRS failed");
  }

  // 3: E = A^-1 * E
  int one = 1;
  dgetrs_(
      /* TRANS */ &LapackNoTrans,
      /* N     */ &m,
      /* NRHS  */ &one,
      /* A     */ A,
      /* LDA   */ &n,
      /* IPIV  */ ipiv,
      /* B     */ E,
      /* LDB   */ &n,
      /* INFO  */ &info);

  if (info != 0) {
    throw lapack_exception("DGETRS failed");
  }

  // 4: D = D - C * B
  cblas_dgemm(CblasColMajor,
              /* TRANSA */ CblasNoTrans,
              /* TRANSB */ CblasNoTrans,
              /* M      */ k,
              /* N      */ k,
              /* K      */ m,
              /* ALPHA  */ -1.0,
              /* A      */ C,
              /* LDA    */ n,
              /* B      */ B,
              /* LDB    */ n,
              /* BETA   */ 1.0,
              /* C      */ D,
              /* LDC    */ n);

  // 5: F = F - C * E
  cblas_dgemv(CblasColMajor,
              /* TRANS */ CblasNoTrans,
              /* M     */ k,
              /* N     */ m,
              /* ALPHA */ -1.0,
              /* A     */ C,
              /* LDA   */ n,
              /* X     */ E,
              /* INCX  */ 1,
              /* BETA  */ 1.0,
              /* Y     */ F,
              /* INCY  */ 1);

  // 6.1: Zero matrix A
  for (int i = 0; i < m; i++)
    for (int j = 0; j < m; j++)
      matrix[j * n + i] = 0.0;

  // 6.2: Set 1 on diagonal of A
  for (int i = 0; i < m; i++)
    matrix[i * n + i] = 1.0;

  // 7: Zero matrix C
  for (int i = m; i < n; i++)
    for (int j = 0; j < m; j++)
      matrix[j * n + i] = 0.0;
}

eliminate.h

#include <stdexcept>

class lapack_exception: public std::runtime_error
{
    public:
        lapack_exception(std::string const& msg):
            std::runtime_error(msg)
        {}
};

void eliminate(double * matrix, double * rhs, int n, int m);

partial_elimination.tex

\documentclass[12pt]{article}
\usepackage{amsmath}
\usepackage[T1]{fontenc}
\usepackage[utf8x]{inputenc}

\title{How to do partial elimination with BLAS and LAPACK}
\begin{document}
\maketitle

\section{Problem statement}
Lets assume that we have equation system of size $n$. We want to partially eliminate first $m$ variables from the system, so that $k = n - m$ variables may still be not eliminated. What we need to do is to calculate Shur complement of the equation system.
\begin{equation}
    \begin{bmatrix}
    A_{m × m} & B_{m × k} \\
    C_{k × m} & D_{k × k}
    \end{bmatrix}_{n × n} x = \begin{bmatrix}
    E_{m × 1} \\
    F_{m × 1} 
    \end{bmatrix}_{n × 1}
\longrightarrow
    \begin{bmatrix}
    \boldsymbol{I}_{m × m} & B'_{m × k} \\
    \boldsymbol{0}_{k × m} & D'_{k × k}
    \end{bmatrix}_{n × n} x = \begin{bmatrix}
    E'_{m × 1} \\
    F'_{m × 1} 
    \end{bmatrix}_{n × 1}
\end{equation}

\section{Algorithm}
\begin{align}
\label{step1}    B' &= A^{-1} \times B \\ 
\label{step2}    E' &= A^{-1} \times E \\
\label{step3}    D' &= D - C \times B' \\
\label{step4}    F' &= F - C \times E'
\end{align}


\section{BLAS/LAPACK implementation}

Let's assume that
\begin{itemize}
    \item matrices have continuous columns (like in Foratran), which is not standard behavior in C,
    \item \textbf{G} is LHS matrix of the system that is composed of submatrices \textbf{A} to \textbf{D},
    \item \textbf{H} is RHS vector of the system that is composed of subvectors \textbf{E} and \textbf{F},
    \item \textbf{A} to \textbf{F} are pointers to the begining of the submatrices,
    \item \textbf{n}, \textbf{m} and \textbf{k} are defined as in introduction,
    \item we check if operation succeeded after operations that may fail (e.g. LU decomposition),
    \item the operation will be in-place and will modify input matrix.
\end{itemize}

\paragraph{Algorithm}
\begin{enumerate}
    \item Do the LU decomposition for submatrix A: \\ \texttt{dgetrf(m, m, A, n, ipiv, status)} \\ The row pivot vector returned in \texttt{ipiv} will be needed in next two steps.
    \item $B = (PLU)^{-1} \times B$: \\ \texttt{dgetrs(NoTrans, m, k, A, n, ipiv, B, n, status)}
    \item $E = (PLU)^{-1} \times E$: \\ \texttt{dgetrs(NoTrans, m, 1, A, n, ipiv, E, n, status)}
    \item $D = D - C \times B$: \\ \texttt{dgemm(NoTrans, NoTrans, k, k, m, -1.0, C, n, B, n, 1.0, D, n)}
    \item $F = F - C \times E$: \\ \texttt{dgemv(NoTrans, k, m, -1.0, C, n, E, 1, 1.0, F, 1)}
    \item\label{diag} $A = \boldsymbol{I}_{m × m}$ – $A$ is identity matrix
    \item\label{zero} $C = \boldsymbol{0}_{k × m}$ – $C$ is zero matrix
\end{enumerate}

Steps \ref{diag} and \ref{zero} may be omited if you only need Schur complement (matrices D and F).

\end{document}