// matrix/kaldi-matrix.h
  
  // Copyright 2009-2011  Ondrej Glembek;  Microsoft Corporation;  Lukas Burget;
  //                      Saarland University;  Petr Schwarz;  Yanmin Qian;
  //                      Karel Vesely;  Go Vivace Inc.;  Haihua Xu
  //           2017       Shiyin Kang
  //           2019       Yiwen Shao
  
  // 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.
  
  #ifndef KALDI_MATRIX_KALDI_MATRIX_H_
  #define KALDI_MATRIX_KALDI_MATRIX_H_ 1
  
  #include <algorithm>
  
  #include "matrix/matrix-common.h"
  
  namespace kaldi {
  
  /// @{ \addtogroup matrix_funcs_scalar
  
  /// We need to declare this here as it will be a friend function.
  /// tr(A B), or tr(A B^T).
  template<typename Real>
  Real TraceMatMat(const MatrixBase<Real> &A, const MatrixBase<Real> &B,
                   MatrixTransposeType trans = kNoTrans);
  /// @}
  
  /// \addtogroup matrix_group
  /// @{
  
  /// Base class which provides matrix operations not involving resizing
  /// or allocation.   Classes Matrix and SubMatrix inherit from it and take care
  /// of allocation and resizing.
  template<typename Real>
  class MatrixBase {
   public:
    // so this child can access protected members of other instances.
    friend class Matrix<Real>;
    // friend declarations for CUDA matrices (see ../cudamatrix/)
    friend class CuMatrixBase<Real>;
    friend class CuMatrix<Real>;
    friend class CuSubMatrix<Real>;
    friend class CuPackedMatrix<Real>;
    friend class PackedMatrix<Real>;
    friend class SparseMatrix<Real>;
    friend class SparseMatrix<float>;
    friend class SparseMatrix<double>;
  
    /// Returns number of rows (or zero for empty matrix).
    inline MatrixIndexT  NumRows() const { return num_rows_; }
  
    /// Returns number of columns (or zero for empty matrix).
    inline MatrixIndexT NumCols() const { return num_cols_; }
  
    /// Stride (distance in memory between each row).  Will be >= NumCols.
    inline MatrixIndexT Stride() const {  return stride_; }
  
    /// Returns size in bytes of the data held by the matrix.
    size_t  SizeInBytes() const {
      return static_cast<size_t>(num_rows_) * static_cast<size_t>(stride_) *
          sizeof(Real);
    }
  
    /// Gives pointer to raw data (const).
    inline const Real* Data() const {
      return data_;
    }
  
    /// Gives pointer to raw data (non-const).
    inline Real* Data() { return data_; }
  
    /// Returns pointer to data for one row (non-const)
    inline  Real* RowData(MatrixIndexT i) {
      KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(i) <
                   static_cast<UnsignedMatrixIndexT>(num_rows_));
      return data_ + i * stride_;
    }
  
    /// Returns pointer to data for one row (const)
    inline const Real* RowData(MatrixIndexT i) const {
      KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(i) <
                   static_cast<UnsignedMatrixIndexT>(num_rows_));
      return data_ + i * stride_;
    }
  
    /// Indexing operator, non-const
    /// (only checks sizes if compiled with -DKALDI_PARANOID)
    inline Real&  operator() (MatrixIndexT r, MatrixIndexT c) {
      KALDI_PARANOID_ASSERT(static_cast<UnsignedMatrixIndexT>(r) <
                            static_cast<UnsignedMatrixIndexT>(num_rows_) &&
                            static_cast<UnsignedMatrixIndexT>(c) <
                            static_cast<UnsignedMatrixIndexT>(num_cols_));
      return *(data_ + r * stride_ + c);
    }
    /// Indexing operator, provided for ease of debugging (gdb doesn't work
    /// with parenthesis operator).
    Real &Index (MatrixIndexT r, MatrixIndexT c) {  return (*this)(r, c); }
  
    /// Indexing operator, const
    /// (only checks sizes if compiled with -DKALDI_PARANOID)
    inline const Real operator() (MatrixIndexT r, MatrixIndexT c) const {
      KALDI_PARANOID_ASSERT(static_cast<UnsignedMatrixIndexT>(r) <
                            static_cast<UnsignedMatrixIndexT>(num_rows_) &&
                            static_cast<UnsignedMatrixIndexT>(c) <
                            static_cast<UnsignedMatrixIndexT>(num_cols_));
      return *(data_ + r * stride_ + c);
    }
  
    /*   Basic setting-to-special values functions. */
  
    /// Sets matrix to zero.
    void SetZero();
    /// Sets all elements to a specific value.
    void Set(Real);
    /// Sets to zero, except ones along diagonal [for non-square matrices too]
    void SetUnit();
    /// Sets to random values of a normal distribution
    void SetRandn();
    /// Sets to numbers uniformly distributed on (0, 1)
    void SetRandUniform();
  
    /*  Copying functions.  These do not resize the matrix! */
  
  
    /// Copy given matrix. (no resize is done).
    template<typename OtherReal>
    void CopyFromMat(const MatrixBase<OtherReal> & M,
                     MatrixTransposeType trans = kNoTrans);
  
    /// Copy from compressed matrix.
    void CopyFromMat(const CompressedMatrix &M);
  
    /// Copy given spmatrix. (no resize is done).
    template<typename OtherReal>
    void CopyFromSp(const SpMatrix<OtherReal> &M);
  
    /// Copy given tpmatrix. (no resize is done).
    template<typename OtherReal>
    void CopyFromTp(const TpMatrix<OtherReal> &M,
                    MatrixTransposeType trans = kNoTrans);
  
    /// Copy from CUDA matrix.  Implemented in ../cudamatrix/cu-matrix.h
    template<typename OtherReal>
    void CopyFromMat(const CuMatrixBase<OtherReal> &M,
                     MatrixTransposeType trans = kNoTrans);
  
    /// This function has two modes of operation.  If v.Dim() == NumRows() *
    /// NumCols(), then treats the vector as a row-by-row concatenation of a
    /// matrix and copies to *this.
    /// if v.Dim() == NumCols(), it sets each row of *this to a copy of v.
    void CopyRowsFromVec(const VectorBase<Real> &v);
  
    /// This version of CopyRowsFromVec is implemented in ../cudamatrix/cu-vector.cc
    void CopyRowsFromVec(const CuVectorBase<Real> &v);
  
    template<typename OtherReal>
    void CopyRowsFromVec(const VectorBase<OtherReal> &v);
  
    /// Copies vector into matrix, column-by-column.
    /// Note that rv.Dim() must either equal NumRows()*NumCols() or NumRows();
    /// this has two modes of operation.
    void CopyColsFromVec(const VectorBase<Real> &v);
  
    /// Copy vector into specific column of matrix.
    void CopyColFromVec(const VectorBase<Real> &v, const MatrixIndexT col);
    /// Copy vector into specific row of matrix.
    void CopyRowFromVec(const VectorBase<Real> &v, const MatrixIndexT row);
    /// Copy vector into diagonal of matrix.
    void CopyDiagFromVec(const VectorBase<Real> &v);
  
    /* Accessing of sub-parts of the matrix. */
  
    /// Return specific row of matrix [const].
    inline const SubVector<Real> Row(MatrixIndexT i) const {
      KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(i) <
                   static_cast<UnsignedMatrixIndexT>(num_rows_));
      return SubVector<Real>(data_ + (i * stride_), NumCols());
    }
  
    /// Return specific row of matrix.
    inline SubVector<Real> Row(MatrixIndexT i) {
      KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(i) <
                   static_cast<UnsignedMatrixIndexT>(num_rows_));
      return SubVector<Real>(data_ + (i * stride_), NumCols());
    }
  
    /// Return a sub-part of matrix.
    inline SubMatrix<Real> Range(const MatrixIndexT row_offset,
                                 const MatrixIndexT num_rows,
                                 const MatrixIndexT col_offset,
                                 const MatrixIndexT num_cols) const {
      return SubMatrix<Real>(*this, row_offset, num_rows,
                             col_offset, num_cols);
    }
    inline SubMatrix<Real> RowRange(const MatrixIndexT row_offset,
                                    const MatrixIndexT num_rows) const {
      return SubMatrix<Real>(*this, row_offset, num_rows, 0, num_cols_);
    }
    inline SubMatrix<Real> ColRange(const MatrixIndexT col_offset,
                                    const MatrixIndexT num_cols) const {
      return SubMatrix<Real>(*this, 0, num_rows_, col_offset, num_cols);
    }
  
    /* Various special functions. */
    /// Returns sum of all elements in matrix.
    Real Sum() const;
    /// Returns trace of matrix.
    Real Trace(bool check_square = true) const;
    // If check_square = true, will crash if matrix is not square.
  
    /// Returns maximum element of matrix.
    Real Max() const;
    /// Returns minimum element of matrix.
    Real Min() const;
  
    /// Element by element multiplication with a given matrix.
    void MulElements(const MatrixBase<Real> &A);
  
    /// Divide each element by the corresponding element of a given matrix.
    void DivElements(const MatrixBase<Real> &A);
  
    /// Multiply each element with a scalar value.
    void Scale(Real alpha);
  
    /// Set, element-by-element, *this = max(*this, A)
    void Max(const MatrixBase<Real> &A);
    /// Set, element-by-element, *this = min(*this, A)
    void Min(const MatrixBase<Real> &A);
  
    /// Equivalent to (*this) = (*this) * diag(scale).  Scaling
    /// each column by a scalar taken from that dimension of the vector.
    void MulColsVec(const VectorBase<Real> &scale);
  
    /// Equivalent to (*this) = diag(scale) * (*this).  Scaling
    /// each row by a scalar taken from that dimension of the vector.
    void MulRowsVec(const VectorBase<Real> &scale);
  
    /// Divide each row into src.NumCols() equal groups, and then scale i'th row's
    /// j'th group of elements by src(i, j).  Requires src.NumRows() ==
    /// this->NumRows() and this->NumCols() % src.NumCols() == 0.
    void MulRowsGroupMat(const MatrixBase<Real> &src);
  
    /// Returns logdet of matrix.
    Real LogDet(Real *det_sign = NULL) const;
  
    /// matrix inverse.
    /// if inverse_needed = false, will fill matrix with garbage.
    /// (only useful if logdet wanted).
    void Invert(Real *log_det = NULL, Real *det_sign = NULL,
                bool inverse_needed = true);
    /// matrix inverse [double].
    /// if inverse_needed = false, will fill matrix with garbage
    /// (only useful if logdet wanted).
    /// Does inversion in double precision even if matrix was not double.
    void InvertDouble(Real *LogDet = NULL, Real *det_sign = NULL,
                        bool inverse_needed = true);
  
    /// Inverts all the elements of the matrix
    void InvertElements();
  
    /// Transpose the matrix.  This one is only
    /// applicable to square matrices (the one in the
    /// Matrix child class works also for non-square.
    void Transpose();
  
    /// Copies column r from column indices[r] of src.
    /// As a special case, if indexes[i] == -1, sets column i to zero.
    /// all elements of "indices" must be in [-1, src.NumCols()-1],
    /// and src.NumRows() must equal this.NumRows()
    void CopyCols(const MatrixBase<Real> &src,
                  const MatrixIndexT *indices);
  
    /// Copies row r from row indices[r] of src (does nothing
    /// As a special case, if indexes[i] == -1, sets row i to zero.
    /// all elements of "indices" must be in [-1, src.NumRows()-1],
    /// and src.NumCols() must equal this.NumCols()
    void CopyRows(const MatrixBase<Real> &src,
                  const MatrixIndexT *indices);
  
    /// Add column indices[r] of src to column r.
    /// As a special case, if indexes[i] == -1, skip column i
    /// indices.size() must equal this->NumCols(),
    /// all elements of "reorder" must be in [-1, src.NumCols()-1],
    /// and src.NumRows() must equal this.NumRows()
    void AddCols(const MatrixBase<Real> &src,
                 const MatrixIndexT *indices);
  
    /// Copies row r of this matrix from an array of floats at the location given
    /// by src[r]. If any src[r] is NULL then this.Row(r) will be set to zero.
    /// Note: we are using "pointer to const pointer to const object" for "src",
    ///       because we may create "src" by calling Data() of const CuArray
    void CopyRows(const Real *const *src);
  
    /// Copies row r of this matrix to the array of floats at the location given
    /// by dst[r]. If dst[r] is NULL, does not copy anywhere.  Requires that none
    /// of the memory regions pointed to by the pointers in "dst" overlap (e.g.
    /// none of the pointers should be the same).
    void CopyToRows(Real *const *dst) const;
  
    /// Does for each row r, this.Row(r) += alpha * src.row(indexes[r]).
    /// If indexes[r] < 0, does not add anything. all elements of "indexes" must
    /// be in [-1, src.NumRows()-1], and src.NumCols() must equal this.NumCols().
    void AddRows(Real alpha,
                 const MatrixBase<Real> &src,
                 const MatrixIndexT *indexes);
  
    /// Does for each row r, this.Row(r) += alpha * src[r], treating src[r] as the
    /// beginning of a region of memory representing a vector of floats, of the
    /// same length as this.NumCols(). If src[r] is NULL, does not add anything.
    void AddRows(Real alpha, const Real *const *src);
  
    /// For each row r of this matrix, adds it (times alpha) to the array of
    /// floats at the location given by dst[r]. If dst[r] is NULL, does not do
    /// anything for that row. Requires that none of the memory regions pointed
    /// to by the pointers in "dst" overlap (e.g. none of the pointers should be
    /// the same).
    void AddToRows(Real alpha, Real *const *dst) const;
  
    /// For each row i of *this, adds this->Row(i) to
    /// dst->Row(indexes(i)) if indexes(i) >= 0, else do nothing.
    /// Requires that all the indexes[i] that are >= 0
    /// be distinct, otherwise the behavior is undefined.
    void AddToRows(Real alpha,
                   const MatrixIndexT *indexes,
                   MatrixBase<Real> *dst) const;
  
    inline void ApplyPow(Real power) {
      this -> Pow(*this, power);
    }
  
  
    inline void ApplyPowAbs(Real power, bool include_sign=false) {
      this -> PowAbs(*this, power, include_sign);
    }
  
    inline void ApplyHeaviside() {
      this -> Heaviside(*this);
    }
  
    inline void ApplyFloor(Real floor_val) {
      this -> Floor(*this, floor_val);
    }
  
    inline void ApplyCeiling(Real ceiling_val) {
      this -> Ceiling(*this, ceiling_val);
    }
  
    inline void ApplyExp() {
      this -> Exp(*this);
    }
  
    inline void ApplyExpSpecial() {
      this -> ExpSpecial(*this);
    }
  
    inline void ApplyExpLimited(Real lower_limit, Real upper_limit) {
      this -> ExpLimited(*this, lower_limit, upper_limit);
    }
  
    inline void ApplyLog() {
      this -> Log(*this);
    }
  
    /// Eigenvalue Decomposition of a square NxN matrix into the form (*this) = P D
    /// P^{-1}.  Be careful: the relationship of D to the eigenvalues we output is
    /// slightly complicated, due to the need for P to be real.  In the symmetric
    /// case D is diagonal and real, but in
    /// the non-symmetric case there may be complex-conjugate pairs of eigenvalues.
    /// In this case, for the equation (*this) = P D P^{-1} to hold, D must actually
    /// be block diagonal, with 2x2 blocks corresponding to any such pairs.  If a
    /// pair is lambda +- i*mu, D will have a corresponding 2x2 block
    /// [lambda, mu; -mu, lambda].
    /// Note that if the input matrix (*this) is non-invertible, P may not be invertible
    /// so in this case instead of the equation (*this) = P D P^{-1} holding, we have
    /// instead (*this) P = P D.
    ///
    /// The non-member function CreateEigenvalueMatrix creates D from eigs_real and eigs_imag.
    void Eig(MatrixBase<Real> *P,
             VectorBase<Real> *eigs_real,
             VectorBase<Real> *eigs_imag) const;
  
    /// The Power method attempts to take the matrix to a power using a method that
    /// works in general for fractional and negative powers.  The input matrix must
    /// be invertible and have reasonable condition (or we don't guarantee the
    /// results.  The method is based on the eigenvalue decomposition.  It will
    /// return false and leave the matrix unchanged, if at entry the matrix had
    /// real negative eigenvalues (or if it had zero eigenvalues and the power was
    /// negative).
    bool Power(Real pow);
  
    /** Singular value decomposition
       Major limitations:
       For nonsquare matrices, we assume m>=n (NumRows >= NumCols), and we return
       the "skinny" Svd, i.e. the matrix in the middle is diagonal, and the
       one on the left is rectangular.
  
       In Svd, *this = U*diag(S)*Vt.
       Null pointers for U and/or Vt at input mean we do not want that output.  We
       expect that S.Dim() == m, U is either NULL or m by n,
       and v is either NULL or n by n.
       The singular values are not sorted (use SortSvd for that).  */
    void DestructiveSvd(VectorBase<Real> *s, MatrixBase<Real> *U,
                        MatrixBase<Real> *Vt);  // Destroys calling matrix.
  
    /// Compute SVD (*this) = U diag(s) Vt.   Note that the V in the call is already
    /// transposed; the normal formulation is U diag(s) V^T.
    /// Null pointers for U or V mean we don't want that output (this saves
    /// compute).  The singular values are not sorted (use SortSvd for that).
    void Svd(VectorBase<Real> *s, MatrixBase<Real> *U,
             MatrixBase<Real> *Vt) const;
    /// Compute SVD but only retain the singular values.
    void Svd(VectorBase<Real> *s) const { Svd(s, NULL, NULL); }
  
  
    /// Returns smallest singular value.
    Real MinSingularValue() const {
      Vector<Real> tmp(std::min(NumRows(), NumCols()));
      Svd(&tmp);
      return tmp.Min();
    }
  
    void TestUninitialized() const; // This function is designed so that if any element
    // if the matrix is uninitialized memory, valgrind will complain.
  
    /// Returns condition number by computing Svd.  Works even if cols > rows.
    /// Returns infinity if all singular values are zero.
    Real Cond() const;
  
    /// Returns true if matrix is Symmetric.
    bool IsSymmetric(Real cutoff = 1.0e-05) const;  // replace magic number
  
    /// Returns true if matrix is Diagonal.
    bool IsDiagonal(Real cutoff = 1.0e-05) const;  // replace magic number
  
    /// Returns true if the matrix is all zeros, except for ones on diagonal.  (it
    /// does not have to be square).  More specifically, this function returns
    /// false if for any i, j, (*this)(i, j) differs by more than cutoff from the
    /// expression (i == j ? 1 : 0).
    bool IsUnit(Real cutoff = 1.0e-05) const;     // replace magic number
  
    /// Returns true if matrix is all zeros.
    bool IsZero(Real cutoff = 1.0e-05) const;     // replace magic number
  
    /// Frobenius norm, which is the sqrt of sum of square elements.  Same as Schatten 2-norm,
    /// or just "2-norm".
    Real FrobeniusNorm() const;
  
    /// Returns true if ((*this)-other).FrobeniusNorm()
    /// <= tol * (*this).FrobeniusNorm().
    bool ApproxEqual(const MatrixBase<Real> &other, float tol = 0.01) const;
  
    /// Tests for exact equality.  It's usually preferable to use ApproxEqual.
    bool Equal(const MatrixBase<Real> &other) const;
  
    /// largest absolute value.
    Real LargestAbsElem() const;  // largest absolute value.
  
    /// Returns log(sum(exp())) without exp overflow
    /// If prune > 0.0, it uses a pruning beam, discarding
    /// terms less than (max - prune).  Note: in future
    /// we may change this so that if prune = 0.0, it takes
    /// the max, so use -1 if you don't want to prune.
    Real LogSumExp(Real prune = -1.0) const;
  
    /// Apply soft-max to the collection of all elements of the
    /// matrix and return normalizer (log sum of exponentials).
    Real ApplySoftMax();
  
    /// Set each element to the sigmoid of the corresponding element of "src".
    void Sigmoid(const MatrixBase<Real> &src);
  
    /// Sets each element to the Heaviside step function (x > 0 ? 1 : 0) of the
    /// corresponding element in "src".  Note: in general you can make different
    /// choices for x = 0, but for now please leave it as it (i.e. returning zero)
    /// because it affects the RectifiedLinearComponent in the neural net code.
    void Heaviside(const MatrixBase<Real> &src);
  
    void Exp(const MatrixBase<Real> &src);
  
    void Pow(const MatrixBase<Real> &src, Real power);
  
    void Log(const MatrixBase<Real> &src);
  
    /// Apply power to the absolute value of each element.
    /// If include_sign is true, the result will be multiplied with
    /// the sign of the input value.
    /// If the power is negative and the input to the power is zero,
    /// The output will be set zero. If include_sign is true, it will
    /// multiply the result by the sign of the input.
    void PowAbs(const MatrixBase<Real> &src, Real power, bool include_sign=false);
  
    void Floor(const MatrixBase<Real> &src, Real floor_val);
  
    void Ceiling(const MatrixBase<Real> &src, Real ceiling_val);
  
    /// For each element x of the matrix, set it to
    /// (x < 0 ? exp(x) : x + 1).  This function is used
    /// in our RNNLM training.
    void ExpSpecial(const MatrixBase<Real> &src);
  
    /// This is equivalent to running:
    /// Floor(src, lower_limit);
    /// Ceiling(src, upper_limit);
    /// Exp(src)
    void ExpLimited(const MatrixBase<Real> &src, Real lower_limit, Real upper_limit);
  
    /// Set each element to y = log(1 + exp(x))
    void SoftHinge(const MatrixBase<Real> &src);
  
    /// Apply the function y(i) = (sum_{j = i*G}^{(i+1)*G-1} x_j^(power))^(1 / p).
    /// Requires src.NumRows() == this->NumRows() and  src.NumCols() % this->NumCols() == 0.
    void GroupPnorm(const MatrixBase<Real> &src, Real power);
  
    /// Calculate derivatives for the GroupPnorm function above...
    /// if "input" is the input to the GroupPnorm function above (i.e. the "src" variable),
    /// and "output" is the result of the computation (i.e. the "this" of that function
    /// call), and *this has the same dimension as "input", then it sets each element
    /// of *this to the derivative d(output-elem)/d(input-elem) for each element of "input", where
    /// "output-elem" is whichever element of output depends on that input element.
    void GroupPnormDeriv(const MatrixBase<Real> &input, const MatrixBase<Real> &output,
                         Real power);
  
    /// Apply the function y(i) = (max_{j = i*G}^{(i+1)*G-1} x_j
    /// Requires src.NumRows() == this->NumRows() and  src.NumCols() % this->NumCols() == 0.
    void GroupMax(const MatrixBase<Real> &src);
  
    /// Calculate derivatives for the GroupMax function above, where
    /// "input" is the input to the GroupMax function above (i.e. the "src" variable),
    /// and "output" is the result of the computation (i.e. the "this" of that function
    /// call), and *this must have the same dimension as "input". Each element
    /// of *this will be set to 1 if the corresponding input equals the output of
    /// the group, and 0 otherwise. The equals the function derivative where it is
    /// defined (it's not defined where multiple inputs in the group are equal to the output).
    void GroupMaxDeriv(const MatrixBase<Real> &input, const MatrixBase<Real> &output);
  
    /// Set each element to the tanh of the corresponding element of "src".
    void Tanh(const MatrixBase<Real> &src);
  
    // Function used in backpropagating derivatives of the sigmoid function:
    // element-by-element, set *this = diff * value * (1.0 - value).
    void DiffSigmoid(const MatrixBase<Real> &value,
                     const MatrixBase<Real> &diff);
  
    // Function used in backpropagating derivatives of the tanh function:
    // element-by-element, set *this = diff * (1.0 - value^2).
    void DiffTanh(const MatrixBase<Real> &value,
                  const MatrixBase<Real> &diff);
  
    /** Uses Svd to compute the eigenvalue decomposition of a symmetric positive
     * semi-definite matrix: (*this) = rP * diag(rS) * rP^T, with rP an
     * orthogonal matrix so rP^{-1} = rP^T.   Throws exception if input was not
     * positive semi-definite (check_thresh controls how stringent the check is;
     * set it to 2 to ensure it won't ever complain, but it will zero out negative
     * dimensions in your matrix.
     *
     * Caution: if you want the eigenvalues, it may make more sense to convert to
     * SpMatrix and use Eig() function there, which uses eigenvalue decomposition
     * directly rather than SVD.
    */
    void SymPosSemiDefEig(VectorBase<Real> *s, MatrixBase<Real> *P,
                          Real check_thresh = 0.001);
  
    friend Real kaldi::TraceMatMat<Real>(const MatrixBase<Real> &A,
        const MatrixBase<Real> &B, MatrixTransposeType trans);  // tr (A B)
  
    // so it can get around const restrictions on the pointer to data_.
    friend class SubMatrix<Real>;
  
    /// Add a scalar to each element
    void Add(const Real alpha);
  
    /// Add a scalar to each diagonal element.
    void AddToDiag(const Real alpha);
  
    /// *this += alpha * a * b^T
    template<typename OtherReal>
    void AddVecVec(const Real alpha, const VectorBase<OtherReal> &a,
                   const VectorBase<OtherReal> &b);
  
    /// [each row of *this] += alpha * v
    template<typename OtherReal>
    void AddVecToRows(const Real alpha, const VectorBase<OtherReal> &v);
  
    /// [each col of *this] += alpha * v
    template<typename OtherReal>
    void AddVecToCols(const Real alpha, const VectorBase<OtherReal> &v);
  
    /// *this += alpha * M [or M^T]
    void AddMat(const Real alpha, const MatrixBase<Real> &M,
                MatrixTransposeType transA = kNoTrans);
  
    /// *this += alpha * A [or A^T].
    void AddSmat(Real alpha, const SparseMatrix<Real> &A,
                 MatrixTransposeType trans = kNoTrans);
  
    /// (*this) = alpha * op(A) * B + beta * (*this), where A is sparse.
    /// Multiplication of sparse with dense matrix.  See also AddMatSmat.
    void AddSmatMat(Real alpha, const SparseMatrix<Real> &A,
                    MatrixTransposeType transA, const MatrixBase<Real> &B,
                    Real beta);
  
    /// (*this) = alpha * A * op(B) + beta * (*this), where B is sparse
    /// and op(B) is either B or trans(B) depending on the 'transB' argument.
    /// This is multiplication of a dense by a sparse matrix.  See also
    /// AddSmatMat.
    void AddMatSmat(Real alpha, const MatrixBase<Real> &A,
                    const SparseMatrix<Real> &B, MatrixTransposeType transB,
                    Real beta);
  
    /// *this = beta * *this + alpha * M M^T, for symmetric matrices.  It only
    /// updates the lower triangle of *this.  It will leave the matrix asymmetric;
    /// if you need it symmetric as a regular matrix, do CopyLowerToUpper().
    void SymAddMat2(const Real alpha, const MatrixBase<Real> &M,
                    MatrixTransposeType transA, Real beta);
  
    /// *this = beta * *this + alpha * diag(v) * M [or M^T].
    /// The same as adding M but scaling each row M_i by v(i).
    void AddDiagVecMat(const Real alpha, const VectorBase<Real> &v,
                       const MatrixBase<Real> &M, MatrixTransposeType transM,
                       Real beta = 1.0);
  
    /// *this = beta * *this + alpha * M [or M^T] * diag(v)
    /// The same as adding M but scaling each column M_j by v(j).
    void AddMatDiagVec(const Real alpha,
                       const MatrixBase<Real> &M, MatrixTransposeType transM,
                       VectorBase<Real> &v,
                       Real beta = 1.0);
  
    /// *this = beta * *this + alpha * A .* B (.* element by element multiplication)
    void AddMatMatElements(const Real alpha,
                           const MatrixBase<Real>& A,
                           const MatrixBase<Real>& B,
                           const Real beta);
  
    /// *this += alpha * S
    template<typename OtherReal>
    void AddSp(const Real alpha, const SpMatrix<OtherReal> &S);
  
    void AddMatMat(const Real alpha,
                   const MatrixBase<Real>& A, MatrixTransposeType transA,
                   const MatrixBase<Real>& B, MatrixTransposeType transB,
                   const Real beta);
  
    /// *this = a * b / c (by element; when c = 0, *this = a)
    void SetMatMatDivMat(const MatrixBase<Real>& A,
                         const MatrixBase<Real>& B,
                         const MatrixBase<Real>& C);
  
    /// A version of AddMatMat specialized for when the second argument
    /// contains a lot of zeroes.
    void AddMatSmat(const Real alpha,
                    const MatrixBase<Real>& A, MatrixTransposeType transA,
                    const MatrixBase<Real>& B, MatrixTransposeType transB,
                    const Real beta);
  
    /// A version of AddMatMat specialized for when the first argument
    /// contains a lot of zeroes.
    void AddSmatMat(const Real alpha,
                    const MatrixBase<Real>& A, MatrixTransposeType transA,
                    const MatrixBase<Real>& B, MatrixTransposeType transB,
                    const Real beta);
  
    /// this <-- beta*this + alpha*A*B*C.
    void AddMatMatMat(const Real alpha,
                      const MatrixBase<Real>& A, MatrixTransposeType transA,
                      const MatrixBase<Real>& B, MatrixTransposeType transB,
                      const MatrixBase<Real>& C, MatrixTransposeType transC,
                      const Real beta);
  
    /// this <-- beta*this + alpha*SpA*B.
    // This and the routines below are really
    // stubs that need to be made more efficient.
    void AddSpMat(const Real alpha,
                  const SpMatrix<Real>& A,
                  const MatrixBase<Real>& B, MatrixTransposeType transB,
                  const Real beta) {
      Matrix<Real> M(A);
      return AddMatMat(alpha, M, kNoTrans, B, transB, beta);
    }
    /// this <-- beta*this + alpha*A*B.
    void AddTpMat(const Real alpha,
                  const TpMatrix<Real>& A, MatrixTransposeType transA,
                  const MatrixBase<Real>& B, MatrixTransposeType transB,
                  const Real beta) {
      Matrix<Real> M(A);
      return AddMatMat(alpha, M, transA, B, transB, beta);
    }
    /// this <-- beta*this + alpha*A*B.
    void AddMatSp(const Real alpha,
                  const MatrixBase<Real>& A, MatrixTransposeType transA,
                  const SpMatrix<Real>& B,
                  const Real beta) {
      Matrix<Real> M(B);
      return AddMatMat(alpha, A, transA, M, kNoTrans, beta);
    }
    /// this <-- beta*this + alpha*A*B*C.
    void AddSpMatSp(const Real alpha,
                    const SpMatrix<Real> &A,
                    const MatrixBase<Real>& B, MatrixTransposeType transB,
                    const SpMatrix<Real>& C,
                  const Real beta) {
      Matrix<Real> M(A), N(C);
      return AddMatMatMat(alpha, M, kNoTrans, B, transB, N, kNoTrans, beta);
    }
    /// this <-- beta*this + alpha*A*B.
    void AddMatTp(const Real alpha,
                  const MatrixBase<Real>& A, MatrixTransposeType transA,
                  const TpMatrix<Real>& B, MatrixTransposeType transB,
                  const Real beta) {
      Matrix<Real> M(B);
      return AddMatMat(alpha, A, transA, M, transB, beta);
    }
  
    /// this <-- beta*this + alpha*A*B.
    void AddTpTp(const Real alpha,
                 const TpMatrix<Real>& A, MatrixTransposeType transA,
                 const TpMatrix<Real>& B, MatrixTransposeType transB,
                 const Real beta) {
      Matrix<Real> M(A), N(B);
      return AddMatMat(alpha, M, transA, N, transB, beta);
    }
  
    /// this <-- beta*this + alpha*A*B.
    // This one is more efficient, not like the others above.
    void AddSpSp(const Real alpha,
                 const SpMatrix<Real>& A, const SpMatrix<Real>& B,
                 const Real beta);
  
    /// Copy lower triangle to upper triangle (symmetrize)
    void CopyLowerToUpper();
  
    /// Copy upper triangle to lower triangle (symmetrize)
    void CopyUpperToLower();
  
    /// This function orthogonalizes the rows of a matrix using the Gram-Schmidt
    /// process.  It is only applicable if NumRows() <= NumCols().  It will use
    /// random number generation to fill in rows with something nonzero, in cases
    /// where the original matrix was of deficient row rank.
    void OrthogonalizeRows();
  
    /// stream read.
    /// Use instead of stream<<*this, if you want to add to existing contents.
    // Will throw exception on failure.
    void Read(std::istream & in, bool binary, bool add = false);
    /// write to stream.
    void Write(std::ostream & out, bool binary) const;
  
    // Below is internal methods for Svd, user does not have to know about this.
  #if !defined(HAVE_ATLAS) && !defined(USE_KALDI_SVD)
    // protected:
    // Should be protected but used directly in testing routine.
    // destroys *this!
    void LapackGesvd(VectorBase<Real> *s, MatrixBase<Real> *U,
                       MatrixBase<Real> *Vt);
  #else
   protected:
    // destroys *this!
    bool JamaSvd(VectorBase<Real> *s, MatrixBase<Real> *U,
                 MatrixBase<Real> *V);
  
  #endif
   protected:
  
    ///  Initializer, callable only from child.
    explicit MatrixBase(Real *data, MatrixIndexT cols, MatrixIndexT rows, MatrixIndexT stride) :
      data_(data), num_cols_(cols), num_rows_(rows), stride_(stride) {
      KALDI_ASSERT_IS_FLOATING_TYPE(Real);
    }
  
    ///  Initializer, callable only from child.
    /// Empty initializer, for un-initialized matrix.
    explicit MatrixBase(): data_(NULL) {
      KALDI_ASSERT_IS_FLOATING_TYPE(Real);
    }
  
    // Make sure pointers to MatrixBase cannot be deleted.
    ~MatrixBase() { }
  
    /// A workaround that allows SubMatrix to get a pointer to non-const data
    /// for const Matrix. Unfortunately C++ does not allow us to declare a
    /// "public const" inheritance or anything like that, so it would require
    /// a lot of work to make the SubMatrix class totally const-correct--
    /// we would have to override many of the Matrix functions.
    inline Real*  Data_workaround() const {
      return data_;
    }
  
    /// data memory area
    Real*   data_;
  
    /// these attributes store the real matrix size as it is stored in memory
    /// including memalignment
    MatrixIndexT    num_cols_;   /// < Number of columns
    MatrixIndexT    num_rows_;   /// < Number of rows
    /** True number of columns for the internal matrix. This number may differ
     * from num_cols_ as memory alignment might be used. */
    MatrixIndexT    stride_;
   private:
    KALDI_DISALLOW_COPY_AND_ASSIGN(MatrixBase);
  };
  
  /// A class for storing matrices.
  template<typename Real>
  class Matrix : public MatrixBase<Real> {
   public:
  
    /// Empty constructor.
    Matrix();
  
    /// Basic constructor.
    Matrix(const MatrixIndexT r, const MatrixIndexT c,
           MatrixResizeType resize_type = kSetZero,
           MatrixStrideType stride_type = kDefaultStride):
        MatrixBase<Real>() { Resize(r, c, resize_type, stride_type); }
  
    /// Copy constructor from CUDA matrix
    /// This is defined in ../cudamatrix/cu-matrix.h
    template<typename OtherReal>
    explicit Matrix(const CuMatrixBase<OtherReal> &cu,
                    MatrixTransposeType trans = kNoTrans);
  
  
    /// Swaps the contents of *this and *other.  Shallow swap.
    void Swap(Matrix<Real> *other);
  
    /// Defined in ../cudamatrix/cu-matrix.cc
    void Swap(CuMatrix<Real> *mat);
  
    /// Constructor from any MatrixBase. Can also copy with transpose.
    /// Allocates new memory.
    explicit Matrix(const MatrixBase<Real> & M,
                    MatrixTransposeType trans = kNoTrans);
  
    /// Same as above, but need to avoid default copy constructor.
    Matrix(const Matrix<Real> & M);  //  (cannot make explicit)
  
    /// Copy constructor: as above, but from another type.
    template<typename OtherReal>
    explicit Matrix(const MatrixBase<OtherReal> & M,
                      MatrixTransposeType trans = kNoTrans);
  
    /// Copy constructor taking SpMatrix...
    /// It is symmetric, so no option for transpose, and NumRows == Cols
    template<typename OtherReal>
    explicit Matrix(const SpMatrix<OtherReal> & M) : MatrixBase<Real>() {
      Resize(M.NumRows(), M.NumRows(), kUndefined);
      this->CopyFromSp(M);
    }
  
    /// Constructor from CompressedMatrix
    explicit Matrix(const CompressedMatrix &C);
  
    /// Copy constructor taking TpMatrix...
    template <typename OtherReal>
    explicit Matrix(const TpMatrix<OtherReal> & M,
                    MatrixTransposeType trans = kNoTrans) : MatrixBase<Real>() {
      if (trans == kNoTrans) {
        Resize(M.NumRows(), M.NumCols(), kUndefined);
        this->CopyFromTp(M);
      } else {
        Resize(M.NumCols(), M.NumRows(), kUndefined);
        this->CopyFromTp(M, kTrans);
      }
    }
  
    /// read from stream.
    // Unlike one in base, allows resizing.
    void Read(std::istream & in, bool binary, bool add = false);
  
    /// Remove a specified row.
    void RemoveRow(MatrixIndexT i);
  
    /// Transpose the matrix.  Works for non-square
    /// matrices as well as square ones.
    void Transpose();
  
    /// Distructor to free matrices.
    ~Matrix() { Destroy(); }
  
    /// Sets matrix to a specified size (zero is OK as long as both r and c are
    /// zero).  The value of the new data depends on resize_type:
    ///   -if kSetZero, the new data will be zero
    ///   -if kUndefined, the new data will be undefined
    ///   -if kCopyData, the new data will be the same as the old data in any
    ///      shared positions, and zero elsewhere.
    ///
    /// You can set stride_type to kStrideEqualNumCols to force the stride
    /// to equal the number of columns; by default it is set so that the stride
    /// in bytes is a multiple of 16.
    ///
    /// This function takes time proportional to the number of data elements.
    void Resize(const MatrixIndexT r,
                const MatrixIndexT c,
                MatrixResizeType resize_type = kSetZero,
                MatrixStrideType stride_type = kDefaultStride);
  
    /// Assignment operator that takes MatrixBase.
    Matrix<Real> &operator = (const MatrixBase<Real> &other) {
      if (MatrixBase<Real>::NumRows() != other.NumRows() ||
          MatrixBase<Real>::NumCols() != other.NumCols())
        Resize(other.NumRows(), other.NumCols(), kUndefined);
      MatrixBase<Real>::CopyFromMat(other);
      return *this;
    }
  
    /// Assignment operator. Needed for inclusion in std::vector.
    Matrix<Real> &operator = (const Matrix<Real> &other) {
      if (MatrixBase<Real>::NumRows() != other.NumRows() ||
          MatrixBase<Real>::NumCols() != other.NumCols())
        Resize(other.NumRows(), other.NumCols(), kUndefined);
      MatrixBase<Real>::CopyFromMat(other);
      return *this;
    }
  
  
   private:
    /// Deallocates memory and sets to empty matrix (dimension 0, 0).
    void Destroy();
  
    /// Init assumes the current class contents are invalid (i.e. junk or have
    /// already been freed), and it sets the matrix to newly allocated memory with
    /// the specified number of rows and columns.  r == c == 0 is acceptable.  The data
    /// memory contents will be undefined.
    void Init(const MatrixIndexT r,
              const MatrixIndexT c,
              const MatrixStrideType stride_type);
  
  };
  /// @} end "addtogroup matrix_group"
  
  /// \addtogroup matrix_funcs_io
  /// @{
  
  /// A structure containing the HTK header.
  /// [TODO: change the style of the variables to Kaldi-compliant]
  struct HtkHeader {
    /// Number of samples.
    int32    mNSamples;
    /// Sample period.
    int32    mSamplePeriod;
    /// Sample size
    int16    mSampleSize;
    /// Sample kind.
    uint16   mSampleKind;
  };
  
  // Read HTK formatted features from file into matrix.
  template<typename Real>
  bool ReadHtk(std::istream &is, Matrix<Real> *M, HtkHeader *header_ptr);
  
  // Write (HTK format) features to file from matrix.
  template<typename Real>
  bool WriteHtk(std::ostream &os, const MatrixBase<Real> &M, HtkHeader htk_hdr);
  
  // Write (CMUSphinx format) features to file from matrix.
  template<typename Real>
  bool WriteSphinx(std::ostream &os, const MatrixBase<Real> &M);
  
  /// @} end of "addtogroup matrix_funcs_io"
  
  /**
    Sub-matrix representation.
    Can work with sub-parts of a matrix using this class.
    Note that SubMatrix is not very const-correct-- it allows you to
    change the contents of a const Matrix.  Be careful!
  */
  
  template<typename Real>
  class SubMatrix : public MatrixBase<Real> {
   public:
    // Initialize a SubMatrix from part of a matrix; this is
    // a bit like A(b:c, d:e) in Matlab.
    // This initializer is against the proper semantics of "const", since
    // SubMatrix can change its contents.  It would be hard to implement
    // a "const-safe" version of this class.
    SubMatrix(const MatrixBase<Real>& T,
              const MatrixIndexT ro,  // row offset, 0 < ro < NumRows()
              const MatrixIndexT r,   // number of rows, r > 0
              const MatrixIndexT co,  // column offset, 0 < co < NumCols()
              const MatrixIndexT c);   // number of columns, c > 0
  
    // This initializer is mostly intended for use in CuMatrix and related
    // classes.  Be careful!
    SubMatrix(Real *data,
              MatrixIndexT num_rows,
              MatrixIndexT num_cols,
              MatrixIndexT stride);
  
    ~SubMatrix<Real>() {}
  
    /// This type of constructor is needed for Range() to work [in Matrix base
    /// class]. Cannot make it explicit.
    SubMatrix<Real> (const SubMatrix &other):
    MatrixBase<Real> (other.data_, other.num_cols_, other.num_rows_,
                      other.stride_) {}
  
   private:
    /// Disallow assignment.
    SubMatrix<Real> &operator = (const SubMatrix<Real> &other);
  };
  /// @} End of "addtogroup matrix_funcs_io".
  
  /// \addtogroup matrix_funcs_scalar
  /// @{
  
  // Some declarations.  These are traces of products.
  
  
  template<typename Real>
  bool ApproxEqual(const MatrixBase<Real> &A,
                   const MatrixBase<Real> &B, Real tol = 0.01) {
    return A.ApproxEqual(B, tol);
  }
  
  template<typename Real>
  inline void AssertEqual(const MatrixBase<Real> &A, const MatrixBase<Real> &B,
                          float tol = 0.01) {
    KALDI_ASSERT(A.ApproxEqual(B, tol));
  }
  
  /// Returns trace of matrix.
  template <typename Real>
  double TraceMat(const MatrixBase<Real> &A) { return A.Trace(); }
  
  
  /// Returns tr(A B C)
  template <typename Real>
  Real TraceMatMatMat(const MatrixBase<Real> &A, MatrixTransposeType transA,
                        const MatrixBase<Real> &B, MatrixTransposeType transB,
                        const MatrixBase<Real> &C, MatrixTransposeType transC);
  
  /// Returns tr(A B C D)
  template <typename Real>
  Real TraceMatMatMatMat(const MatrixBase<Real> &A, MatrixTransposeType transA,
                           const MatrixBase<Real> &B, MatrixTransposeType transB,
                           const MatrixBase<Real> &C, MatrixTransposeType transC,
                           const MatrixBase<Real> &D, MatrixTransposeType transD);
  
  /// @} end "addtogroup matrix_funcs_scalar"
  
  
  /// \addtogroup matrix_funcs_misc
  /// @{
  
  
  /// Function to ensure that SVD is sorted.  This function is made as generic as
  /// possible, to be applicable to other types of problems.  s->Dim() should be
  /// the same as U->NumCols(), and we sort s from greatest to least absolute
  /// value (if sort_on_absolute_value == true) or greatest to least value
  /// otherwise, moving the columns of U, if it exists, and the rows of Vt, if it
  /// exists, around in the same way.  Note: the "absolute value" part won't matter
  /// if this is an actual SVD, since singular values are non-negative.
  template<typename Real> void SortSvd(VectorBase<Real> *s, MatrixBase<Real> *U,
                                       MatrixBase<Real>* Vt = NULL,
                                       bool sort_on_absolute_value = true);
  
  /// Creates the eigenvalue matrix D that is part of the decomposition used Matrix::Eig.
  /// D will be block-diagonal with blocks of size 1 (for real eigenvalues) or 2x2
  /// for complex pairs.  If a complex pair is lambda +- i*mu, D will have a corresponding
  /// 2x2 block [lambda, mu; -mu, lambda].
  /// This function will throw if any complex eigenvalues are not in complex conjugate
  /// pairs (or the members of such pairs are not consecutively numbered).
  template<typename Real>
  void CreateEigenvalueMatrix(const VectorBase<Real> &real, const VectorBase<Real> &imag,
                              MatrixBase<Real> *D);
  
  /// The following function is used in Matrix::Power, and separately tested, so we
  /// declare it here mainly for the testing code to see.  It takes a complex value to
  /// a power using a method that will work for noninteger powers (but will fail if the
  /// complex value is real and negative).
  template<typename Real>
  bool AttemptComplexPower(Real *x_re, Real *x_im, Real power);
  
  
  
  /// @} end of addtogroup matrix_funcs_misc
  
  /// \addtogroup matrix_funcs_io
  /// @{
  template<typename Real>
  std::ostream & operator << (std::ostream & Out, const MatrixBase<Real> & M);
  
  template<typename Real>
  std::istream & operator >> (std::istream & In, MatrixBase<Real> & M);
  
  // The Matrix read allows resizing, so we override the MatrixBase one.
  template<typename Real>
  std::istream & operator >> (std::istream & In, Matrix<Real> & M);
  
  
  template<typename Real>
  bool SameDim(const MatrixBase<Real> &M, const MatrixBase<Real> &N) {
    return (M.NumRows() == N.NumRows() && M.NumCols() == N.NumCols());
  }
  
  /// @} end of \addtogroup matrix_funcs_io
  
  
  }  // namespace kaldi
  
  
  
  // we need to include the implementation and some
  // template specializations.
  #include "matrix/kaldi-matrix-inl.h"
  
  
  #endif  // KALDI_MATRIX_KALDI_MATRIX_H_