uselapack.h

/* MLPACK 0.2
 *
 * Copyright (c) 2008, 2009 Alexander Gray,
 *                          Garry Boyer,
 *                          Ryan Riegel,
 *                          Nikolaos Vasiloglou,
 *                          Dongryeol Lee,
 *                          Chip Mappus, 
 *                          Nishant Mehta,
 *                          Hua Ouyang,
 *                          Parikshit Ram,
 *                          Long Tran,
 *                          Wee Chin Wong
 *
 * Copyright (c) 2008, 2009 Georgia Institute of Technology
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
 * 02110-1301, USA.
 */
#ifndef LA_USELAPACK_H
#define LA_USELAPACK_H

#include "fastlib/la/matrix.h"
//#include "matrix.h"
#include "fastlib/col/arraylist.h"

#define DEBUG_VECSIZE(a, b) \
    DEBUG_SAME_SIZE((a).length(), (b).length())
#define DEBUG_MATSIZE(a, b) \
    DEBUG_SAME_SIZE((a).n_rows(), (b).n_rows());\
    DEBUG_SAME_SIZE((a).n_cols(), (b).n_cols())
#define DEBUG_MATSQUARE(a) \
    DEBUG_SAME_SIZE((a).n_rows(), (a).n_cols())

/*
 * TODO: this could an overhaul; LAPACK failures may mean either input
 * problems or just some linear algebra result (e.g. singular matrix).
 *
 * Perhaps extend success_t?
 */
#define SUCCESS_FROM_LAPACK(info) \
    (likely((info) == 0) ? SUCCESS_PASS : SUCCESS_FAIL)

#ifndef NO_LAPACK
#define USE_LAPACK
#define USE_BLAS_L1
#endif

#if defined(DOXYGEN) && !defined(USE_LAPACK)
/* Use the doxygen from the blas level 1 code. */
#define USE_LAPACK
#endif

#ifdef USE_LAPACK
#include "fastlib/la/blas.h"
#include "fastlib/la/clapack.h"
//#include "blas.h"
//#include "clapack.h"
#endif

namespace la {
  inline void ScaleRows(index_t n_rows, index_t n_cols,
      const double *scales, double *matrix) {
    do {
      for (index_t i = 0; i < n_rows; i++) {
        matrix[i] *= scales[i];
      }
      matrix += n_rows;
    } while (--n_cols);
  }

#if defined(USE_LAPACK) && defined(USE_BLAS_L1)
  /* --- Wrappers for BLAS level 1 --- */
  /*
   * TODO: Add checking to ensure arrays are non-overlapping.
   */

  inline double LengthEuclidean(index_t length, const double *x) {
    return F77_FUNC(dnrm2)(length, x, 1);
  }

  inline double Dot(index_t length, const double *x, const double *y) {
    return F77_FUNC(ddot)(length, x, 1, y, 1);
  }

  inline void Scale(index_t length, double alpha, double *x) {
    F77_FUNC(dscal)(length, alpha, x, 1);
  }
  inline void ScaleOverwrite(index_t length,
      double alpha, const double *x, double *y) {
    mem::Copy(y, x, length);
    Scale(length, alpha, y);
  }

  inline void AddExpert(index_t length,
      double alpha, const double *x, double *y) {
    F77_FUNC(daxpy)(length, alpha, x, 1, y, 1);
  }

  inline void AddTo(index_t length, const double *x, double *y) {
    AddExpert(length, 1.0, x, y);
  }
  inline void AddOverwrite(index_t length,
      const double *x, const double *y, double *z) {
    mem::Copy(z, y, length);
    AddTo(length, x, z);
  }

  inline void SubFrom(index_t length, const double *x, double *y) {
    AddExpert(length, -1.0, x, y);
  }
  inline void SubOverwrite(index_t length,
      const double *x, const double *y, double *z) {
    mem::Copy(z, y, length);
    SubFrom(length, x, z);
  }
#else
  /* --- Equivalent to BLAS level 1 --- */
  /* These functions steal their comments from the previous versons */
  
  inline double LengthEuclidean(index_t length, const double *x) {
    double sum = 0;
    do {
      sum += *x * *x;
      x++;
    } while (--length);
    return sqrt(sum);
  }

  inline double Dot(index_t length, const double *x, const double *y) {
    double sum = 0;
    const double *end = x + length;
    while (x != end) {
      sum += *x++ * *y++;
    }
    return sum;
  }

  inline void Scale(index_t length, double alpha, double *x) {
    const double *end = x + length;
    while (x != end) {
      *x++ *= alpha;
    }
  }

  inline void ScaleOverwrite(index_t length,
      double alpha, const double *x, double *y) {
    const double *end = x + length;
    while (x != end) {
      *y++ = alpha * *x++;
    }
  }

  inline void AddExpert(index_t length,
      double alpha, const double *x, double *y) {
    const double *end = x + length;
    while (x != end) {
      *y++ += alpha * *x++;
    }
  }

  inline void AddTo(index_t length, const double *x, double *y) {
    const double *end = x + length;
    while (x != end) {
      *y++ += *x++;
    }
  }

  inline void AddOverwrite(index_t length,
      const double *x, const double *y, double *z) {
    for (index_t i = 0; i < length; i++) {
      z[i] = y[i] + x[i];
    }
  }

  inline void SubFrom(index_t length, const double *x, double *y) {
    const double *end = x + length;
    while (x != end) {
      *y++ -= *x++;
    }
  }

  inline void SubOverwrite(index_t length,
      const double *x, const double *y, double *z) {
    for (index_t i = 0; i < length; i++) {
      z[i] = y[i] - x[i];
    }
  }
#endif

  inline double LengthEuclidean(const Vector &x) {
    return LengthEuclidean(x.length(), x.ptr());
  }

  inline double Dot(const Vector &x, const Vector &y) {
    DEBUG_SAME_SIZE(x.length(), y.length());
    return Dot(x.length(), x.ptr(), y.ptr());
  }
  inline double Dot(const Matrix &x, const Matrix &y) {
    DEBUG_SAME_SIZE(x.n_rows(), y.n_rows());
    DEBUG_SAME_SIZE(x.n_cols(), y.n_cols());
    return Dot(x.n_elements(), x.ptr(), y.ptr());
  }

  /* --- Matrix/Vector Scaling --- */

  inline void Scale(double alpha, Vector *x) {
    Scale(x->length(), alpha, x->ptr());
  }
  inline void Scale(double alpha, Matrix *X) {
    Scale(X->n_elements(), alpha, X->ptr());
  }
  inline void ScaleRows(const Vector& d, Matrix *X) {
    DEBUG_SAME_SIZE(d.length(), X->n_rows());
    ScaleRows(d.length(), X->n_cols(), d.ptr(), X->ptr());
  }
  inline void ScaleRows(const Matrix& d, Matrix *X) {
    DEBUG_SAME_SIZE(min(d.n_cols(), d.n_rows()), 1);
    DEBUG_SAME_SIZE(max(d.n_cols(), d.n_rows()), X->n_rows());
    ScaleRows(X->n_rows(), X->n_cols(), d.ptr(), X->ptr());
  }
  inline void ScaleOverwrite(double alpha, const Vector &x, Vector *y) {
    DEBUG_SAME_SIZE(x.length(), y->length());
    ScaleOverwrite(x.length(), alpha, x.ptr(), y->ptr());
  }
  inline void ScaleOverwrite(double alpha, const Matrix &X, Matrix *Y) {
    DEBUG_SAME_SIZE(X.n_rows(), Y->n_rows());
    DEBUG_SAME_SIZE(X.n_cols(), Y->n_cols());
    ScaleOverwrite(X.n_elements(), alpha, X.ptr(), Y->ptr());
  }

  inline void ScaleInit(double alpha, const Vector &x, Vector *y) {
    y->Init(x.length());
    ScaleOverwrite(alpha, x, y);
  }
  inline void ScaleInit(double alpha, const Matrix &X, Matrix *Y) {
    Y->Init(X.n_rows(), X.n_cols());
    ScaleOverwrite(alpha, X, Y);
  }

  /* --- Scaled Matrix/Vector Addition --- */

  inline void AddExpert(double alpha, const Vector &x, Vector *y) {
    DEBUG_SAME_SIZE(x.length(), y->length());
    AddExpert(x.length(), alpha, x.ptr(), y->ptr());
  }
  inline void AddExpert(double alpha, const Matrix &X, Matrix *Y) {
    DEBUG_SAME_SIZE(X.n_rows(), Y->n_rows());
    DEBUG_SAME_SIZE(X.n_cols(), Y->n_cols());
    AddExpert(X.n_elements(), alpha, X.ptr(), Y->ptr());
  }

  /* --- Matrix/Vector Addition --- */

  inline void AddTo(const Vector &x, Vector *y) {
    DEBUG_SAME_SIZE(x.length(), y->length());
    AddTo(x.length(), x.ptr(), y->ptr());
  }
  inline void AddTo(const Matrix &X, Matrix *Y) {
    DEBUG_SAME_SIZE(X.n_rows(), Y->n_rows());
    DEBUG_SAME_SIZE(X.n_cols(), Y->n_cols());
    AddTo(X.n_elements(), X.ptr(), Y->ptr());
  }

  inline void AddOverwrite(const Vector &x, const Vector &y, Vector *z) {
    DEBUG_SAME_SIZE(x.length(), y.length());
    DEBUG_SAME_SIZE(z->length(), y.length());
    AddOverwrite(x.length(), x.ptr(), y.ptr(), z->ptr());
  }
  inline void AddOverwrite(const Matrix &X, const Matrix &Y, Matrix *Z) {
    DEBUG_SAME_SIZE(X.n_rows(), Y.n_rows());
    DEBUG_SAME_SIZE(X.n_cols(), Y.n_cols());
    DEBUG_SAME_SIZE(Z->n_rows(), Y.n_rows());
    DEBUG_SAME_SIZE(Z->n_cols(), Y.n_cols());
    AddOverwrite(X.n_elements(), X.ptr(), Y.ptr(), Z->ptr());
  }

  inline void AddInit(const Vector &x, const Vector &y, Vector *z) {
    z->Init(x.length());
    AddOverwrite(x, y, z);
  }
  inline void AddInit(const Matrix &X, const Matrix &Y, Matrix *Z) {
    Z->Init(X.n_rows(), X.n_cols());
    AddOverwrite(X, Y, Z);
  }

  /* --- Matrix/Vector Subtraction --- */

  inline void SubFrom(const Vector &x, Vector *y) {
    DEBUG_SAME_SIZE(x.length(), y->length());
    SubFrom(x.length(), x.ptr(), y->ptr());
  }
  inline void SubFrom(const Matrix &X, Matrix *Y) {
    DEBUG_SAME_SIZE(X.n_rows(), Y->n_rows());
    DEBUG_SAME_SIZE(X.n_cols(), Y->n_cols());
    SubFrom(X.n_elements(), X.ptr(), Y->ptr());
  }

  inline void SubOverwrite(const Vector &x, const Vector &y, Vector *z) {
    DEBUG_SAME_SIZE(x.length(), y.length());
    DEBUG_SAME_SIZE(z->length(), y.length());
    SubOverwrite(x.length(), x.ptr(), y.ptr(), z->ptr());
  }
  inline void SubOverwrite(const Matrix &X, const Matrix &Y, Matrix *Z) {
    DEBUG_SAME_SIZE(X.n_rows(), Y.n_rows());
    DEBUG_SAME_SIZE(X.n_cols(), Y.n_cols());
    DEBUG_SAME_SIZE(Z->n_rows(), Y.n_rows());
    DEBUG_SAME_SIZE(Z->n_cols(), Y.n_cols());
    SubOverwrite(X.n_elements(), X.ptr(), Y.ptr(), Z->ptr());
  }

  inline void SubInit(const Vector &x, const Vector &y, Vector *z) {
    z->Init(x.length());
    SubOverwrite(x, y, z);
  }
  inline void SubInit(const Matrix &X, const Matrix &Y, Matrix *Z) {
    Z->Init(X.n_rows(), X.n_cols());
    SubOverwrite(X, Y, Z);
  }

  /* --- Matrix Transpose --- */

  /*
   * TODO: These could be an order of magnitude faster if we use the
   * cache-efficient Morton layout matrix transpose algoritihm.
   */

  inline void TransposeSquare(Matrix *X) {
    DEBUG_MATSQUARE(*X);
    index_t nr = X->n_rows();
    for (index_t r = 1; r < nr; r++) {
      for (index_t c = 0; c < r; c++) {
        double temp = X->get(r, c);
        X->set(r, c, X->get(c, r));
        X->set(c, r, temp);        
      }
    }
  }

  inline void TransposeOverwrite(const Matrix &X, Matrix *Y) {
    DEBUG_SAME_SIZE(X.n_rows(), Y->n_cols());
    DEBUG_SAME_SIZE(X.n_cols(), Y->n_rows());
    index_t nr = X.n_rows();
    index_t nc = X.n_cols();
    for (index_t r = 0; r < nr; r++) {
      for (index_t c = 0; c < nc; c++) {
        Y->set(c, r, X.get(r, c));
      }
    }
  }
  inline void TransposeInit(const Matrix &X, Matrix *Y) {
    Y->Init(X.n_cols(), X.n_rows());
    TransposeOverwrite(X, Y);
  }



#ifdef USE_LAPACK
  /* --- Wrappers for BLAS level 2 --- */

  inline void MulExpert(
      double alpha, const Matrix &A, const Vector &x,
      double beta, Vector *y) {
    DEBUG_ASSERT(x.ptr() != y->ptr());
    DEBUG_SAME_SIZE(A.n_cols(), x.length());
    DEBUG_SAME_SIZE(A.n_rows(), y->length());
    F77_FUNC(dgemv)("N", A.n_rows(), A.n_cols(),
        alpha, A.ptr(), A.n_rows(), x.ptr(), 1,
        beta, y->ptr(), 1);
  }

  inline void MulOverwrite(const Matrix &A, const Vector &x, Vector *y) {
    MulExpert(1.0, A, x, 0.0, y);
  }
  inline void MulInit(const Matrix &A, const Vector &x, Vector *y) {
    y->Init(A.n_rows());
    MulOverwrite(A, x, y);
  }

  inline void MulExpert(
      double alpha, const Vector &x, const Matrix &A,
      double beta, Vector *y) {
    DEBUG_ASSERT(x.ptr() != y->ptr());
    DEBUG_SAME_SIZE(A.n_rows(), x.length());
    DEBUG_SAME_SIZE(A.n_cols(), y->length());
    F77_FUNC(dgemv)("T", A.n_rows(), A.n_cols(),
        alpha, A.ptr(), A.n_rows(), x.ptr(), 1,
        beta, y->ptr(), 1);
  }

  inline void MulOverwrite(const Vector &x, const Matrix &A, Vector *y) {
    MulExpert(1.0, x, A, 0.0, y);
  }
  inline void MulInit(const Vector &x, const Matrix &A, Vector *y) {
    y->Init(A.n_cols());
    MulOverwrite(x, A, y);
  }



  /* --- Wrappers for BLAS level 3 --- */

  inline void MulExpert(
      double alpha,
      bool trans_A, const Matrix &A,
      bool trans_B, const Matrix &B,
      double beta, Matrix *C) {
    DEBUG_ASSERT(A.ptr() != C->ptr());
    DEBUG_ASSERT(B.ptr() != C->ptr());
    DEBUG_SAME_SIZE(trans_A ? A.n_rows() : A.n_cols(),
                   trans_B ? B.n_cols() : B.n_rows());
    DEBUG_SAME_SIZE(trans_A ? A.n_cols() : A.n_rows(), C->n_rows());
    DEBUG_SAME_SIZE(trans_B ? B.n_rows() : B.n_cols(), C->n_cols());
    F77_FUNC(dgemm)(trans_A ? "T" : "N", trans_B ? "T" : "N",
        C->n_rows(), C->n_cols(), trans_A ? A.n_rows() : A.n_cols(),
        alpha, A.ptr(), A.n_rows(), B.ptr(), B.n_rows(),
        beta, C->ptr(), C->n_rows());
  }
  inline void MulExpert(
      double alpha, const Matrix &A, const Matrix &B,
      double beta, Matrix *C) {
    MulExpert(alpha, false, A, false, B, beta, C);
  }

  inline void MulOverwrite(const Matrix &A, const Matrix &B, Matrix *C) {
    MulExpert(1.0, A, B, 0.0, C);
  }
  inline void MulInit(const Matrix &A, const Matrix &B, Matrix *C) {
    C->Init(A.n_rows(), B.n_cols());
    MulOverwrite(A, B, C);
  }

  inline void MulTransAOverwrite(
      const Matrix &A, const Matrix &B, Matrix *C) {
    MulExpert(1.0, true, A, false, B, 0.0, C);
  }
  inline void MulTransAInit(
      const Matrix &A, const Matrix &B, Matrix *C) {
    C->Init(A.n_cols(), B.n_cols());
    MulTransAOverwrite(A, B, C);
  }

  inline void MulTransBOverwrite(
      const Matrix &A, const Matrix &B, Matrix *C) {
    MulExpert(1.0, false, A, true, B, 0.0, C);
  }
  inline void MulTransBInit(
      const Matrix &A, const Matrix &B, Matrix *C) {
    C->Init(A.n_rows(), B.n_rows());
    MulTransBOverwrite(A, B, C);
  }



  /* --- Wrappers for LAPACK --- */

  extern int dgetri_block_size;
  extern int dgeqrf_block_size;
  extern int dorgqr_block_size;
  extern int dgeqrf_dorgqr_block_size;

  inline success_t PLUExpert(f77_integer *pivots, Matrix *A_in_LU_out) {
    f77_integer info;
    F77_FUNC(dgetrf)(A_in_LU_out->n_rows(), A_in_LU_out->n_cols(),
        A_in_LU_out->ptr(), A_in_LU_out->n_rows(), pivots, &info);
    return SUCCESS_FROM_LAPACK(info);
  }
  success_t PLUInit(const Matrix &A,
      ArrayList<f77_integer> *pivots, Matrix *L, Matrix *U);

  inline success_t InverseExpert(f77_integer *pivots, Matrix *LU_in_B_out) {
    f77_integer info;
    f77_integer n = LU_in_B_out->n_rows();
    f77_integer lwork = dgetri_block_size * n;
    double work[lwork];
    DEBUG_MATSQUARE(*LU_in_B_out);
    F77_FUNC(dgetri)(n, LU_in_B_out->ptr(), n, pivots, work, lwork, &info);
    return SUCCESS_FROM_LAPACK(info);
  }
  success_t Inverse(Matrix *A);
  success_t InverseOverwrite(const Matrix &A, Matrix *B);
  inline success_t InverseInit(const Matrix &A, Matrix *B) {
    B->Init(A.n_rows(), A.n_cols());
    return InverseOverwrite(A, B);
  }

  long double Determinant(const Matrix &A);
  double DeterminantLog(const Matrix &A, int *sign_out);

  inline success_t SolveExpert(
      f77_integer *pivots, Matrix *A_in_LU_out,
      index_t k, double *B_in_X_out) {
    DEBUG_MATSQUARE(*A_in_LU_out);
    f77_integer info;
    f77_integer n = A_in_LU_out->n_rows();
    F77_FUNC(dgesv)(n, k, A_in_LU_out->ptr(), n, pivots,
        B_in_X_out, n, &info);
    return SUCCESS_FROM_LAPACK(info);
  }
  inline success_t SolveInit(const Matrix &A, const Matrix &B, Matrix *X) {
    DEBUG_MATSQUARE(A);
    DEBUG_SAME_SIZE(A.n_rows(), B.n_rows());
    Matrix tmp;
    index_t n = B.n_rows();
    f77_integer pivots[n];
    tmp.Copy(A);
    X->Copy(B);
    return SolveExpert(pivots, &tmp, B.n_cols(), X->ptr());
  }
  inline success_t SolveInit(const Matrix &A, const Vector &b, Vector *x) {
    DEBUG_MATSQUARE(A);
    DEBUG_SAME_SIZE(A.n_rows(), b.length());
    Matrix tmp;
    index_t n = b.length();
    f77_integer pivots[n];
    tmp.Copy(A);
    x->Copy(b);
    return SolveExpert(pivots, &tmp, 1, x->ptr());
  }

  success_t QRExpert(Matrix *A_in_Q_out, Matrix *R);
  success_t QRInit(const Matrix &A, Matrix *Q, Matrix *R);

  success_t SchurExpert(Matrix *A_in_T_out,
      double *w_real, double *w_imag, double *Z);
  inline success_t SchurInit(const Matrix &A,
      Vector *w_real, Vector *w_imag, Matrix *T, Matrix *Z) {
    index_t n = A.n_rows();
    T->Copy(A);
    w_real->Init(n);
    w_imag->Init(n);
    Z->Init(n, n);
    return SchurExpert(T, w_real->ptr(), w_imag->ptr(), Z->ptr());
  }

  success_t EigenExpert(Matrix *A_garbage,
      double *w_real, double *w_imag, double *V_raw);

  inline success_t EigenvaluesInit(const Matrix &A,
      Vector *w_real, Vector *w_imag) {
    DEBUG_MATSQUARE(A);
    int n = A.n_rows();
    w_real->Init(n);
    w_imag->Init(n);
    Matrix tmp;
    tmp.Copy(A);
    return SchurExpert(&tmp, w_real->ptr(), w_imag->ptr(), NULL);
  }
  success_t EigenvaluesInit(const Matrix &A, Vector *w);

  success_t EigenvectorsInit(const Matrix &A,
      Vector *w_real, Vector *w_imag, Matrix *V_real, Matrix *V_imag);

  success_t EigenvectorsInit(const Matrix &A, Vector *w, Matrix *V);

  success_t GenEigenSymmetric(int itype, Matrix *A_eigenvec, Matrix *B_chol, double *w);

  success_t GenEigenNonSymmetric(Matrix *A_garbage, Matrix *B_garbage,
      double *alpha_real, double *alpha_imag, double *beta, double *V_raw);

  success_t SVDExpert(Matrix* A_garbage, double *s, double *U, double *VT);

  inline success_t SVDInit(const Matrix &A, Vector *s) {
    s->Init(std::min(A.n_rows(), A.n_cols()));
    Matrix tmp;
    tmp.Copy(A);
    return SVDExpert(&tmp, s->ptr(), NULL, NULL);
  }

  inline success_t SVDInit(const Matrix &A, Vector *s, Matrix *U, Matrix *VT) {
    index_t k = std::min(A.n_rows(), A.n_cols());
    s->Init(k);
    U->Init(A.n_rows(), k);
    VT->Init(k, A.n_cols());
    Matrix tmp;
    tmp.Copy(A);
    return SVDExpert(&tmp, s->ptr(), U->ptr(), VT->ptr());
  }

  success_t Cholesky(Matrix *A_in_U_out);

  inline success_t CholeskyInit(const Matrix &A, Matrix *U) {
    U->Copy(A);
    return Cholesky(U);
  }
#endif

  /*
   * TODO:
   *   symmetric eigenvalue
   *   http://www.netlib.org/lapack/lug/node71.html
   *   dgels for least-squares problems
   *   dgesv for linear equations
   */
};

#endif