cu-matrix.h

Go to the documentation of this file.

// cudamatrix/cu-matrix.h

// Copyright 2009-2012  Karel Vesely
//                2013  Johns Hopkins University (author: Daniel Povey)
//                2013  Hainan Xu
//                2013  Xiaohui Zhang
//           2013-2015  Guoguo Chen
//                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_CUDAMATRIX_CU_MATRIX_H_
#define KALDI_CUDAMATRIX_CU_MATRIX_H_

#include <sstream>
#include <vector>

#include "cudamatrix/cu-matrixdim.h"
#include "cudamatrix/cu-common.h"
#include "cudamatrix/cu-value.h"
#include "matrix/matrix-common.h"
#include "matrix/kaldi-matrix.h"
#include "cudamatrix/cu-array.h"
#include "cudamatrix/cu-math.h"
#include "cudamatrix/cu-rand.h"
#include "cudamatrix/cu-sparse-matrix.h"

namespace kaldi {

template<typename Real>
Real TraceMatMat(const CuMatrixBase<Real> &A, const CuMatrixBase<Real> &B,
                 MatrixTransposeType trans = kNoTrans);

template<typename Real>
void AddMatMatBatched(const Real alpha, std::vector<CuSubMatrix<Real>* > &C,
                      const std::vector<CuSubMatrix<Real>* > &A,
                      MatrixTransposeType transA,
                      const std::vector<CuSubMatrix<Real>* > &B,
                      MatrixTransposeType transB,
                      const Real beta);

/*
template<typename Real>
struct MatrixElement {
  int row;
  int column;
  Real weight;
};
// */

template<typename Real>
class CuMatrixBase {
 public:
  friend class CuMatrixBase<float>;
  friend class CuMatrixBase<double>;
  friend class CuVectorBase<float>;
  friend class CuVectorBase<double>;
  friend class VectorBase<Real>;
  friend class CuSpMatrix<Real>;
  friend class CuTpMatrix<float>;
  friend class CuTpMatrix<double>;
  friend class CuVectorBase<Real>;
  friend class CuSubMatrix<Real>;
  friend class CuRand<Real>;
  friend class CuSubVector<Real>;
  friend class CuBlockMatrix<Real>;
  friend class CuSparseMatrix<float>;
  friend class CuSparseMatrix<double>;
  friend class CuSparseMatrix<Real>;

  void CopyCols(const CuMatrixBase<Real> &src,
                const CuArrayBase<MatrixIndexT> &indexes);


  void AddCols(const CuMatrixBase<Real> &src,
               const CuArrayBase<MatrixIndexT> &indices);

  void CopyRows(const CuMatrixBase<Real> &src,
                const CuArrayBase<MatrixIndexT> &indexes);

  void CopyRows(const CuArrayBase<const Real*> &src);

  void CopyToRows(const CuArrayBase<Real*> &dst) const;

  void AddRows(Real alpha,
               const CuMatrixBase<Real> &src,
               const CuArrayBase<MatrixIndexT> &indexes);


  void MulRows(const CuMatrixBase<Real> &src,
               const CuArrayBase<MatrixIndexT> &indexes);


  void AddRows(Real alpha,
               const CuArrayBase<const Real*> &src);


  void AddToRows(Real alpha,
                 const CuArrayBase<MatrixIndexT> &indexes,
                 CuMatrixBase<Real> *dst) const;


  void AddToRows(Real alpha, const CuArrayBase<Real*> &dst) const;


  void SumColumnRanges(const CuMatrixBase<Real> &src,
                       const CuArrayBase<Int32Pair> &indexes);


  void AddRowRanges(const CuMatrixBase<Real> &src,
                    const CuArrayBase<Int32Pair> &indexes);


  friend Real TraceMatMat<Real>(const CuMatrixBase<Real> &A,
                                const CuMatrixBase<Real> &B,
                                MatrixTransposeType trans);

  friend Real TraceMatSmat<Real>(const CuMatrixBase<Real> &A,
                                 const CuSparseMatrix<Real> &B,
                                 MatrixTransposeType trans);

  friend void AddMatMatBatched<Real>(const Real alpha,
                                     std::vector<CuSubMatrix<Real>* > &C,
                                     const std::vector<CuSubMatrix<Real>* > &A,
                                     MatrixTransposeType transA,
                                     const std::vector<CuSubMatrix<Real>* > &B,
                                     MatrixTransposeType transB,
                                     const Real beta);

  void AddToDiag(Real value);

  MatrixIndexT NumRows() const { return num_rows_;  }
  MatrixIndexT NumCols() const { return num_cols_;  }
  MatrixIndexT Stride() const { return stride_; }

  // MatrixDim is a struct containing "rows", "cols" and "stride",
  // that is an argument of most CUDA kernels.
  ::MatrixDim Dim() const {
    ::MatrixDim d = { num_rows_, num_cols_, stride_ };
    return d;
  }

  Real FrobeniusNorm() const { return sqrt(TraceMatMat(*this, *this, kTrans)); }

  bool IsUnit(Real tol = 0.001) const;

  bool ApproxEqual(const CuMatrixBase<Real> &other, float tol = 0.01) const;

  MatrixIndexT SizeInBytes() const { return num_rows_*stride_*sizeof(Real); }

  // Copy functions.  These do not resize.
  template<typename OtherReal>
  void CopyFromMat(const MatrixBase<OtherReal> &src,
                   MatrixTransposeType trans = kNoTrans);


  void CopyFromGeneralMat(const GeneralMatrix &src,
                          MatrixTransposeType trans = kNoTrans);

  void CopyFromMat(const MatrixBase<Real> &src,
                   MatrixTransposeType trans = kNoTrans);

  void CopyFromSp(const CuSpMatrix<Real> &M);

  template<typename OtherReal>
  void CopyFromTp(const CuTpMatrix<OtherReal> &M,
                  MatrixTransposeType trans = kNoTrans);
  
  // This function will copy from source rows (start_range, end_range]
  // if the range is outside of the clamped region then the clamped
  // row will be replicated across the out of range areas
  void CopyRangeFromMatClamped(const CuMatrixBase<Real> & src,
      int32_t start_range, int32_t end_range,
      int32_t clamp_low, int32_t clamp_high);

  template<typename OtherReal>
  void CopyFromMat(const CuMatrixBase<OtherReal> &M,
                   MatrixTransposeType trans = kNoTrans);

  template<typename OtherReal>
  void CopyToMat(MatrixBase<OtherReal> *dst,
                 MatrixTransposeType trans = kNoTrans) const;

  void CopyRowsFromVec(const CuVectorBase<Real> &v);

  void CopyRowsFromVec(const VectorBase<Real> &v);

  void CopyColsFromVec(const CuVectorBase<Real> &v);

  void CopyColFromVec(const CuVectorBase<Real> &v, const MatrixIndexT col);

  void Sigmoid(const CuMatrixBase<Real> &src);

  void Heaviside(const CuMatrixBase<Real> &src);

  void Exp(const CuMatrixBase<Real> &src);

  void Log(const CuMatrixBase<Real> &src);

  void Pow(const CuMatrixBase<Real> &src, Real power);

  void PowAbs(const CuMatrixBase<Real> &src, Real power, bool include_sign=false);

  void Floor(const CuMatrixBase<Real> &src, Real floor_val);
  
  void Ceiling(const CuMatrixBase<Real> &src, Real ceiling_val);
  
  void ExpLimited(const CuMatrixBase<Real> &src, Real lower_limit, Real upper_limit);

  void ExpSpecial(const CuMatrixBase<Real> &src);
  
  void SoftMaxPerRow(const CuMatrixBase<Real> &src);

  void LogSoftMaxPerRow(const CuMatrixBase<Real> &src);

  
  void SoftHinge(const CuMatrixBase<Real> &src);

  void GroupPnorm(const CuMatrixBase<Real> &src, Real pow);

  void DiffGroupPnorm(const CuMatrixBase<Real> &in_value,
                      const CuMatrixBase<Real> &out_value,
                      const CuMatrixBase<Real> &out_deriv, Real power);

  void GroupMax(const CuMatrixBase<Real> &src);

  void GroupMaxDeriv(const CuMatrixBase<Real> &input,
                     const CuMatrixBase<Real> &output);

  void ParametricRelu(const CuMatrixBase<Real> &src,
                      const CuVectorBase<Real> &alpha,
                      const CuVectorBase<Real> &beta);

  void DiffParametricRelu(const CuMatrixBase<Real> &value,
                          const CuMatrixBase<Real> &diff,
                          const CuVectorBase<Real> &alpha,
                          const CuVectorBase<Real> &beta);

  void Tanh(const CuMatrixBase<Real> &src);

  void DiffSigmoid(const CuMatrixBase<Real> &value,
                   const CuMatrixBase<Real> &diff);

  void DiffTanh(const CuMatrixBase<Real> &value,
                const CuMatrixBase<Real> &diff);

  void DiffSoftmaxPerRow(const CuMatrixBase<Real> &value,
                         const CuMatrixBase<Real> &diff);

  void DiffLogSoftmaxPerRow(const CuMatrixBase<Real> &out_value,
                            const CuMatrixBase<Real> &out_deriv);

  void DiffXent(const CuArrayBase<int32> &tgt,
                CuVector<Real> *log_post_tgt);

  void Cholesky(CuMatrixBase<Real> *inv_cholesky = NULL);


  void SymInvertPosDef();

  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 ApplyExpLimited(Real lower_limit, Real upper_limit) {
    this -> ExpLimited(*this, lower_limit, upper_limit);
  };

  inline void ApplyExpSpecial() {
    this -> ExpSpecial(*this);
  };

  inline void ApplySoftMaxPerRow() {
    this -> SoftMaxPerRow(*this);
  };

  inline void ApplyLogSoftMaxPerRow() {
    this -> LogSoftMaxPerRow(*this);
  };

  inline void ApplyLog() {
    this -> Log(*this);
  };

  void FindRowMaxId(CuArray<int32> *id) const;

  void SetZero();
  void Set(Real value);
  void Add(Real value);
  void SetZeroAboveDiag();
  void Scale(Real value);

  void MulElements(const CuMatrixBase<Real> &A);
  void DivElements(const CuMatrixBase<Real> &A);
  void Max(const CuMatrixBase<Real> &A);
  void Min(const CuMatrixBase<Real> &A);
  void MulColsVec(const CuVectorBase<Real> &scale);
  void MulRowsVec(const CuVectorBase<Real> &scale);
  void MulRowsGroupMat(const CuMatrixBase<Real> &src);
  void DivRowsVec(const CuVectorBase<Real> &div);
  void InvertElements();
  void AddMat(Real alpha, const CuMatrixBase<Real> &A,
              MatrixTransposeType trans = kNoTrans);

  void AddSmat(Real alpha, const CuSparseMatrix<Real> &A,
              MatrixTransposeType trans = kNoTrans);

  void AddSmatMat(Real alpha, const CuSparseMatrix<Real> &A,
                  MatrixTransposeType transA, const CuMatrixBase<Real> &B,
                  Real beta);

  void AddMatSmat(Real alpha, const CuMatrixBase<Real> &A,
                  const CuSparseMatrix<Real> &B, MatrixTransposeType transB,
                  Real beta);


  void AddToElements(Real alpha, const CuArrayBase<int32> &elements);


  void AddMatBlocks(Real alpha, const CuMatrixBase<Real> &A,
                    MatrixTransposeType trans = kNoTrans);

  void AddVecToCols(Real alpha, const CuVectorBase<Real> &col, Real beta = 1.0);
  void AddVecToRows(Real alpha, const CuVectorBase<Real> &row, Real beta = 1.0);
  void AddMatMat(Real alpha, const CuMatrixBase<Real> &A, MatrixTransposeType transA,
                 const CuMatrixBase<Real> &B, MatrixTransposeType transB, Real beta);
  void AddVecVec(Real alpha, const CuVectorBase<Real> &x, const CuVectorBase<Real> &y);
  void SetMatMatDivMat(const CuMatrixBase<Real> &A, const CuMatrixBase<Real> &B, const CuMatrixBase<Real> &C);

  void SymAddMat2(const Real alpha, const CuMatrixBase<Real> &M,
                  MatrixTransposeType transA, Real beta);


  void AddMatBlock(Real alpha, const CuMatrixBase<Real> &A, MatrixTransposeType transA,
                   const CuBlockMatrix<Real> &B, MatrixTransposeType transB, Real beta);

  void AddDiagVecMat(const Real alpha, const CuVectorBase<Real> &v,
                     const CuMatrixBase<Real> &M, MatrixTransposeType transM,
                     Real beta = 1.0);

  // *this = beta * *this + alpha * M  * diag(v) [or M^T].
  // The same as adding M but scaling each column M_j by v(j).
  void AddMatDiagVec(const Real alpha,
                     const CuMatrixBase<Real> &M, MatrixTransposeType transM,
                     CuVectorBase<Real> &v,
                     Real beta = 1.0);

  void AddMatMatElements(const Real alpha,
                         const CuMatrixBase<Real>& A,
                         const CuMatrixBase<Real>& B,
                         const Real beta);

  void AddMatSp(const Real alpha,
                const CuMatrixBase<Real> &A, MatrixTransposeType transA,
                const CuSpMatrix<Real> &B,
                const Real beta) {
    CuMatrix<Real> M(B);
    return AddMatMat(alpha, A, transA, M, kNoTrans, beta);
  }

  void AddSpMat(const Real alpha,
                const CuSpMatrix<Real> &A,
                const CuMatrixBase<Real> &B, MatrixTransposeType transB,
                const Real beta) {
    CuMatrix<Real> M(A);
    return AddMatMat(alpha, M, kNoTrans, B, transB, beta);
  }

  void AddTpMat(const Real alpha,
                const CuTpMatrix<Real> &A, MatrixTransposeType transA,
                const CuMatrixBase<Real> &B, MatrixTransposeType transB,
                const Real beta) {
    CuMatrix<Real> M(A);
    return AddMatMat(alpha, M, transA, B, transB, beta);
  }

  void AddMatTp(const Real alpha,
                const CuMatrixBase<Real> &A, MatrixTransposeType transA,
                const CuTpMatrix<Real> &B, MatrixTransposeType transB,
                const Real beta) {
    CuMatrix<Real> M(B);
    return AddMatMat(alpha, A, transA, M, transB, beta);
  }

  void CopyFromBlock(const CuBlockMatrix<Real> &B,
                     MatrixTransposeType trans = kNoTrans);
  void CopyLowerToUpper();
  void CopyUpperToLower();
  inline CuSubMatrix<Real> Range(const MatrixIndexT row_offset,
                                 const MatrixIndexT num_rows,
                                 const MatrixIndexT col_offset,
                                 const MatrixIndexT num_cols) const {
    return CuSubMatrix<Real>(*this, row_offset, num_rows,
                             col_offset, num_cols);
  }
  inline CuSubMatrix<Real> RowRange(const MatrixIndexT row_offset,
                                    const MatrixIndexT num_rows) const {
    return CuSubMatrix<Real>(*this, row_offset, num_rows,
                             0, num_cols_);
  }
  inline CuSubMatrix<Real> ColRange(const MatrixIndexT col_offset,
                                    const MatrixIndexT num_cols) const {
    return CuSubMatrix<Real>(*this, 0, num_rows_, col_offset, num_cols);
  }

  inline const CuSubVector<Real> Row(MatrixIndexT i) const {
    KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(i) <
                 static_cast<UnsignedMatrixIndexT>(num_rows_));
    return CuSubVector<Real>(data_ + (i * stride_), NumCols());
  }

  inline CuSubVector<Real> Row(MatrixIndexT i) {
    KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(i) <
                 static_cast<UnsignedMatrixIndexT>(num_rows_));
    return CuSubVector<Real>(data_ + (i * stride_), NumCols());
  }

  inline CuValue<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 CuValue<Real>(data_ + r * stride_ + c);
  }

  inline 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 CuValue<Real>(data_ + r * stride_ + c);  // will be casted to Real.
  }

  Real Sum() const;
  Real Max() const;
  Real Min() const;

  Real Trace(bool check_square = true) const;

  void SetRandn();

  void SetRandUniform();

  void Write(std::ostream &os, bool binary) const;

  // This function, adds a list of MatrixElements (scaled by alpha) to corresponding locations to
  // (*this).
  void AddElements(Real alpha, const std::vector<MatrixElement<Real> >& input);

  // For each i, with indexes[i] = (j, k), does (*this)(j, k) += input[i].
  // Requires, but does not check, that the vector of indexes does not contrain
  // repeated elements, 'input' is the start of an array of length equal to
  // indexes.Dim(), which is located on GPU memory if we are using the GPU.
  void AddElements(Real alpha, const CuArrayBase<Int32Pair> &indexes,
                   const Real *input);

  // This function requires that 'output' is a host array and is allocated with size
  // of indexes.size(), and for each element of 'indexes' it interprets it as
  // a (row, column) index into *this, and puts (*this)(row, column) into
  // the corresponding element of 'output'.
  void Lookup(const std::vector<Int32Pair> &indexes,
              Real *output) const;

  // CUDA version of Lookup, would be called internally by the above function.
  void Lookup(const CuArrayBase<Int32Pair> &indexes,
              Real *output) const;

  // Creates binary mask with per-element equality predicates of *this, mat.
  // Output stored to 'mask', values : 1.0 = equal, 0.0 = not-equal.
  void EqualElementMask(const CuMatrixBase<Real> &mat, CuMatrix<Real> *mask) const;


  inline const Real* RowData(MatrixIndexT r) const { return data_ + r * stride_; }
  inline Real* RowData(MatrixIndexT r) { return data_ + r * stride_; }
  inline const Real *Data() const { return data_; }
  inline Real *Data() { return data_; }

  // The following two functions should only be called if we did not compile
  // with CUDA or could not get a CUDA card; in that case the contents are
  // interpreted the same as a regular matrix.  DON'T USE THESE UNLESS YOU KNOW
  // WHAT YOU ARE DOING!
  inline const MatrixBase<Real> &Mat() const {
    return *(reinterpret_cast<const MatrixBase<Real>* >(this));
  }
  inline MatrixBase<Real> &Mat() {
    return *(reinterpret_cast<MatrixBase<Real>* >(this));
  }

 protected:

  // The constructors are protected to prevent the user creating an instance of
  // this class (you should create a child class CuMatrix or CuSubMatrix.

  CuMatrixBase(): data_(NULL), num_cols_(0), num_rows_(0), stride_(0) { }

  CuMatrixBase(Real *data,
               MatrixIndexT num_rows,
               MatrixIndexT num_cols,
               MatrixIndexT stride):
  data_(data), num_cols_(num_cols), num_rows_(num_rows), stride_(stride) { }

  Real *data_;       
  // Note: it might seem a bit backwards that we have the number of columns
  // first here; it's necessary because we need the data to be laid out the same
  // as for MatrixBase so the Mat() function call will work.  We don't want to
  // change the layout of MatrixBase at this point, or there will be crashes if
  // people don't thoroughly recompile.
  MatrixIndexT num_cols_;
  MatrixIndexT num_rows_;
  MatrixIndexT stride_;

 private:
  KALDI_DISALLOW_COPY_AND_ASSIGN(CuMatrixBase);
}; // class CuMatrixBase

template<typename Real>
class CuMatrix: public CuMatrixBase<Real> {
 public:

  CuMatrix() { }

  CuMatrix(MatrixIndexT rows, MatrixIndexT cols,
           MatrixResizeType resize_type = kSetZero,
           MatrixStrideType stride_type = kDefaultStride) {
    Resize(rows, cols, resize_type, stride_type);
  }

  // Note: we had to remove the "explicit" keyword due
  // to problems with STL vectors of CuMatrixBase.
  CuMatrix(const CuMatrix<Real> &other,
           MatrixTransposeType trans = kNoTrans);

  explicit CuMatrix(const CuBlockMatrix<Real> &other,
                    MatrixTransposeType trans = kNoTrans);

  explicit CuMatrix(const CuMatrixBase<Real> &other,
                    MatrixTransposeType trans = kNoTrans);

  template<typename OtherReal>
  explicit CuMatrix(const MatrixBase<OtherReal> &other,
                    MatrixTransposeType trans = kNoTrans);

  explicit CuMatrix(const CuSpMatrix<Real> &M) : CuMatrixBase<Real>() {
    Resize(M.NumRows(), M.NumRows(), kUndefined);
    this->CopyFromSp(M);
  }

  template <typename OtherReal>
  explicit CuMatrix(const CuTpMatrix<OtherReal> & M,
                    MatrixTransposeType trans = kNoTrans) : CuMatrixBase<Real>() {
    Resize(M.NumCols(), M.NumRows(), kUndefined);
    this->CopyFromTp(M, trans);
  }

  template<typename OtherReal>
  explicit CuMatrix(const CuMatrixBase<OtherReal> &M,
                    MatrixTransposeType trans = kNoTrans);

  CuMatrix<Real> &operator = (const CuMatrixBase<Real> &other) {
    this->Resize(other.NumRows(), other.NumCols(), kUndefined);
    this->CopyFromMat(other);
    return *this;
  }

  CuMatrix<Real> &operator = (const CuMatrix<Real> &other) {
    this->Resize(other.NumRows(), other.NumCols(), kUndefined);
    this->CopyFromMat(other);
    return *this;
  }

  CuMatrix<Real> &operator = (const MatrixBase<Real> &other) {
    this->Resize(other.NumRows(), other.NumCols(), kUndefined);
    this->CopyFromMat(other);
    return *this;
  }

  void Transpose();

  void Resize(MatrixIndexT rows, MatrixIndexT cols,
              MatrixResizeType resize_type = kSetZero,
              MatrixStrideType stride_type = kDefaultStride);

  void Swap(Matrix<Real> *mat);
  void Swap(CuMatrix<Real> *mat);

  template<typename OtherReal>
  void Swap(CuMatrix<OtherReal> *mat);

  void Read(std::istream &is, bool binary);

  ~CuMatrix() { Destroy(); }

  inline const Matrix<Real> &Mat() const {
    return *(reinterpret_cast<const Matrix<Real>* >(this));
  }
  inline Matrix<Real> &Mat() {
    return *(reinterpret_cast<Matrix<Real>* >(this));
  }

  void CompObjfAndDeriv(const std::vector<MatrixElement<Real> > &elements,
                        const CuMatrix<Real> &A,
                        Real *tot_objf,
                        Real *tot_weight);

 private:
  void Destroy();
};


template<typename Real>
class CuSubMatrix: public CuMatrixBase<Real> {
 public:
  inline CuSubMatrix(const CuMatrixBase<Real> &mat,
                     const MatrixIndexT row_offset,
                     const MatrixIndexT num_rows,
                     const MatrixIndexT col_offset,
                     const MatrixIndexT num_cols);

  // This constructor should be used with caution; it can be used for
  // constructing 'fake' submatrices if you want to play with
  // the stride. 'data' should point to GPU data if you're using the
  // GPU.
  inline CuSubMatrix(const Real *data,
                     const MatrixIndexT num_rows,
                     const MatrixIndexT num_cols,
                     const MatrixIndexT stride);

  inline CuSubMatrix<Real> (const CuSubMatrix &other):
  CuMatrixBase<Real> (other.data_, other.num_rows_, other.num_cols_,
                      other.stride_) {}
 private:
  CuSubMatrix<Real> &operator = (const CuSubMatrix<Real> &other);
};


template<typename Real>
bool ApproxEqual(const CuMatrixBase<Real> &A,
                 const CuMatrixBase<Real> &B, Real tol = 0.01) {
  return A.ApproxEqual(B, tol);
}

template<typename Real>
inline void AssertEqual(const CuMatrixBase<Real> &A,
                        const CuMatrixBase<Real> &B, float tol = 0.01) {
  KALDI_ASSERT(A.ApproxEqual(B, tol));
}

template<typename Real>
bool SameDim(const CuMatrixBase<Real> &M, const CuMatrixBase<Real> &N) {
  return (M.NumRows() == N.NumRows() && M.NumCols() == N.NumCols());
}

template<typename Real>
bool SameDimAndStride(const CuMatrixBase<Real> &M, const CuMatrixBase<Real> &N) {
  return (M.NumRows() == N.NumRows() && M.NumCols() == N.NumCols()
          && M.Stride() == N.Stride());
}

template<typename Real>
std::ostream &operator << (std::ostream &out, const CuMatrixBase<Real> &mat);


template<typename Real>
template<typename OtherReal>
Matrix<Real>::Matrix(const CuMatrixBase<OtherReal> &M,
                     MatrixTransposeType trans) {
  if (trans == kNoTrans) Init(M.NumRows(), M.NumCols(), kDefaultStride);
  else Init(M.NumCols(), M.NumRows(), kDefaultStride);
  M.CopyToMat(this, trans);
}

template<typename Real>
template<typename OtherReal>
void MatrixBase<Real>::CopyFromMat(const CuMatrixBase<OtherReal> &cu,
                                   MatrixTransposeType trans) {
  cu.CopyToMat(this, trans);
}


}  // namespace


#include "cudamatrix/cu-matrix-inl.h"

#endif