// matrix/tp-matrix.cc
  
  // Copyright 2009-2011  Ondrej Glembek;  Lukas Burget;  Microsoft Corporation
  //                      Saarland University;  Yanmin Qian;   Haihua Xu
  
  // See ../../COPYING for clarification regarding multiple authors
  //
  // Licensed under the Apache License, Version 2.0 (the "License");
  // you may not use this file except in compliance with the License.
  // You may obtain a copy of the License at
  
  //  http://www.apache.org/licenses/LICENSE-2.0
  
  // THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
  // KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
  // WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
  // MERCHANTABLITY OR NON-INFRINGEMENT.
  // See the Apache 2 License for the specific language governing permissions and
  // limitations under the License.
  
  #include "matrix/tp-matrix.h"
  #include "matrix/sp-matrix.h"
  #include "matrix/kaldi-matrix.h"
  #include "matrix/cblas-wrappers.h"
  
  
  namespace kaldi {
  
  #ifndef HAVE_ATLAS
  template<typename Real>
  void TpMatrix<Real>::Invert() {
    // these are CLAPACK types
    KaldiBlasInt result;
    KaldiBlasInt rows = static_cast<int>(this->num_rows_);
  
    // clapack call
    // NOTE: Even though "U" is for upper, lapack assumes column-wise storage
    // of the data. We have a row-wise storage, therefore, we need to "invert"
    clapack_Xtptri(&rows, this->data_, &result);
  
    if (result < 0) {
      KALDI_ERR << "Call to CLAPACK stptri_ function failed";
    } else if (result > 0) {
      KALDI_ERR << "Matrix is singular";
    }
  }
  #else
  template<typename Real>
  void TpMatrix<Real>::Invert() {
    // ATLAS doesn't implement triangular matrix inversion in packed
    // format, so we temporarily put in non-packed format.
    Matrix<Real> tmp(*this);
    int rows = static_cast<int>(this->num_rows_);
  
    // ATLAS call.  It's really row-major ordering and a lower triangular matrix,
    // but there is some weirdness with Fortran-style indexing that we need to
    // take account of, so everything gets swapped.
    int result = clapack_Xtrtri( rows, tmp.Data(), tmp.Stride());
    // Let's hope ATLAS has the same return value conventions as clapack.
    // I couldn't find any documentation online.
    if (result < 0) {
      KALDI_ERR << "Call to ATLAS strtri function failed";
    } else if (result > 0) {
      KALDI_ERR << "Matrix is singular";
    }
    (*this).CopyFromMat(tmp);
  }
  #endif
  
  template<typename Real>
  Real TpMatrix<Real>::Determinant() {
    double   det = 1.0;
    for (MatrixIndexT i = 0; i<this->NumRows(); i++) {
      det *= (*this)(i, i);
    }
    return static_cast<Real>(det);
  }
  
  
  template<typename Real>
  void TpMatrix<Real>::Swap(TpMatrix<Real> *other) {
    std::swap(this->data_, other->data_);
    std::swap(this->num_rows_, other->num_rows_);
  }
  
  
  template<typename Real>
  void TpMatrix<Real>::Cholesky(const SpMatrix<Real> &orig) {
    KALDI_ASSERT(orig.NumRows() == this->NumRows());
    MatrixIndexT n = this->NumRows();
    this->SetZero();
    Real *data = this->data_, *jdata = data;  // start of j'th row of matrix.
    const Real *orig_jdata = orig.Data(); // start of j'th row of matrix.
    for (MatrixIndexT j = 0; j < n; j++, jdata += j, orig_jdata += j) {
      Real *kdata = data; // start of k'th row of matrix.
      Real d(0.0);
      for (MatrixIndexT k = 0; k < j; k++, kdata += k) {
        Real s = cblas_Xdot(k, kdata, 1, jdata, 1);
        // (*this)(j, k) = s = (orig(j, k) - s)/(*this)(k, k);
        jdata[k] = s = (orig_jdata[k] - s)/kdata[k];
        d = d + s*s;
      }
      // d = orig(j, j) - d;
      d = orig_jdata[j] - d;
  
      if (d >= 0.0) {
        // (*this)(j, j) = std::sqrt(d);
        jdata[j] = std::sqrt(d);
      } else {
        KALDI_ERR << "Cholesky decomposition failed. Maybe matrix "
            "is not positive definite.";
      }
    }
  }
  
  template<typename Real>
  void TpMatrix<Real>::CopyFromMat(const MatrixBase<Real> &M,
                                   MatrixTransposeType Trans) {
    if (Trans == kNoTrans) {
      KALDI_ASSERT(this->NumRows() == M.NumRows() && M.NumRows() == M.NumCols());
      MatrixIndexT D = this->NumRows();
      const Real *in_i = M.Data();
      MatrixIndexT stride = M.Stride();
      Real *out_i = this->data_;
      for (MatrixIndexT i = 0; i < D; i++, in_i += stride, out_i += i)
        for (MatrixIndexT j = 0; j <= i; j++)
          out_i[j] = in_i[j];
    } else {
      KALDI_ASSERT(this->NumRows() == M.NumRows() && M.NumRows() == M.NumCols());
      MatrixIndexT D = this->NumRows();
      const Real *in_i = M.Data();
      MatrixIndexT stride = M.Stride();
      Real *out_i = this->data_;
      for (MatrixIndexT i = 0; i < D; i++, in_i++, out_i += i) {
        for (MatrixIndexT j = 0; j <= i; j++)
          out_i[j] = in_i[stride*j];
      }
    }
  }
  
  
  template class TpMatrix<float>;
  template class TpMatrix<double>;
  
  }  // namespace kaldi