lapack/experimental/ilaxly.cpp

#include "../../Eigen/Core"

extern "C" {
int ilaslc_(int *m, int *n, float *a, int *lda);
int iladlc_(int *m, int *n, double *a, int *lda);
int ilaclc_(int *m, int *n, std::complex<float> *a, int *lda);
int ilazlc_(int *m, int *n, std::complex<double> *a, int *lda);
int ilaslr_(int *m, int *n, float *a, int *lda);
int iladlr_(int *m, int *n, double *a, int *lda);
int ilaclr_(int *m, int *n, std::complex<float> *a, int *lda);
int ilazlr_(int *m, int *n, std::complex<double> *a, int *lda);
}

using Eigen::ColMajor;
using Eigen::Dynamic;
using Eigen::Map;
using Eigen::Matrix;
using Eigen::OuterStride;

// Scan a matrix for its last non-zero column.
template <typename T>
int ilaxlc(const int m, const int n, T *a, const int lda) {
  if (n == 0) return 0;
  using MatrixType =
      Map<Matrix<T, Dynamic, Dynamic, ColMajor>, 0, OuterStride<>>;
  MatrixType a_matrix(a, m, n, OuterStride<>(lda));

  // Quick test for the common case where one corner is non-zero.
  if ((a_matrix(0, n - 1) != T(0)) || (a_matrix(m - 1, n - 1) != T(0))) {
    return n;
  }

  // Now scan each column from the end, returning with the first non-zero.
  int c = n - 1;
  for (; c >= 0; --c) {
    if ((a_matrix.array().col(c) != T(0)).any()) {
      break;
    }
  }
  return c + 1;
}

int ilaslc_(int *m, int *n, float *a, int *lda) {
  return ilaxlc(*m, *n, a, *lda);
}

int iladlc_(int *m, int *n, double *a, int *lda) {
  return ilaxlc(*m, *n, a, *lda);
}

int ilaclc_(int *m, int *n, std::complex<float> *a, int *lda) {
  return ilaxlc(*m, *n, a, *lda);
}

int ilazlc_(int *m, int *n, std::complex<double> *a, int *lda) {
  return ilaxlc(*m, *n, a, *lda);
}

// Scan a matrix for its last non-zero row.
template <typename T>
int ilaxlr(const int m, const int n, T *a, const int lda) {
  if (m == 0) return m;
  using MatrixType =
      Map<Matrix<T, Dynamic, Dynamic, ColMajor>, 0, OuterStride<>>;
  MatrixType a_matrix(a, m, n, OuterStride<>(lda));

  // Quick test for the common case where one corner is non-zero.
  if ((a_matrix(m - 1, 0) != T(0)) || (a_matrix(m - 1, n - 1) != T(0))) {
    return m;
  }

  // Now scan each rowfrom the end, returning with the first non-zero.
  int r = m - 1;
  for (; r >= 0; --r) {
    if ((a_matrix.array().row(r) == T(0)).any()) {
      break;
    }
  }
  return r + 1;
}

int ilaslr_(int *m, int *n, float *a, int *lda) {
  return ilaxlr(*m, *n, a, *lda);
}

int iladlr_(int *m, int *n, double *a, int *lda) {
  return ilaxlr(*m, *n, a, *lda);
}

int ilaclr_(int *m, int *n, std::complex<float> *a, int *lda) {
  return ilaxlr(*m, *n, a, *lda);
}

int ilazlr_(int *m, int *n, std::complex<double> *a, int *lda) {
  return ilaxlr(*m, *n, a, *lda);
}