// cudamatrix/cu-matrix.cc

// Copyright 2009-2012 Karel Vesely, Lucas Ondel

// 2013 Ehsan Variani

// 2013 Johns Hopkins University (author: Daniel Povey)

// 2013 Hainan Xu

// 2013 Xiaohui Zhang

// 2013-2015 Guoguo Chen

// 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.

#if HAVE_CUDA == 1

#include <cuda_runtime_api.h>

#include <cublas_v2.h>

#endif

#include "base/timer.h"

#include "cudamatrix/cu-common.h"

#include "cudamatrix/cu-vector.h"

#include "cudamatrix/cu-device.h"

#include "cudamatrix/cu-kernels.h"

#include "cudamatrix/cu-randkernels.h"

#include "cudamatrix/cu-array.h"

#include "cudamatrix/cu-math.h"

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

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

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

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

#include "cudamatrix/cublas-wrappers.h"

namespace kaldi {

template<typename Real>

void CuMatrix<Real>::Resize(MatrixIndexT rows, MatrixIndexT cols,

MatrixResizeType resize_type,

MatrixStrideType stride_type) {

// This code does not currently support the other resize_type options.

KALDI_ASSERT(resize_type == kSetZero || resize_type == kUndefined);

if (rows * cols == 0) KALDI_ASSERT(rows == 0 && cols == 0);

if (this->num_rows_ == rows && this->num_cols_ == cols) {

if (resize_type == kSetZero) this->SetZero();

return;

}

if (this->num_rows_ != 0)

this->Destroy();

if (rows == 0) return;

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

Timer tim;

MatrixIndexT row_bytes = cols * sizeof(Real);

size_t pitch;

if (stride_type == kDefaultStride) {

this->data_ = static_cast<Real*>(CuDevice::Instantiate().MallocPitch(

row_bytes, rows, &pitch));

this->num_rows_ = rows;

this->num_cols_ = cols;

this->stride_ = pitch / sizeof(Real);

} else { // kStrideEqualNumCols

size_t bytes = rows * cols * sizeof(Real);

this->data_ = static_cast<Real*>(CuDevice::Instantiate().Malloc(bytes));

this->num_rows_ = rows;

this->num_cols_ = cols;

this->stride_ = cols;

}

if (resize_type == kSetZero) this->SetZero();

CuDevice::Instantiate().AccuProfile("CuMatrix::Resize", tim.Elapsed());

} else

#endif

{ // Let the initializer of Matrix<Real> handle the allocation,

// and then just do Swap which will switch the pointers.

// This wastes a few instructions but is simple to code.

Matrix<Real> mat(rows, cols, resize_type, stride_type);

this->Swap(&mat);

}

}

template<typename Real>

void CuMatrix<Real>::Destroy() {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

if (this->data_ != NULL) {

Timer tim;

CuDevice::Instantiate().Free(this->data_);

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

}

} else

#endif

{

if (this->data_ != NULL) KALDI_MEMALIGN_FREE(this->data_);

}

this->data_ = NULL;

this->num_rows_ = 0;

this->num_cols_ = 0;

this->stride_ = 0;

}

template<typename Real>

void CuMatrix<Real>::Swap(CuMatrix<Real> *mat) {

std::swap(mat->data_, this->data_);

std::swap(mat->num_cols_, this->num_cols_);

std::swap(mat->num_rows_, this->num_rows_);

std::swap(mat->stride_, this->stride_);

}

template<typename Real>

void CuMatrix<Real>::Swap(Matrix<Real> *mat) {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

if (this->num_rows_ == 0) {

if (mat->num_rows_ != 0) {

// *this is empty, but mat is nonempty.

this->Resize(mat->num_rows_, mat->num_cols_, kUndefined);

this->CopyFromMat(*mat);

mat->Resize(0, 0);

}

// else both are empty.

} else { // *this is nonempty.

if (mat->num_rows_ != 0) {

// Both *this and *mat are nonempty. Recurse to simpler cases.

// this could be done more efficiently in the case where

// the size does not change.

Matrix<Real> temp;

this->Swap(&temp); // now temp is full, *this is empty.

mat->Swap(&temp); // now mat has data from *this, temp has

// data from mat.

this->Swap(&temp); // copy data in mat to *this, which is now empty.

} else { // *this is full but *mat is empty.

mat->Resize(this->num_rows_, this->num_cols_, kUndefined);

this->CopyToMat(mat);

this->Destroy();

}

}

} else

#endif

{

std::swap(mat->data_, this->data_);

std::swap(mat->num_cols_, this->num_cols_);

std::swap(mat->num_rows_, this->num_rows_);

std::swap(mat->stride_, this->stride_);

}

}

template <typename Real>

void CuMatrixBase<Real>::CopyFromBlock(const CuBlockMatrix<Real> &B,

MatrixTransposeType trans) {

this->SetZero();

if (trans == kNoTrans) {

KALDI_ASSERT(NumRows() == B.NumRows() && NumCols() == B.NumCols());

int32 row_offset = 0, col_offset = 0;

for (int32 b = 0; b < B.NumBlocks(); b++) {

const CuMatrixBase<Real> &block = B.Block(b);

int32 num_rows = block.NumRows(), num_cols = block.NumCols();

CuSubMatrix<Real> this_block(*this, row_offset, num_rows,

col_offset, num_cols);

this_block.CopyFromMat(block);

row_offset += num_rows;

col_offset += num_cols;

}

KALDI_ASSERT(row_offset == NumRows() && col_offset == NumCols());

} else {

KALDI_ASSERT(NumRows() == B.NumCols() && NumCols() == B.NumRows());

int32 row_offset = 0, col_offset = 0;

for (int32 b = 0; b < B.NumBlocks(); b++) {

const CuMatrixBase<Real> &block = B.Block(b);

int32 num_rows = block.NumCols(), num_cols = block.NumRows();

CuSubMatrix<Real> this_block(*this, row_offset, num_rows,

col_offset, num_cols);

this_block.CopyFromMat(block, kTrans);

row_offset += num_rows;

col_offset += num_cols;

}

KALDI_ASSERT(row_offset == NumRows() && col_offset == NumCols());

}

}

template <typename Real>

CuMatrix<Real>::CuMatrix(const CuBlockMatrix<Real> &B,

MatrixTransposeType trans): CuMatrixBase<Real>() {

if (trans == kNoTrans) {

Resize(B.NumRows(), B.NumCols(), kUndefined);

this->CopyFromBlock(B);

} else {

Resize(B.NumCols(), B.NumRows(), kUndefined);

this->CopyFromBlock(B, kTrans);

}

}

template<class Real>

template<class OtherReal>

void CuMatrixBase<Real>::CopyFromMat(const CuMatrixBase<OtherReal> &M,

MatrixTransposeType trans) {

if (sizeof(Real) == sizeof(OtherReal) &&

static_cast<const void*>(M.Data()) ==

static_cast<const void*>(this->Data())) {

if (M.Data() == NULL)

return;

// CopyFromMat called on same data. Nothing to do (except sanity checks)

KALDI_ASSERT(trans == kNoTrans && M.NumRows() == NumRows() &&

M.NumCols() == NumCols() && M.Stride() == Stride());

return;

}

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

if (trans == kNoTrans) {

KALDI_ASSERT(M.NumRows() == num_rows_ && M.NumCols() == num_cols_);

} else {

KALDI_ASSERT(M.NumCols() == num_rows_ && M.NumRows() == num_cols_);

}

if (M.num_rows_ == 0) return; // Nothing to do.

Timer tim;

if (sizeof(Real) == sizeof(OtherReal) && trans == kNoTrans ) {

MatrixIndexT dst_pitch = stride_ * sizeof(Real);

MatrixIndexT src_pitch = M.Stride() * sizeof(Real);

MatrixIndexT width = M.NumCols() * sizeof(Real);

CU_SAFE_CALL(cudaMemcpy2D(data_, dst_pitch, M.data_, src_pitch,

width, M.num_rows_, cudaMemcpyDeviceToDevice));

} else {

if (trans == kNoTrans) {

dim3 dimGrid, dimBlock;

GetBlockSizesForSimpleMatrixOperation(NumRows(), NumCols(),

&dimGrid, &dimBlock);

cuda_copy_from_mat(dimGrid, dimBlock, data_, M.data_, Dim(), M.Dim());

} else {

// 2D thread block with warps (blockDim.x) along the row-dim of input M.

// Each (8x32) thread block will transpose (32x32) data

const int32 warpSize = 32;

dim3 dimBlock(warpSize, CU1DBLOCK / warpSize);

dim3 dimGrid(n_blocks(M.NumCols(), warpSize),

n_blocks(M.NumRows(), warpSize));

cuda_copy_from_mat_trans(dimGrid, dimBlock, data_, M.data_, Dim(),

M.Dim());

}

CU_SAFE_CALL(cudaGetLastError());

}

CuDevice::Instantiate().AccuProfile("CuMatrixBase::CopyFromMat(from other CuMatrixBase)", tim.Elapsed());

} else

#endif

{

Mat().CopyFromMat(M.Mat(), trans);

}

}

// Instantiate the template above.

template

void CuMatrixBase<float>::CopyFromMat<float>(const CuMatrixBase<float> &M,

MatrixTransposeType trans);

template

void CuMatrixBase<float>::CopyFromMat<double>(const CuMatrixBase<double> &M,

MatrixTransposeType trans);

template

void CuMatrixBase<double>::CopyFromMat<float>(const CuMatrixBase<float> &M,

MatrixTransposeType trans);

template

void CuMatrixBase<double>::CopyFromMat<double>(const CuMatrixBase<double> &M,

MatrixTransposeType trans);

template<typename Real>

template<typename OtherReal>

void CuMatrixBase<Real>::CopyFromTp(const CuTpMatrix<OtherReal> &M,

MatrixTransposeType trans) {

KALDI_ASSERT(num_rows_ == M.NumRows() && num_cols_ == num_rows_);

if (num_rows_ == 0)

return;

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

Timer tim;

dim3 dimBlock(CU2DBLOCK, CU2DBLOCK);

dim3 dimGrid(n_blocks(num_rows_, CU2DBLOCK),

n_blocks(num_rows_, CU2DBLOCK));

if (trans == kNoTrans) {

cuda_copy_from_tp(dimGrid, dimBlock, data_, M.Data(), Dim());

} else {

cuda_copy_from_tp_trans(dimGrid, dimBlock, data_, M.Data(), Dim());

}

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

Mat().CopyFromTp(M.Mat(), trans);

}

}

// instantiate the template above.

template void CuMatrixBase<float>::CopyFromTp(const CuTpMatrix<float> &M,

MatrixTransposeType trans);

template void CuMatrixBase<float>::CopyFromTp(const CuTpMatrix<double> &M,

MatrixTransposeType trans);

template void CuMatrixBase<double>::CopyFromTp(const CuTpMatrix<float> &M,

MatrixTransposeType trans);

template void CuMatrixBase<double>::CopyFromTp(const CuTpMatrix<double> &M,

MatrixTransposeType trans);

template<typename Real>

void CuMatrixBase<Real>::CopyFromMat(const MatrixBase<Real> &src,

MatrixTransposeType trans) {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

if (trans == kNoTrans) {

KALDI_ASSERT(src.NumRows() == num_rows_ && src.NumCols() == num_cols_);

Timer tim;

MatrixIndexT dst_pitch = stride_*sizeof(Real);

MatrixIndexT src_pitch = src.Stride()*sizeof(Real);

MatrixIndexT width = src.NumCols()*sizeof(Real);

CU_SAFE_CALL(cudaMemcpy2D(data_, dst_pitch, src.Data(), src_pitch,

width, src.NumRows(), cudaMemcpyHostToDevice));

CuDevice::Instantiate().AccuProfile("CuMatrixBase::CopyFromMat(from CPU)",tim.Elapsed());

} else {

CuMatrix<Real> trans_mat(src); // Do the transpose on the GPU board.

this->CopyFromMat(trans_mat, kTrans);

}

} else

#endif

{

Mat().CopyFromMat(src, trans);

}

}

template<typename Real>

template<typename OtherReal>

void CuMatrixBase<Real>::CopyFromMat(const MatrixBase<OtherReal> &src,

MatrixTransposeType trans) {

CuMatrix<OtherReal> temp(src);

this->CopyFromMat(temp, trans);

}

// instantiate the template above.

template

void CuMatrixBase<float>::CopyFromMat(const MatrixBase<double> &src,

MatrixTransposeType trans);

template

void CuMatrixBase<double>::CopyFromMat(const MatrixBase<float> &src,

MatrixTransposeType trans);

template<typename Real>

void CuMatrixBase<Real>::CopyFromSp(const CuSpMatrix<Real> &M) {

KALDI_ASSERT(num_rows_ == M.NumRows() && num_cols_ == num_rows_);

if (num_rows_ == 0)

return;

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

Timer tim;

dim3 dimBlock(CU2DBLOCK, CU2DBLOCK);

dim3 dimGrid(n_blocks(NumRows(), CU2DBLOCK),

n_blocks(NumRows(), CU2DBLOCK));

cuda_copy_from_sp(dimGrid, dimBlock, M.Data(), data_, Dim());

CuDevice::Instantiate().AccuProfile("CuMatrix::CopyFromSp",tim.Elapsed());

} else

#endif

{

Mat().CopyFromSp(M.Mat());

}

}

template<typename Real>

CuMatrix<Real>::CuMatrix(const CuMatrix<Real> &other, MatrixTransposeType trans) {

if (trans == kNoTrans)

this->Resize(other.NumRows(), other.NumCols(), kUndefined);

else

this->Resize(other.NumCols(), other.NumRows(), kUndefined);

this->CopyFromMat(other, trans);

}

template<typename Real>

CuMatrix<Real>::CuMatrix(const CuMatrixBase<Real> &other, MatrixTransposeType trans) {

if (trans == kNoTrans)

this->Resize(other.NumRows(), other.NumCols(), kUndefined);

else

this->Resize(other.NumCols(), other.NumRows(), kUndefined);

this->CopyFromMat(other, trans);

}

template<typename Real>

template<typename OtherReal>

CuMatrix<Real>::CuMatrix(const MatrixBase<OtherReal> &other, MatrixTransposeType trans) {

if (trans == kNoTrans)

this->Resize(other.NumRows(), other.NumCols(), kUndefined);

else

this->Resize(other.NumCols(), other.NumRows(), kUndefined);

this->CopyFromMat(other, trans);

}

// Instantiate the template above.

template

CuMatrix<float>::CuMatrix(const MatrixBase<float> &other, MatrixTransposeType trans);

template

CuMatrix<double>::CuMatrix(const MatrixBase<float> &other, MatrixTransposeType trans);

template

CuMatrix<float>::CuMatrix(const MatrixBase<double> &other, MatrixTransposeType trans);

template

CuMatrix<double>::CuMatrix(const MatrixBase<double> &other, MatrixTransposeType trans);

template<typename Real>

template<typename OtherReal>

void CuMatrixBase<Real>::CopyToMat(MatrixBase<OtherReal> *dst,

MatrixTransposeType trans) const {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

if (trans == kTrans || sizeof(OtherReal) != sizeof(Real)) {

CuMatrix<OtherReal> this_trans(*this, trans);

this_trans.CopyToMat(dst, kNoTrans);

} else {

KALDI_ASSERT(dst->NumRows() == NumRows() && dst->NumCols() == NumCols());

if (num_rows_ == 0) return;

Timer tim;

MatrixIndexT src_pitch = stride_*sizeof(Real);

MatrixIndexT dst_pitch = dst->Stride()*sizeof(Real);

MatrixIndexT width = NumCols()*sizeof(Real);

CU_SAFE_CALL(cudaMemcpy2D(dst->Data(), dst_pitch, this->data_, src_pitch,

width, this->num_rows_, cudaMemcpyDeviceToHost));

CuDevice::Instantiate().AccuProfile("CuMatrix::CopyToMatD2H",tim.Elapsed());

}

} else

#endif

{

dst->CopyFromMat(Mat(), trans);

}

}

// instantiate the template above.

template

void CuMatrixBase<float>::CopyToMat(MatrixBase<float> *dst,

MatrixTransposeType trans) const;

template

void CuMatrixBase<double>::CopyToMat(MatrixBase<float> *dst,

MatrixTransposeType trans) const;

template

void CuMatrixBase<float>::CopyToMat(MatrixBase<double> *dst,

MatrixTransposeType trans) const;

template

void CuMatrixBase<double>::CopyToMat(MatrixBase<double> *dst,

MatrixTransposeType trans) const;

template<typename Real>

void CuMatrix<Real>::Read(std::istream &is, bool binary) {

Matrix<Real> temp;

temp.Read(is, binary);

Destroy();

Swap(&temp);

}

template<typename Real>

void CuMatrixBase<Real>::Write(std::ostream &os, bool binary) const {

Matrix<Real> temp(this->num_rows_, this->num_cols_, kUndefined);

this->CopyToMat(&temp);

temp.Write(os, binary);

}

template<typename Real>

void CuMatrixBase<Real>::SetZero() {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

Timer tim;

CU_SAFE_CALL(cudaMemset2D(data_, stride_ * sizeof(Real), 0,

num_cols_ * sizeof(Real), num_rows_ ));

CuDevice::Instantiate().AccuProfile("CuMatrix::SetZero", tim.Elapsed());

} else

#endif

{

Mat().SetZero();

}

}

/*

* Methods wrapping the ANSI-C CUDA kernels

*/

template<typename Real>

void CuMatrixBase<Real>::Set(Real value) {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

if (num_rows_ == 0) return;

Timer tim;

dim3 dimGrid, dimBlock;

GetBlockSizesForSimpleMatrixOperation(NumRows(), NumCols(),

&dimGrid, &dimBlock);

cuda_set_const(dimGrid, dimBlock, data_, value, Dim());

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

Mat().Set(value);

}

}

// set zero the elements above the diagonal.

template<typename Real>

void CuMatrixBase<Real>::SetZeroAboveDiag() {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

if (num_rows_ == 0) return;

Timer tim;

dim3 dimGrid, dimBlock;

GetBlockSizesForSimpleMatrixOperation(NumRows(), NumCols(),

&dimGrid, &dimBlock);

cuda_set_zero_above_diag(dimGrid, dimBlock, data_, Dim());

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

MatrixBase<Real> &mat = Mat();

int32 num_rows = mat.NumRows(), num_cols = mat.NumCols();

for (int32 r = 0; r + 1 < num_rows; r++) {

SubVector<Real> vec(mat, r),

vec_part(vec, r + 1, num_cols - (r + 1));

vec_part.SetZero();

}

}

}

template<typename Real>

void CuMatrixBase<Real>::Add(Real value) {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

if (num_rows_ == 0) return;

Timer tim;

dim3 dimGrid, dimBlock;

GetBlockSizesForSimpleMatrixOperation(NumRows(), NumCols(),

&dimGrid, &dimBlock);

cuda_add(dimGrid, dimBlock, data_, value, Dim());

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

Mat().Add(value);

}

}

template<typename Real>

void CuMatrixBase<Real>::AddToDiag(Real value) {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

if (num_rows_ == 0) return;

Timer tim;

// We'll create a fake matrix with "num_diag" rows, one

// columnn, and a stride of "this_stride". The y-value of

// the grid/blocks corresponds to the row, in this kernel.

MatrixIndexT num_diag = std::min(num_rows_, num_cols_),

this_stride = stride_ + 1;

dim3 dimBlock(1, CU1DBLOCK);

dim3 dimGrid(1, n_blocks(num_diag, CU1DBLOCK));

::MatrixDim d = { num_diag, 1, this_stride };

cuda_add(dimGrid, dimBlock, data_, value, d);

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

Mat().AddToDiag(value);

}

}

template<typename Real>

bool CuMatrixBase<Real>::IsUnit(Real tol) const {

// want to return:

//FrobeniusNorm(*this - I) <= tol * NumRows(), i.e.:

//sqrt (trace((*this - I)(*this-I)) <= tol * NumRows()

// trace((*this - I)(*this - I)) <= tol * NumRows()

// trace(*this * *this) + trace(I) - 2 * trace(*this) <= tol * NumRows()

// trace(*this * *this) + dim - 2*this.Trace() <= tol * NumRows()

KALDI_ASSERT(this->NumRows() == this->NumCols());

return (TraceMatMat(*this, *this, kTrans) + this->NumRows() - 2.0 * this->Trace() <=

tol * this->NumRows());

}

template<typename Real>

void CuMatrixBase<Real>::Scale(Real value) {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

if (num_rows_ == 0) return;

Timer tim;

dim3 dimGrid, dimBlock;

GetBlockSizesForSimpleMatrixOperation(NumRows(), NumCols(),

&dimGrid, &dimBlock);

cuda_scale(dimGrid, dimBlock, data_, value, Dim());

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

Mat().Scale(value);

}

}

template<typename Real>

void CuMatrixBase<Real>::ApplyLog() {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

if (num_rows_ == 0) return;

Timer tim;

dim3 dimGrid, dimBlock;

GetBlockSizesForSimpleMatrixOperation(NumRows(), NumCols(),

&dimGrid, &dimBlock);

cuda_apply_log(dimGrid, dimBlock, data_, Dim());

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

Mat().ApplyLog();

}

}

template<typename Real>

void CuMatrixBase<Real>::MulElements(const CuMatrixBase<Real>& A) {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

Timer tim;

KALDI_ASSERT(num_cols_ == A.NumCols());

KALDI_ASSERT(num_rows_ == A.NumRows());

dim3 dimGrid, dimBlock;

GetBlockSizesForSimpleMatrixOperation(NumRows(), NumCols(),

&dimGrid, &dimBlock);

cuda_mul_elements(dimGrid, dimBlock, data_, A.data_, Dim(), A.Stride());

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

Mat().MulElements(A.Mat());

}

}

template<typename Real>

void CuMatrixBase<Real>::DivElements(const CuMatrixBase<Real>& A) {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

Timer tim;

KALDI_ASSERT(num_cols_ == A.NumCols());

KALDI_ASSERT(num_rows_ == A.NumRows());

dim3 dimGrid, dimBlock;

GetBlockSizesForSimpleMatrixOperation(NumRows(), NumCols(),

&dimGrid, &dimBlock);

cuda_div_elements(dimGrid, dimBlock, data_, A.data_, Dim(), A.Stride());

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

Mat().DivElements(A.Mat());

}

}

template<typename Real>

void CuMatrixBase<Real>::Max(const CuMatrixBase<Real>& A) {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

Timer tim;

KALDI_ASSERT(num_cols_ == A.NumCols());

KALDI_ASSERT(num_rows_ == A.NumRows());

dim3 dimGrid, dimBlock;

GetBlockSizesForSimpleMatrixOperation(NumRows(), NumCols(),

&dimGrid, &dimBlock);

cuda_max(dimGrid, dimBlock, data_, A.data_, Dim(), A.Stride());

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

Mat().Max(A.Mat());

}

}

template<typename Real>

void CuMatrixBase<Real>::MulColsVec(const CuVectorBase<Real> &scale) {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

Timer tim;

KALDI_ASSERT(scale.Dim() == NumCols());

dim3 dimGrid, dimBlock;

GetBlockSizesForSimpleMatrixOperation(NumRows(), NumCols(),

&dimGrid, &dimBlock);

cuda_mul_cols_vec(dimGrid, dimBlock, data_, scale.data_, Dim());

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

Mat().MulColsVec(scale.Vec());

}

}

template<typename Real>

void CuMatrixBase<Real>::MulRowsVec(const CuVectorBase<Real> &scale) {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

Timer tim;

KALDI_ASSERT(scale.Dim() == NumRows());

dim3 dimGrid, dimBlock;

GetBlockSizesForSimpleMatrixOperation(NumRows(), NumCols(),

&dimGrid, &dimBlock);

cuda_mul_rows_vec(dimGrid, dimBlock, data_, scale.data_, Dim());

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

Mat().MulRowsVec(scale.Vec());

}

}

template<typename Real>

void CuMatrixBase<Real>::MulRowsGroupMat(const CuMatrixBase<Real> &src) {

KALDI_ASSERT(src.NumCols() > 0);

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

Timer tim;

int group_size = this->NumCols() / src.NumCols();

dim3 dimGrid, dimBlock;

GetBlockSizesForSimpleMatrixOperation(NumRows(), NumCols(),

&dimGrid, &dimBlock);

cuda_mul_rows_group_mat(dimGrid, dimBlock, this->data_, src.data_,

this->Dim(), src.Stride(), group_size);

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

Mat().MulRowsGroupMat(src.Mat());

}

}

template<typename Real>

void CuMatrixBase<Real>::GroupPnormDeriv(const CuMatrixBase<Real> &src1,

const CuMatrixBase<Real> &src2,

Real power) {

KALDI_ASSERT(src2.NumCols() > 0);

int group_size = this->NumCols() / src2.NumCols();

KALDI_ASSERT(this->NumCols() == src2.NumCols() * group_size);

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

Timer tim;

dim3 dimBlock(CU2DBLOCK, CU2DBLOCK);

dim3 dimGrid(n_blocks(NumCols(), CU2DBLOCK),

n_blocks(NumRows(), CU2DBLOCK));

cuda_calc_pnorm_deriv(dimGrid, dimBlock, this->data_, src1.Data(),

src2.Data(), Dim(), src2.Stride(), group_size, power);

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

Mat().GroupPnormDeriv(src1.Mat(), src2.Mat(), power);

}

}

template<typename Real>

void CuMatrixBase<Real>::DiffGroupPnorm(const CuMatrixBase<Real> &in_value,

const CuMatrixBase<Real> &out_value,

const CuMatrixBase<Real> &out_deriv,

Real power) {

KALDI_ASSERT(out_value.NumCols() > 0);

KALDI_ASSERT(out_value.NumCols() == out_deriv.NumCols());

int group_size = this->NumCols() / out_value.NumCols();

KALDI_ASSERT(this->NumCols() == out_value.NumCols() * group_size);

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

Timer tim;

const int kWarpSize = 32;

dim3 dimBlock(kWarpSize, CU1DBLOCK / kWarpSize);

dim3 dimGrid(n_blocks(NumCols(), dimBlock.x),

n_blocks(NumRows(), dimBlock.y));

if (dimGrid.x * dimGrid.y > 1024) {

dimGrid.y = std::max(1024 / dimGrid.x, unsigned(1));

}

cuda_diff_group_pnorm(dimGrid, dimBlock, this->data_, in_value.Data(),

out_value.Data(), out_deriv.Data(), Dim(),

in_value.Stride(), out_value.Stride(),

out_deriv.Stride(), group_size, power);

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

GroupPnormDeriv(in_value, out_value, power);

MulRowsGroupMat(out_deriv);

}

}

template<typename Real>

void CuMatrixBase<Real>::GroupMaxDeriv(const CuMatrixBase<Real> &src1,

const CuMatrixBase<Real> &src2) {

KALDI_ASSERT(src2.NumCols() > 0);

int group_size = this->NumCols() / src2.NumCols();

KALDI_ASSERT(this->NumCols() == src2.NumCols() * group_size);

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

Timer tim;

dim3 dimBlock(CU2DBLOCK, CU2DBLOCK);

dim3 dimGrid(n_blocks(NumCols(), CU2DBLOCK),

n_blocks(NumRows(), CU2DBLOCK));

cuda_calc_group_max_deriv(dimGrid, dimBlock, this->data_,

src1.Data(), src2.Data(), Dim(),

src2.Stride(), group_size);

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

Mat().GroupMaxDeriv(src1.Mat(), src2.Mat());

}

}

template<typename Real>

void CuMatrixBase<Real>::DivRowsVec(const CuVectorBase<Real> &div) {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

Timer tim;

KALDI_ASSERT(div.Dim() == NumRows());

dim3 dimGrid, dimBlock;

GetBlockSizesForSimpleMatrixOperation(NumRows(), NumCols(),

&dimGrid, &dimBlock);

// For large matrix we do more work per thread by limiting the

// the grid size to reduce the block launching overhead.

if (dimGrid.x * dimGrid.y > 1024) {

dimGrid.x = 1024 / dimGrid.y;

if (dimGrid.x == 0) {

dimGrid.x = 1;

}

}

cuda_div_rows_vec(dimGrid, dimBlock, data_, div.data_, Dim());

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

Vector<Real> temp(div.Vec()); // will copy.

temp.InvertElements();

Mat().MulRowsVec(temp);

}

}

template<typename Real>

void CuMatrixBase<Real>::InvertElements() {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

Timer tim;

dim3 dimGrid, dimBlock;

GetBlockSizesForSimpleMatrixOperation(NumRows(), NumCols(),

&dimGrid, &dimBlock);

cuda_invert_elements(dimGrid, dimBlock, data_, Dim());

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

Mat().InvertElements();

}

}

template<typename Real>

void CuMatrixBase<Real>::AddMat(Real alpha, const CuMatrixBase<Real>& A,

MatrixTransposeType transA) {

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

if (transA == kNoTrans) {

KALDI_ASSERT(A.NumRows() == num_rows_ && A.NumCols() == num_cols_);

} else {

KALDI_ASSERT(A.NumCols() == num_rows_ && A.NumRows() == num_cols_);

}

if (num_rows_ == 0) return;

Timer tim;

// This block dimension seems to work better than the

// one from GetBlockSizesForSimpleMatrixOperation().

dim3 dimBlock(CU2DBLOCK, CU2DBLOCK);

dim3 dimGrid(n_blocks(NumCols(), CU2DBLOCK),

n_blocks(NumRows(), CU2DBLOCK));

cuda_add_mat(dimGrid, dimBlock, alpha, A.data_,

data_, Dim(), A.Stride(),

(transA == kTrans ? 1 : 0));

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

Mat().AddMat(alpha, A.Mat(), transA);

}

}

template<typename Real>

void CuMatrixBase<Real>::AddMatBlocks(Real alpha, const CuMatrixBase<Real> &A,

MatrixTransposeType transA) {

if (num_rows_ == 0 || num_cols_ == 0) return;

int32 num_row_blocks, num_col_blocks;

if (transA == kNoTrans) {

KALDI_ASSERT(A.NumRows() % num_rows_ == 0 && A.NumCols() % num_cols_ == 0);

num_row_blocks = A.Mat().NumRows() / num_rows_;

num_col_blocks = A.Mat().NumCols() / num_cols_;

} else {

KALDI_ASSERT(A.NumRows() % num_cols_ == 0 && A.NumCols() % num_rows_ == 0);

num_row_blocks = A.Mat().NumRows() / num_cols_;

num_col_blocks = A.Mat().NumCols() / num_rows_;

}

#if HAVE_CUDA == 1

if (CuDevice::Instantiate().Enabled()) {

Timer tim;

dim3 dimGrid, dimBlock;

GetBlockSizesForSimpleMatrixOperation(NumRows(), NumCols(),

&dimGrid, &dimBlock);

cuda_add_mat_blocks(dimGrid, dimBlock, alpha, A.data_, num_row_blocks,

num_col_blocks, data_, Dim(), A.Stride(),

(transA == kTrans ? 1 : 0));

CU_SAFE_CALL(cudaGetLastError());

CuDevice::Instantiate().AccuProfile(__func__, tim.Elapsed());

} else

#endif

{

int32 nr, nc;

if (transA == kNoTrans) {

nr = num_rows_;

nc = num_cols_;