You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
3074 lines
111 KiB
3074 lines
111 KiB
// matrix/kaldi-matrix.cc
|
|
|
|
// Copyright 2009-2011 Lukas Burget; Ondrej Glembek; Go Vivace Inc.;
|
|
// Microsoft Corporation; Saarland University;
|
|
// Yanmin Qian; Petr Schwarz; Jan Silovsky;
|
|
// 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.
|
|
|
|
#include "matrix/kaldi-matrix.h"
|
|
|
|
namespace kaldi {
|
|
/*
|
|
template<typename Real>
|
|
void MatrixBase<Real>::Invert(Real *log_det, Real *det_sign,
|
|
bool inverse_needed) {
|
|
KALDI_ASSERT(num_rows_ == num_cols_);
|
|
if (num_rows_ == 0) {
|
|
if (det_sign) *det_sign = 1;
|
|
if (log_det) *log_det = 0.0;
|
|
return;
|
|
}
|
|
#ifndef HAVE_ATLAS
|
|
KaldiBlasInt *pivot = new KaldiBlasInt[num_rows_];
|
|
KaldiBlasInt M = num_rows_;
|
|
KaldiBlasInt N = num_cols_;
|
|
KaldiBlasInt LDA = stride_;
|
|
KaldiBlasInt result = -1;
|
|
KaldiBlasInt l_work = std::max<KaldiBlasInt>(1, N);
|
|
Real *p_work;
|
|
void *temp;
|
|
if ((p_work = static_cast<Real*>(
|
|
KALDI_MEMALIGN(16, sizeof(Real)*l_work, &temp))) == NULL) {
|
|
delete[] pivot;
|
|
throw std::bad_alloc();
|
|
}
|
|
|
|
clapack_Xgetrf2(&M, &N, data_, &LDA, pivot, &result);
|
|
const int pivot_offset = 1;
|
|
#else
|
|
int *pivot = new int[num_rows_];
|
|
int result;
|
|
clapack_Xgetrf(num_rows_, num_cols_, data_, stride_, pivot, &result);
|
|
const int pivot_offset = 0;
|
|
#endif
|
|
KALDI_ASSERT(result >= 0 && "Call to CLAPACK sgetrf_ or ATLAS clapack_sgetrf "
|
|
"called with wrong arguments");
|
|
if (result > 0) {
|
|
if (inverse_needed) {
|
|
KALDI_ERR << "Cannot invert: matrix is singular";
|
|
} else {
|
|
if (log_det) *log_det = -std::numeric_limits<Real>::infinity();
|
|
if (det_sign) *det_sign = 0;
|
|
delete[] pivot;
|
|
#ifndef HAVE_ATLAS
|
|
KALDI_MEMALIGN_FREE(p_work);
|
|
#endif
|
|
return;
|
|
}
|
|
}
|
|
if (det_sign != NULL) {
|
|
int sign = 1;
|
|
for (MatrixIndexT i = 0; i < num_rows_; i++)
|
|
if (pivot[i] != static_cast<int>(i) + pivot_offset) sign *= -1;
|
|
*det_sign = sign;
|
|
}
|
|
if (log_det != NULL || det_sign != NULL) { // Compute log determinant.
|
|
if (log_det != NULL) *log_det = 0.0;
|
|
Real prod = 1.0;
|
|
for (MatrixIndexT i = 0; i < num_rows_; i++) {
|
|
prod *= (*this)(i, i);
|
|
if (i == num_rows_ - 1 || std::fabs(prod) < 1.0e-10 ||
|
|
std::fabs(prod) > 1.0e+10) {
|
|
if (log_det != NULL) *log_det += kaldi::Log(std::fabs(prod));
|
|
if (det_sign != NULL) *det_sign *= (prod > 0 ? 1.0 : -1.0);
|
|
prod = 1.0;
|
|
}
|
|
}
|
|
}
|
|
#ifndef HAVE_ATLAS
|
|
if (inverse_needed) clapack_Xgetri2(&M, data_, &LDA, pivot, p_work, &l_work,
|
|
&result);
|
|
delete[] pivot;
|
|
KALDI_MEMALIGN_FREE(p_work);
|
|
#else
|
|
if (inverse_needed)
|
|
clapack_Xgetri(num_rows_, data_, stride_, pivot, &result);
|
|
delete [] pivot;
|
|
#endif
|
|
KALDI_ASSERT(result == 0 && "Call to CLAPACK sgetri_ or ATLAS clapack_sgetri "
|
|
"called with wrong arguments");
|
|
}
|
|
|
|
template<>
|
|
template<>
|
|
void MatrixBase<float>::AddVecVec(const float alpha,
|
|
const VectorBase<float> &a,
|
|
const VectorBase<float> &rb) {
|
|
KALDI_ASSERT(a.Dim() == num_rows_ && rb.Dim() == num_cols_);
|
|
cblas_Xger(a.Dim(), rb.Dim(), alpha, a.Data(), 1, rb.Data(),
|
|
1, data_, stride_);
|
|
}
|
|
|
|
template<typename Real>
|
|
template<typename OtherReal>
|
|
void MatrixBase<Real>::AddVecVec(const Real alpha,
|
|
const VectorBase<OtherReal> &a,
|
|
const VectorBase<OtherReal> &b) {
|
|
KALDI_ASSERT(a.Dim() == num_rows_ && b.Dim() == num_cols_);
|
|
if (num_rows_ * num_cols_ > 100) { // It's probably worth it to allocate
|
|
// temporary vectors of the right type and use BLAS.
|
|
Vector<Real> temp_a(a), temp_b(b);
|
|
cblas_Xger(num_rows_, num_cols_, alpha, temp_a.Data(), 1,
|
|
temp_b.Data(), 1, data_, stride_);
|
|
} else {
|
|
const OtherReal *a_data = a.Data(), *b_data = b.Data();
|
|
Real *row_data = data_;
|
|
for (MatrixIndexT i = 0; i < num_rows_; i++, row_data += stride_) {
|
|
BaseFloat alpha_ai = alpha * a_data[i];
|
|
for (MatrixIndexT j = 0; j < num_cols_; j++)
|
|
row_data[j] += alpha_ai * b_data[j];
|
|
}
|
|
}
|
|
}
|
|
|
|
// instantiate the template above.
|
|
template
|
|
void MatrixBase<float>::AddVecVec(const float alpha,
|
|
const VectorBase<double> &a,
|
|
const VectorBase<double> &b);
|
|
template
|
|
void MatrixBase<double>::AddVecVec(const double alpha,
|
|
const VectorBase<float> &a,
|
|
const VectorBase<float> &b);
|
|
|
|
template<>
|
|
template<>
|
|
void MatrixBase<double>::AddVecVec(const double alpha,
|
|
const VectorBase<double> &a,
|
|
const VectorBase<double> &rb) {
|
|
KALDI_ASSERT(a.Dim() == num_rows_ && rb.Dim() == num_cols_);
|
|
if (num_rows_ == 0) return;
|
|
cblas_Xger(a.Dim(), rb.Dim(), alpha, a.Data(), 1, rb.Data(),
|
|
1, data_, stride_);
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::AddMatMat(const Real alpha,
|
|
const MatrixBase<Real>& A,
|
|
MatrixTransposeType transA,
|
|
const MatrixBase<Real>& B,
|
|
MatrixTransposeType transB,
|
|
const Real beta) {
|
|
KALDI_ASSERT((transA == kNoTrans && transB == kNoTrans && A.num_cols_ == B.num_rows_ && A.num_rows_ == num_rows_ && B.num_cols_ == num_cols_)
|
|
|| (transA == kTrans && transB == kNoTrans && A.num_rows_ == B.num_rows_ && A.num_cols_ == num_rows_ && B.num_cols_ == num_cols_)
|
|
|| (transA == kNoTrans && transB == kTrans && A.num_cols_ == B.num_cols_ && A.num_rows_ == num_rows_ && B.num_rows_ == num_cols_)
|
|
|| (transA == kTrans && transB == kTrans && A.num_rows_ == B.num_cols_ && A.num_cols_ == num_rows_ && B.num_rows_ == num_cols_));
|
|
KALDI_ASSERT(&A != this && &B != this);
|
|
if (num_rows_ == 0) return;
|
|
cblas_Xgemm(alpha, transA, A.data_, A.num_rows_, A.num_cols_, A.stride_,
|
|
transB, B.data_, B.stride_, beta, data_, num_rows_, num_cols_, stride_);
|
|
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::SetMatMatDivMat(const MatrixBase<Real>& A,
|
|
const MatrixBase<Real>& B,
|
|
const MatrixBase<Real>& C) {
|
|
KALDI_ASSERT(A.NumRows() == B.NumRows() && A.NumCols() == B.NumCols());
|
|
KALDI_ASSERT(A.NumRows() == C.NumRows() && A.NumCols() == C.NumCols());
|
|
for (int32 r = 0; r < A.NumRows(); r++) { // each frame...
|
|
for (int32 c = 0; c < A.NumCols(); c++) {
|
|
BaseFloat i = C(r, c), o = B(r, c), od = A(r, c),
|
|
id;
|
|
if (i != 0.0) {
|
|
id = od * (o / i); /// o / i is either zero or "scale".
|
|
} else {
|
|
id = od; /// Just imagine the scale was 1.0. This is somehow true in
|
|
/// expectation; anyway, this case should basically never happen so it doesn't
|
|
/// really matter.
|
|
}
|
|
(*this)(r, c) = id;
|
|
}
|
|
}
|
|
}
|
|
*/
|
|
|
|
//template<typename Real>
|
|
//void MatrixBase<Real>::CopyLowerToUpper() {
|
|
//KALDI_ASSERT(num_rows_ == num_cols_);
|
|
//Real *data = data_;
|
|
//MatrixIndexT num_rows = num_rows_, stride = stride_;
|
|
//for (int32 i = 0; i < num_rows; i++)
|
|
//for (int32 j = 0; j < i; j++)
|
|
//data[j * stride + i ] = data[i * stride + j];
|
|
//}
|
|
|
|
|
|
//template<typename Real>
|
|
//void MatrixBase<Real>::CopyUpperToLower() {
|
|
//KALDI_ASSERT(num_rows_ == num_cols_);
|
|
//Real *data = data_;
|
|
//MatrixIndexT num_rows = num_rows_, stride = stride_;
|
|
//for (int32 i = 0; i < num_rows; i++)
|
|
//for (int32 j = 0; j < i; j++)
|
|
//data[i * stride + j] = data[j * stride + i];
|
|
//}
|
|
|
|
/*
|
|
template<typename Real>
|
|
void MatrixBase<Real>::SymAddMat2(const Real alpha,
|
|
const MatrixBase<Real> &A,
|
|
MatrixTransposeType transA,
|
|
Real beta) {
|
|
KALDI_ASSERT(num_rows_ == num_cols_ &&
|
|
((transA == kNoTrans && A.num_rows_ == num_rows_) ||
|
|
(transA == kTrans && A.num_cols_ == num_cols_)));
|
|
KALDI_ASSERT(A.data_ != data_);
|
|
if (num_rows_ == 0) return;
|
|
|
|
/// When the matrix dimension(this->num_rows_) is not less than 56
|
|
/// and the transpose type transA == kTrans, the cblas_Xsyrk(...)
|
|
/// function will produce NaN in the output. This is a bug in the
|
|
/// ATLAS library. To overcome this, the AddMatMat function, which calls
|
|
/// cblas_Xgemm(...) rather than cblas_Xsyrk(...), is used in this special
|
|
/// sitation.
|
|
/// Wei Shi: Note this bug is observerd for single precision matrix
|
|
/// on a 64-bit machine
|
|
#ifdef HAVE_ATLAS
|
|
if (transA == kTrans && num_rows_ >= 56) {
|
|
this->AddMatMat(alpha, A, kTrans, A, kNoTrans, beta);
|
|
return;
|
|
}
|
|
#endif // HAVE_ATLAS
|
|
|
|
MatrixIndexT A_other_dim = (transA == kNoTrans ? A.num_cols_ : A.num_rows_);
|
|
|
|
// This function call is hard-coded to update the lower triangle.
|
|
cblas_Xsyrk(transA, num_rows_, A_other_dim, alpha, A.Data(),
|
|
A.Stride(), beta, this->data_, this->stride_);
|
|
}
|
|
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::AddMatSmat(const Real alpha,
|
|
const MatrixBase<Real> &A,
|
|
MatrixTransposeType transA,
|
|
const MatrixBase<Real> &B,
|
|
MatrixTransposeType transB,
|
|
const Real beta) {
|
|
KALDI_ASSERT((transA == kNoTrans && transB == kNoTrans && A.num_cols_ == B.num_rows_ && A.num_rows_ == num_rows_ && B.num_cols_ == num_cols_)
|
|
|| (transA == kTrans && transB == kNoTrans && A.num_rows_ == B.num_rows_ && A.num_cols_ == num_rows_ && B.num_cols_ == num_cols_)
|
|
|| (transA == kNoTrans && transB == kTrans && A.num_cols_ == B.num_cols_ && A.num_rows_ == num_rows_ && B.num_rows_ == num_cols_)
|
|
|| (transA == kTrans && transB == kTrans && A.num_rows_ == B.num_cols_ && A.num_cols_ == num_rows_ && B.num_rows_ == num_cols_));
|
|
KALDI_ASSERT(&A != this && &B != this);
|
|
|
|
// We iterate over the columns of B.
|
|
|
|
MatrixIndexT Astride = A.stride_, Bstride = B.stride_, stride = this->stride_,
|
|
Arows = A.num_rows_, Acols = A.num_cols_;
|
|
Real *data = this->data_, *Adata = A.data_, *Bdata = B.data_;
|
|
MatrixIndexT num_cols = this->num_cols_;
|
|
if (transB == kNoTrans) {
|
|
// Iterate over the columns of *this and of B.
|
|
for (MatrixIndexT c = 0; c < num_cols; c++) {
|
|
// for each column of *this, do
|
|
// [this column] = [alpha * A * this column of B] + [beta * this column]
|
|
Xgemv_sparsevec(transA, Arows, Acols, alpha, Adata, Astride,
|
|
Bdata + c, Bstride, beta, data + c, stride);
|
|
}
|
|
} else {
|
|
// Iterate over the columns of *this and the rows of B.
|
|
for (MatrixIndexT c = 0; c < num_cols; c++) {
|
|
// for each column of *this, do
|
|
// [this column] = [alpha * A * this row of B] + [beta * this column]
|
|
Xgemv_sparsevec(transA, Arows, Acols, alpha, Adata, Astride,
|
|
Bdata + (c * Bstride), 1, beta, data + c, stride);
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::AddSmatMat(const Real alpha,
|
|
const MatrixBase<Real> &A,
|
|
MatrixTransposeType transA,
|
|
const MatrixBase<Real> &B,
|
|
MatrixTransposeType transB,
|
|
const Real beta) {
|
|
KALDI_ASSERT((transA == kNoTrans && transB == kNoTrans && A.num_cols_ == B.num_rows_ && A.num_rows_ == num_rows_ && B.num_cols_ == num_cols_)
|
|
|| (transA == kTrans && transB == kNoTrans && A.num_rows_ == B.num_rows_ && A.num_cols_ == num_rows_ && B.num_cols_ == num_cols_)
|
|
|| (transA == kNoTrans && transB == kTrans && A.num_cols_ == B.num_cols_ && A.num_rows_ == num_rows_ && B.num_rows_ == num_cols_)
|
|
|| (transA == kTrans && transB == kTrans && A.num_rows_ == B.num_cols_ && A.num_cols_ == num_rows_ && B.num_rows_ == num_cols_));
|
|
KALDI_ASSERT(&A != this && &B != this);
|
|
|
|
MatrixIndexT Astride = A.stride_, Bstride = B.stride_, stride = this->stride_,
|
|
Brows = B.num_rows_, Bcols = B.num_cols_;
|
|
MatrixTransposeType invTransB = (transB == kTrans ? kNoTrans : kTrans);
|
|
Real *data = this->data_, *Adata = A.data_, *Bdata = B.data_;
|
|
MatrixIndexT num_rows = this->num_rows_;
|
|
if (transA == kNoTrans) {
|
|
// Iterate over the rows of *this and of A.
|
|
for (MatrixIndexT r = 0; r < num_rows; r++) {
|
|
// for each row of *this, do
|
|
// [this row] = [alpha * (this row of A) * B^T] + [beta * this row]
|
|
Xgemv_sparsevec(invTransB, Brows, Bcols, alpha, Bdata, Bstride,
|
|
Adata + (r * Astride), 1, beta, data + (r * stride), 1);
|
|
}
|
|
} else {
|
|
// Iterate over the rows of *this and the columns of A.
|
|
for (MatrixIndexT r = 0; r < num_rows; r++) {
|
|
// for each row of *this, do
|
|
// [this row] = [alpha * (this column of A) * B^T] + [beta * this row]
|
|
Xgemv_sparsevec(invTransB, Brows, Bcols, alpha, Bdata, Bstride,
|
|
Adata + r, Astride, beta, data + (r * stride), 1);
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::AddSpSp(const Real alpha, const SpMatrix<Real> &A_in,
|
|
const SpMatrix<Real> &B_in, const Real beta) {
|
|
MatrixIndexT sz = num_rows_;
|
|
KALDI_ASSERT(sz == num_cols_ && sz == A_in.NumRows() && sz == B_in.NumRows());
|
|
|
|
Matrix<Real> A(A_in), B(B_in);
|
|
// CblasLower or CblasUpper would work below as symmetric matrix is copied
|
|
// fully (to save work, we used the matrix constructor from SpMatrix).
|
|
// CblasLeft means A is on the left: C <-- alpha A B + beta C
|
|
if (sz == 0) return;
|
|
cblas_Xsymm(alpha, sz, A.data_, A.stride_, B.data_, B.stride_, beta, data_, stride_);
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::AddMat(const Real alpha, const MatrixBase<Real>& A,
|
|
MatrixTransposeType transA) {
|
|
if (&A == this) {
|
|
if (transA == kNoTrans) {
|
|
Scale(alpha + 1.0);
|
|
} else {
|
|
KALDI_ASSERT(num_rows_ == num_cols_ && "AddMat: adding to self (transposed): not symmetric.");
|
|
Real *data = data_;
|
|
if (alpha == 1.0) { // common case-- handle separately.
|
|
for (MatrixIndexT row = 0; row < num_rows_; row++) {
|
|
for (MatrixIndexT col = 0; col < row; col++) {
|
|
Real *lower = data + (row * stride_) + col, *upper = data + (col
|
|
* stride_) + row;
|
|
Real sum = *lower + *upper;
|
|
*lower = *upper = sum;
|
|
}
|
|
*(data + (row * stride_) + row) *= 2.0; // diagonal.
|
|
}
|
|
} else {
|
|
for (MatrixIndexT row = 0; row < num_rows_; row++) {
|
|
for (MatrixIndexT col = 0; col < row; col++) {
|
|
Real *lower = data + (row * stride_) + col, *upper = data + (col
|
|
* stride_) + row;
|
|
Real lower_tmp = *lower;
|
|
*lower += alpha * *upper;
|
|
*upper += alpha * lower_tmp;
|
|
}
|
|
*(data + (row * stride_) + row) *= (1.0 + alpha); // diagonal.
|
|
}
|
|
}
|
|
}
|
|
} else {
|
|
int aStride = (int) A.stride_, stride = stride_;
|
|
Real *adata = A.data_, *data = data_;
|
|
if (transA == kNoTrans) {
|
|
KALDI_ASSERT(A.num_rows_ == num_rows_ && A.num_cols_ == num_cols_);
|
|
if (num_rows_ == 0) return;
|
|
for (MatrixIndexT row = 0; row < num_rows_; row++, adata += aStride,
|
|
data += stride) {
|
|
cblas_Xaxpy(num_cols_, alpha, adata, 1, data, 1);
|
|
}
|
|
} else {
|
|
KALDI_ASSERT(A.num_cols_ == num_rows_ && A.num_rows_ == num_cols_);
|
|
if (num_rows_ == 0) return;
|
|
for (MatrixIndexT row = 0; row < num_rows_; row++, adata++, data += stride)
|
|
cblas_Xaxpy(num_cols_, alpha, adata, aStride, data, 1);
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::AddSmat(Real alpha, const SparseMatrix<Real> &A,
|
|
MatrixTransposeType trans) {
|
|
if (trans == kNoTrans) {
|
|
KALDI_ASSERT(NumRows() == A.NumRows());
|
|
KALDI_ASSERT(NumCols() == A.NumCols());
|
|
MatrixIndexT a_num_rows = A.NumRows();
|
|
for (MatrixIndexT i = 0; i < a_num_rows; ++i) {
|
|
const SparseVector<Real> &row = A.Row(i);
|
|
MatrixIndexT num_elems = row.NumElements();
|
|
for (MatrixIndexT id = 0; id < num_elems; ++id) {
|
|
(*this)(i, row.GetElement(id).first) += alpha
|
|
* row.GetElement(id).second;
|
|
}
|
|
}
|
|
} else {
|
|
KALDI_ASSERT(NumRows() == A.NumCols());
|
|
KALDI_ASSERT(NumCols() == A.NumRows());
|
|
MatrixIndexT a_num_rows = A.NumRows();
|
|
for (MatrixIndexT i = 0; i < a_num_rows; ++i) {
|
|
const SparseVector<Real> &row = A.Row(i);
|
|
MatrixIndexT num_elems = row.NumElements();
|
|
for (MatrixIndexT id = 0; id < num_elems; ++id) {
|
|
(*this)(row.GetElement(id).first, i) += alpha
|
|
* row.GetElement(id).second;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::AddSmatMat(Real alpha, const SparseMatrix<Real> &A,
|
|
MatrixTransposeType transA,
|
|
const MatrixBase<Real> &B, Real beta) {
|
|
if (transA == kNoTrans) {
|
|
KALDI_ASSERT(NumRows() == A.NumRows());
|
|
KALDI_ASSERT(NumCols() == B.NumCols());
|
|
KALDI_ASSERT(A.NumCols() == B.NumRows());
|
|
|
|
this->Scale(beta);
|
|
MatrixIndexT a_num_rows = A.NumRows(),
|
|
this_num_cols = this->NumCols();
|
|
for (MatrixIndexT i = 0; i < a_num_rows; ++i) {
|
|
Real *this_row_i = this->RowData(i);
|
|
const SparseVector<Real> &A_row_i = A.Row(i);
|
|
MatrixIndexT num_elems = A_row_i.NumElements();
|
|
for (MatrixIndexT e = 0; e < num_elems; ++e) {
|
|
const std::pair<MatrixIndexT, Real> &p = A_row_i.GetElement(e);
|
|
MatrixIndexT k = p.first;
|
|
Real alpha_A_ik = alpha * p.second;
|
|
const Real *b_row_k = B.RowData(k);
|
|
cblas_Xaxpy(this_num_cols, alpha_A_ik, b_row_k, 1,
|
|
this_row_i, 1);
|
|
//for (MatrixIndexT j = 0; j < this_num_cols; ++j)
|
|
// this_row_i[j] += alpha_A_ik * b_row_k[j];
|
|
}
|
|
}
|
|
} else {
|
|
KALDI_ASSERT(NumRows() == A.NumCols());
|
|
KALDI_ASSERT(NumCols() == B.NumCols());
|
|
KALDI_ASSERT(A.NumRows() == B.NumRows());
|
|
|
|
this->Scale(beta);
|
|
Matrix<Real> buf(NumRows(), NumCols(), kSetZero);
|
|
MatrixIndexT a_num_rows = A.NumRows(),
|
|
this_num_cols = this->NumCols();
|
|
for (int k = 0; k < a_num_rows; ++k) {
|
|
const Real *b_row_k = B.RowData(k);
|
|
const SparseVector<Real> &A_row_k = A.Row(k);
|
|
MatrixIndexT num_elems = A_row_k.NumElements();
|
|
for (MatrixIndexT e = 0; e < num_elems; ++e) {
|
|
const std::pair<MatrixIndexT, Real> &p = A_row_k.GetElement(e);
|
|
MatrixIndexT i = p.first;
|
|
Real alpha_A_ki = alpha * p.second;
|
|
Real *this_row_i = this->RowData(i);
|
|
cblas_Xaxpy(this_num_cols, alpha_A_ki, b_row_k, 1,
|
|
this_row_i, 1);
|
|
//for (MatrixIndexT j = 0; j < this_num_cols; ++j)
|
|
// this_row_i[j] += alpha_A_ki * b_row_k[j];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::AddMatSmat(Real alpha, const MatrixBase<Real> &A,
|
|
const SparseMatrix<Real> &B,
|
|
MatrixTransposeType transB, Real beta) {
|
|
if (transB == kNoTrans) {
|
|
KALDI_ASSERT(NumRows() == A.NumRows());
|
|
KALDI_ASSERT(NumCols() == B.NumCols());
|
|
KALDI_ASSERT(A.NumCols() == B.NumRows());
|
|
|
|
this->Scale(beta);
|
|
MatrixIndexT b_num_rows = B.NumRows(),
|
|
this_num_rows = this->NumRows();
|
|
// Iterate over the rows of sparse matrix B and columns of A.
|
|
for (MatrixIndexT k = 0; k < b_num_rows; ++k) {
|
|
const SparseVector<Real> &B_row_k = B.Row(k);
|
|
MatrixIndexT num_elems = B_row_k.NumElements();
|
|
const Real *a_col_k = A.Data() + k;
|
|
for (MatrixIndexT e = 0; e < num_elems; ++e) {
|
|
const std::pair<MatrixIndexT, Real> &p = B_row_k.GetElement(e);
|
|
MatrixIndexT j = p.first;
|
|
Real alpha_B_kj = alpha * p.second;
|
|
Real *this_col_j = this->Data() + j;
|
|
// Add to entire 'j'th column of *this at once using cblas_Xaxpy.
|
|
// pass stride to write a colmun as matrices are stored in row major order.
|
|
cblas_Xaxpy(this_num_rows, alpha_B_kj, a_col_k, A.stride_,
|
|
this_col_j, this->stride_);
|
|
//for (MatrixIndexT i = 0; i < this_num_rows; ++i)
|
|
// this_col_j[i*this->stride_] += alpha_B_kj * a_col_k[i*A.stride_];
|
|
}
|
|
}
|
|
} else {
|
|
KALDI_ASSERT(NumRows() == A.NumRows());
|
|
KALDI_ASSERT(NumCols() == B.NumRows());
|
|
KALDI_ASSERT(A.NumCols() == B.NumCols());
|
|
|
|
this->Scale(beta);
|
|
MatrixIndexT b_num_rows = B.NumRows(),
|
|
this_num_rows = this->NumRows();
|
|
// Iterate over the rows of sparse matrix B and columns of *this.
|
|
for (MatrixIndexT j = 0; j < b_num_rows; ++j) {
|
|
const SparseVector<Real> &B_row_j = B.Row(j);
|
|
MatrixIndexT num_elems = B_row_j.NumElements();
|
|
Real *this_col_j = this->Data() + j;
|
|
for (MatrixIndexT e = 0; e < num_elems; ++e) {
|
|
const std::pair<MatrixIndexT, Real> &p = B_row_j.GetElement(e);
|
|
MatrixIndexT k = p.first;
|
|
Real alpha_B_jk = alpha * p.second;
|
|
const Real *a_col_k = A.Data() + k;
|
|
// Add to entire 'j'th column of *this at once using cblas_Xaxpy.
|
|
// pass stride to write a column as matrices are stored in row major order.
|
|
cblas_Xaxpy(this_num_rows, alpha_B_jk, a_col_k, A.stride_,
|
|
this_col_j, this->stride_);
|
|
//for (MatrixIndexT i = 0; i < this_num_rows; ++i)
|
|
// this_col_j[i*this->stride_] += alpha_B_jk * a_col_k[i*A.stride_];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
template<typename OtherReal>
|
|
void MatrixBase<Real>::AddSp(const Real alpha, const SpMatrix<OtherReal> &S) {
|
|
KALDI_ASSERT(S.NumRows() == NumRows() && S.NumRows() == NumCols());
|
|
Real *data = data_; const OtherReal *sdata = S.Data();
|
|
MatrixIndexT num_rows = NumRows(), stride = Stride();
|
|
for (MatrixIndexT i = 0; i < num_rows; i++) {
|
|
for (MatrixIndexT j = 0; j < i; j++, sdata++) {
|
|
data[i*stride + j] += alpha * *sdata;
|
|
data[j*stride + i] += alpha * *sdata;
|
|
}
|
|
data[i*stride + i] += alpha * *sdata++;
|
|
}
|
|
}
|
|
|
|
// instantiate the template above.
|
|
template
|
|
void MatrixBase<float>::AddSp(const float alpha, const SpMatrix<float> &S);
|
|
template
|
|
void MatrixBase<double>::AddSp(const double alpha, const SpMatrix<double> &S);
|
|
template
|
|
void MatrixBase<float>::AddSp(const float alpha, const SpMatrix<double> &S);
|
|
template
|
|
void MatrixBase<double>::AddSp(const double alpha, const SpMatrix<float> &S);
|
|
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::AddDiagVecMat(
|
|
const Real alpha, const VectorBase<Real> &v,
|
|
const MatrixBase<Real> &M,
|
|
MatrixTransposeType transM,
|
|
Real beta) {
|
|
if (beta != 1.0) this->Scale(beta);
|
|
|
|
if (transM == kNoTrans) {
|
|
KALDI_ASSERT(SameDim(*this, M));
|
|
} else {
|
|
KALDI_ASSERT(M.NumRows() == NumCols() && M.NumCols() == NumRows());
|
|
}
|
|
KALDI_ASSERT(v.Dim() == this->NumRows());
|
|
|
|
MatrixIndexT M_row_stride = M.Stride(), M_col_stride = 1, stride = stride_,
|
|
num_rows = num_rows_, num_cols = num_cols_;
|
|
if (transM == kTrans) std::swap(M_row_stride, M_col_stride);
|
|
Real *data = data_;
|
|
const Real *Mdata = M.Data(), *vdata = v.Data();
|
|
if (num_rows_ == 0) return;
|
|
for (MatrixIndexT i = 0; i < num_rows; i++, data += stride, Mdata += M_row_stride, vdata++)
|
|
cblas_Xaxpy(num_cols, alpha * *vdata, Mdata, M_col_stride, data, 1);
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::AddMatDiagVec(
|
|
const Real alpha,
|
|
const MatrixBase<Real> &M, MatrixTransposeType transM,
|
|
VectorBase<Real> &v,
|
|
Real beta) {
|
|
|
|
if (beta != 1.0) this->Scale(beta);
|
|
|
|
if (transM == kNoTrans) {
|
|
KALDI_ASSERT(SameDim(*this, M));
|
|
} else {
|
|
KALDI_ASSERT(M.NumRows() == NumCols() && M.NumCols() == NumRows());
|
|
}
|
|
KALDI_ASSERT(v.Dim() == this->NumCols());
|
|
|
|
MatrixIndexT M_row_stride = M.Stride(),
|
|
M_col_stride = 1,
|
|
stride = stride_,
|
|
num_rows = num_rows_,
|
|
num_cols = num_cols_;
|
|
|
|
if (transM == kTrans)
|
|
std::swap(M_row_stride, M_col_stride);
|
|
|
|
Real *data = data_;
|
|
const Real *Mdata = M.Data(), *vdata = v.Data();
|
|
if (num_rows_ == 0) return;
|
|
for (MatrixIndexT i = 0; i < num_rows; i++){
|
|
for(MatrixIndexT j = 0; j < num_cols; j ++ ){
|
|
data[i*stride + j] += alpha * vdata[j] * Mdata[i*M_row_stride + j*M_col_stride];
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::AddMatMatElements(const Real alpha,
|
|
const MatrixBase<Real>& A,
|
|
const MatrixBase<Real>& B,
|
|
const Real beta) {
|
|
KALDI_ASSERT(A.NumRows() == B.NumRows() && A.NumCols() == B.NumCols());
|
|
KALDI_ASSERT(A.NumRows() == NumRows() && A.NumCols() == NumCols());
|
|
Real *data = data_;
|
|
const Real *dataA = A.Data();
|
|
const Real *dataB = B.Data();
|
|
|
|
for (MatrixIndexT i = 0; i < num_rows_; i++) {
|
|
for (MatrixIndexT j = 0; j < num_cols_; j++) {
|
|
data[j] = beta*data[j] + alpha*dataA[j]*dataB[j];
|
|
}
|
|
data += Stride();
|
|
dataA += A.Stride();
|
|
dataB += B.Stride();
|
|
}
|
|
}
|
|
|
|
#if !defined(HAVE_ATLAS) && !defined(USE_KALDI_SVD)
|
|
// ****************************************************************************
|
|
// ****************************************************************************
|
|
template<typename Real>
|
|
void MatrixBase<Real>::LapackGesvd(VectorBase<Real> *s, MatrixBase<Real> *U_in,
|
|
MatrixBase<Real> *V_in) {
|
|
KALDI_ASSERT(s != NULL && U_in != this && V_in != this);
|
|
|
|
Matrix<Real> tmpU, tmpV;
|
|
if (U_in == NULL) tmpU.Resize(this->num_rows_, 1); // work-space if U_in empty.
|
|
if (V_in == NULL) tmpV.Resize(1, this->num_cols_); // work-space if V_in empty.
|
|
|
|
/// Impementation notes:
|
|
/// Lapack works in column-order, therefore the dimensions of *this are
|
|
/// swapped as well as the U and V matrices.
|
|
|
|
KaldiBlasInt M = num_cols_;
|
|
KaldiBlasInt N = num_rows_;
|
|
KaldiBlasInt LDA = Stride();
|
|
|
|
KALDI_ASSERT(N>=M); // NumRows >= columns.
|
|
|
|
if (U_in) {
|
|
KALDI_ASSERT((int)U_in->num_rows_ == N && (int)U_in->num_cols_ == M);
|
|
}
|
|
if (V_in) {
|
|
KALDI_ASSERT((int)V_in->num_rows_ == M && (int)V_in->num_cols_ == M);
|
|
}
|
|
KALDI_ASSERT((int)s->Dim() == std::min(M, N));
|
|
|
|
MatrixBase<Real> *U = (U_in ? U_in : &tmpU);
|
|
MatrixBase<Real> *V = (V_in ? V_in : &tmpV);
|
|
|
|
KaldiBlasInt V_stride = V->Stride();
|
|
KaldiBlasInt U_stride = U->Stride();
|
|
|
|
// Original LAPACK recipe
|
|
// KaldiBlasInt l_work = std::max(std::max<long int>
|
|
// (1, 3*std::min(M, N)+std::max(M, N)), 5*std::min(M, N))*2;
|
|
KaldiBlasInt l_work = -1;
|
|
Real work_query;
|
|
KaldiBlasInt result;
|
|
|
|
// query for work space
|
|
char *u_job = const_cast<char*>(U_in ? "s" : "N"); // "s" == skinny, "N" == "none."
|
|
char *v_job = const_cast<char*>(V_in ? "s" : "N"); // "s" == skinny, "N" == "none."
|
|
clapack_Xgesvd(v_job, u_job,
|
|
&M, &N, data_, &LDA,
|
|
s->Data(),
|
|
V->Data(), &V_stride,
|
|
U->Data(), &U_stride,
|
|
&work_query, &l_work,
|
|
&result);
|
|
|
|
KALDI_ASSERT(result >= 0 && "Call to CLAPACK dgesvd_ called with wrong arguments");
|
|
|
|
l_work = static_cast<KaldiBlasInt>(work_query);
|
|
Real *p_work;
|
|
void *temp;
|
|
if ((p_work = static_cast<Real*>(
|
|
KALDI_MEMALIGN(16, sizeof(Real)*l_work, &temp))) == NULL)
|
|
throw std::bad_alloc();
|
|
|
|
// perform svd
|
|
clapack_Xgesvd(v_job, u_job,
|
|
&M, &N, data_, &LDA,
|
|
s->Data(),
|
|
V->Data(), &V_stride,
|
|
U->Data(), &U_stride,
|
|
p_work, &l_work,
|
|
&result);
|
|
|
|
KALDI_ASSERT(result >= 0 && "Call to CLAPACK dgesvd_ called with wrong arguments");
|
|
|
|
if (result != 0) {
|
|
KALDI_WARN << "CLAPACK sgesvd_ : some weird convergence not satisfied";
|
|
}
|
|
KALDI_MEMALIGN_FREE(p_work);
|
|
}
|
|
|
|
#endif
|
|
*/
|
|
// Copy constructor. Copies data to newly allocated memory.
|
|
template<typename Real>
|
|
Matrix<Real>::Matrix (const MatrixBase<Real> & M,
|
|
MatrixTransposeType trans/*=kNoTrans*/)
|
|
: MatrixBase<Real>() {
|
|
if (trans == kNoTrans) {
|
|
Resize(M.num_rows_, M.num_cols_);
|
|
this->CopyFromMat(M);
|
|
} else {
|
|
Resize(M.num_cols_, M.num_rows_);
|
|
this->CopyFromMat(M, kTrans);
|
|
}
|
|
}
|
|
|
|
// Copy constructor. Copies data to newly allocated memory.
|
|
template<typename Real>
|
|
Matrix<Real>::Matrix (const Matrix<Real> & M):
|
|
MatrixBase<Real>() {
|
|
Resize(M.num_rows_, M.num_cols_);
|
|
this->CopyFromMat(M);
|
|
}
|
|
|
|
/// Copy constructor from another type.
|
|
template<typename Real>
|
|
template<typename OtherReal>
|
|
Matrix<Real>::Matrix(const MatrixBase<OtherReal> & M,
|
|
MatrixTransposeType trans) : MatrixBase<Real>() {
|
|
if (trans == kNoTrans) {
|
|
Resize(M.NumRows(), M.NumCols());
|
|
this->CopyFromMat(M);
|
|
} else {
|
|
Resize(M.NumCols(), M.NumRows());
|
|
this->CopyFromMat(M, kTrans);
|
|
}
|
|
}
|
|
|
|
// Instantiate this constructor for float->double and double->float.
|
|
template
|
|
Matrix<float>::Matrix(const MatrixBase<double> & M,
|
|
MatrixTransposeType trans);
|
|
template
|
|
Matrix<double>::Matrix(const MatrixBase<float> & M,
|
|
MatrixTransposeType trans);
|
|
|
|
template<typename Real>
|
|
inline void Matrix<Real>::Init(const MatrixIndexT rows,
|
|
const MatrixIndexT cols,
|
|
const MatrixStrideType stride_type) {
|
|
if (rows * cols == 0) {
|
|
KALDI_ASSERT(rows == 0 && cols == 0);
|
|
this->num_rows_ = 0;
|
|
this->num_cols_ = 0;
|
|
this->stride_ = 0;
|
|
this->data_ = NULL;
|
|
return;
|
|
}
|
|
KALDI_ASSERT(rows > 0 && cols > 0);
|
|
MatrixIndexT skip, stride;
|
|
size_t size;
|
|
void *data; // aligned memory block
|
|
void *temp; // memory block to be really freed
|
|
|
|
// compute the size of skip and real cols
|
|
skip = ((16 / sizeof(Real)) - cols % (16 / sizeof(Real)))
|
|
% (16 / sizeof(Real));
|
|
stride = cols + skip;
|
|
size = static_cast<size_t>(rows) * static_cast<size_t>(stride)
|
|
* sizeof(Real);
|
|
|
|
// allocate the memory and set the right dimensions and parameters
|
|
if (NULL != (data = KALDI_MEMALIGN(16, size, &temp))) {
|
|
MatrixBase<Real>::data_ = static_cast<Real *> (data);
|
|
MatrixBase<Real>::num_rows_ = rows;
|
|
MatrixBase<Real>::num_cols_ = cols;
|
|
MatrixBase<Real>::stride_ = (stride_type == kDefaultStride ? stride : cols);
|
|
} else {
|
|
throw std::bad_alloc();
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void Matrix<Real>::Resize(const MatrixIndexT rows,
|
|
const MatrixIndexT cols,
|
|
MatrixResizeType resize_type,
|
|
MatrixStrideType stride_type) {
|
|
// the next block uses recursion to handle what we have to do if
|
|
// resize_type == kCopyData.
|
|
if (resize_type == kCopyData) {
|
|
if (this->data_ == NULL || rows == 0) resize_type = kSetZero; // nothing to copy.
|
|
else if (rows == this->num_rows_ && cols == this->num_cols_ &&
|
|
(stride_type == kDefaultStride || this->stride_ == this->num_cols_)) { return; } // nothing to do.
|
|
else {
|
|
// set tmp to a matrix of the desired size; if new matrix
|
|
// is bigger in some dimension, zero it.
|
|
MatrixResizeType new_resize_type =
|
|
(rows > this->num_rows_ || cols > this->num_cols_) ? kSetZero : kUndefined;
|
|
Matrix<Real> tmp(rows, cols, new_resize_type, stride_type);
|
|
MatrixIndexT rows_min = std::min(rows, this->num_rows_),
|
|
cols_min = std::min(cols, this->num_cols_);
|
|
tmp.Range(0, rows_min, 0, cols_min).
|
|
CopyFromMat(this->Range(0, rows_min, 0, cols_min));
|
|
tmp.Swap(this);
|
|
// and now let tmp go out of scope, deleting what was in *this.
|
|
return;
|
|
}
|
|
}
|
|
// At this point, resize_type == kSetZero or kUndefined.
|
|
|
|
if (MatrixBase<Real>::data_ != NULL) {
|
|
if (rows == MatrixBase<Real>::num_rows_
|
|
&& cols == MatrixBase<Real>::num_cols_) {
|
|
if (resize_type == kSetZero)
|
|
this->SetZero();
|
|
return;
|
|
}
|
|
else
|
|
Destroy();
|
|
}
|
|
Init(rows, cols, stride_type);
|
|
if (resize_type == kSetZero) MatrixBase<Real>::SetZero();
|
|
}
|
|
|
|
template<typename Real>
|
|
template<typename OtherReal>
|
|
void MatrixBase<Real>::CopyFromMat(const MatrixBase<OtherReal> &M,
|
|
MatrixTransposeType Trans) {
|
|
if (sizeof(Real) == sizeof(OtherReal) &&
|
|
static_cast<const void*>(M.Data()) ==
|
|
static_cast<const void*>(this->Data())) {
|
|
// 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 (Trans == kNoTrans) {
|
|
KALDI_ASSERT(num_rows_ == M.NumRows() && num_cols_ == M.NumCols());
|
|
for (MatrixIndexT i = 0; i < num_rows_; i++)
|
|
(*this).Row(i).CopyFromVec(M.Row(i));
|
|
} else {
|
|
KALDI_ASSERT(num_cols_ == M.NumRows() && num_rows_ == M.NumCols());
|
|
int32 this_stride = stride_, other_stride = M.Stride();
|
|
Real *this_data = data_;
|
|
const OtherReal *other_data = M.Data();
|
|
for (MatrixIndexT i = 0; i < num_rows_; i++)
|
|
for (MatrixIndexT j = 0; j < num_cols_; j++)
|
|
this_data[i * this_stride + j] = other_data[j * other_stride + i];
|
|
}
|
|
}
|
|
|
|
// template instantiations.
|
|
template
|
|
void MatrixBase<float>::CopyFromMat(const MatrixBase<double> & M,
|
|
MatrixTransposeType Trans);
|
|
template
|
|
void MatrixBase<double>::CopyFromMat(const MatrixBase<float> & M,
|
|
MatrixTransposeType Trans);
|
|
template
|
|
void MatrixBase<float>::CopyFromMat(const MatrixBase<float> & M,
|
|
MatrixTransposeType Trans);
|
|
template
|
|
void MatrixBase<double>::CopyFromMat(const MatrixBase<double> & M,
|
|
MatrixTransposeType Trans);
|
|
|
|
/*
|
|
// Specialize the template for CopyFromSp for float, float.
|
|
template<>
|
|
template<>
|
|
void MatrixBase<float>::CopyFromSp(const SpMatrix<float> & M) {
|
|
KALDI_ASSERT(num_rows_ == M.NumRows() && num_cols_ == num_rows_);
|
|
MatrixIndexT num_rows = num_rows_, stride = stride_;
|
|
const float *Mdata = M.Data();
|
|
float *row_data = data_, *col_data = data_;
|
|
for (MatrixIndexT i = 0; i < num_rows; i++) {
|
|
cblas_scopy(i+1, Mdata, 1, row_data, 1); // copy to the row.
|
|
cblas_scopy(i, Mdata, 1, col_data, stride); // copy to the column.
|
|
Mdata += i+1;
|
|
row_data += stride;
|
|
col_data += 1;
|
|
}
|
|
}
|
|
|
|
// Specialize the template for CopyFromSp for double, double.
|
|
template<>
|
|
template<>
|
|
void MatrixBase<double>::CopyFromSp(const SpMatrix<double> & M) {
|
|
KALDI_ASSERT(num_rows_ == M.NumRows() && num_cols_ == num_rows_);
|
|
MatrixIndexT num_rows = num_rows_, stride = stride_;
|
|
const double *Mdata = M.Data();
|
|
double *row_data = data_, *col_data = data_;
|
|
for (MatrixIndexT i = 0; i < num_rows; i++) {
|
|
cblas_dcopy(i+1, Mdata, 1, row_data, 1); // copy to the row.
|
|
cblas_dcopy(i, Mdata, 1, col_data, stride); // copy to the column.
|
|
Mdata += i+1;
|
|
row_data += stride;
|
|
col_data += 1;
|
|
}
|
|
}
|
|
|
|
|
|
template<typename Real>
|
|
template<typename OtherReal>
|
|
void MatrixBase<Real>::CopyFromSp(const SpMatrix<OtherReal> & M) {
|
|
KALDI_ASSERT(num_rows_ == M.NumRows() && num_cols_ == num_rows_);
|
|
// MORE EFFICIENT IF LOWER TRIANGULAR! Reverse code otherwise.
|
|
for (MatrixIndexT i = 0; i < num_rows_; i++) {
|
|
for (MatrixIndexT j = 0; j < i; j++) {
|
|
(*this)(j, i) = (*this)(i, j) = M(i, j);
|
|
}
|
|
(*this)(i, i) = M(i, i);
|
|
}
|
|
}
|
|
|
|
// Instantiate this function
|
|
template
|
|
void MatrixBase<float>::CopyFromSp(const SpMatrix<double> & M);
|
|
template
|
|
void MatrixBase<double>::CopyFromSp(const SpMatrix<float> & M);
|
|
|
|
|
|
template<typename Real>
|
|
template<typename OtherReal>
|
|
void MatrixBase<Real>::CopyFromTp(const TpMatrix<OtherReal> & M,
|
|
MatrixTransposeType Trans) {
|
|
if (Trans == kNoTrans) {
|
|
KALDI_ASSERT(num_rows_ == M.NumRows() && num_cols_ == num_rows_);
|
|
SetZero();
|
|
Real *out_i = data_;
|
|
const OtherReal *in_i = M.Data();
|
|
for (MatrixIndexT i = 0; i < num_rows_; i++, out_i += stride_, in_i += i) {
|
|
for (MatrixIndexT j = 0; j <= i; j++)
|
|
out_i[j] = in_i[j];
|
|
}
|
|
} else {
|
|
SetZero();
|
|
KALDI_ASSERT(num_rows_ == M.NumRows() && num_cols_ == num_rows_);
|
|
MatrixIndexT stride = stride_;
|
|
Real *out_i = data_;
|
|
const OtherReal *in_i = M.Data();
|
|
for (MatrixIndexT i = 0; i < num_rows_; i++, out_i ++, in_i += i) {
|
|
for (MatrixIndexT j = 0; j <= i; j++)
|
|
out_i[j*stride] = in_i[j];
|
|
}
|
|
}
|
|
}
|
|
|
|
template
|
|
void MatrixBase<float>::CopyFromTp(const TpMatrix<float> & M,
|
|
MatrixTransposeType trans);
|
|
template
|
|
void MatrixBase<float>::CopyFromTp(const TpMatrix<double> & M,
|
|
MatrixTransposeType trans);
|
|
template
|
|
void MatrixBase<double>::CopyFromTp(const TpMatrix<float> & M,
|
|
MatrixTransposeType trans);
|
|
template
|
|
void MatrixBase<double>::CopyFromTp(const TpMatrix<double> & M,
|
|
MatrixTransposeType trans);
|
|
|
|
*/
|
|
template<typename Real>
|
|
void MatrixBase<Real>::CopyRowsFromVec(const VectorBase<Real> &rv) {
|
|
if (rv.Dim() == num_rows_*num_cols_) {
|
|
if (stride_ == num_cols_) {
|
|
// one big copy operation.
|
|
const Real *rv_data = rv.Data();
|
|
std::memcpy(data_, rv_data, sizeof(Real)*num_rows_*num_cols_);
|
|
} else {
|
|
const Real *rv_data = rv.Data();
|
|
for (MatrixIndexT r = 0; r < num_rows_; r++) {
|
|
Real *row_data = RowData(r);
|
|
for (MatrixIndexT c = 0; c < num_cols_; c++) {
|
|
row_data[c] = rv_data[c];
|
|
}
|
|
rv_data += num_cols_;
|
|
}
|
|
}
|
|
} else if (rv.Dim() == num_cols_) {
|
|
const Real *rv_data = rv.Data();
|
|
for (MatrixIndexT r = 0; r < num_rows_; r++)
|
|
std::memcpy(RowData(r), rv_data, sizeof(Real)*num_cols_);
|
|
} else {
|
|
KALDI_ERR << "Wrong sized arguments";
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
template<typename OtherReal>
|
|
void MatrixBase<Real>::CopyRowsFromVec(const VectorBase<OtherReal> &rv) {
|
|
if (rv.Dim() == num_rows_*num_cols_) {
|
|
const OtherReal *rv_data = rv.Data();
|
|
for (MatrixIndexT r = 0; r < num_rows_; r++) {
|
|
Real *row_data = RowData(r);
|
|
for (MatrixIndexT c = 0; c < num_cols_; c++) {
|
|
row_data[c] = static_cast<Real>(rv_data[c]);
|
|
}
|
|
rv_data += num_cols_;
|
|
}
|
|
} else if (rv.Dim() == num_cols_) {
|
|
const OtherReal *rv_data = rv.Data();
|
|
Real *first_row_data = RowData(0);
|
|
for (MatrixIndexT c = 0; c < num_cols_; c++)
|
|
first_row_data[c] = rv_data[c];
|
|
for (MatrixIndexT r = 1; r < num_rows_; r++)
|
|
std::memcpy(RowData(r), first_row_data, sizeof(Real)*num_cols_);
|
|
} else {
|
|
KALDI_ERR << "Wrong sized arguments.";
|
|
}
|
|
}
|
|
|
|
|
|
template
|
|
void MatrixBase<float>::CopyRowsFromVec(const VectorBase<double> &rv);
|
|
template
|
|
void MatrixBase<double>::CopyRowsFromVec(const VectorBase<float> &rv);
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::CopyColsFromVec(const VectorBase<Real> &rv) {
|
|
if (rv.Dim() == num_rows_*num_cols_) {
|
|
const Real *v_inc_data = rv.Data();
|
|
Real *m_inc_data = data_;
|
|
|
|
for (MatrixIndexT c = 0; c < num_cols_; c++) {
|
|
for (MatrixIndexT r = 0; r < num_rows_; r++) {
|
|
m_inc_data[r * stride_] = v_inc_data[r];
|
|
}
|
|
v_inc_data += num_rows_;
|
|
m_inc_data ++;
|
|
}
|
|
} else if (rv.Dim() == num_rows_) {
|
|
const Real *v_inc_data = rv.Data();
|
|
Real *m_inc_data = data_;
|
|
for (MatrixIndexT r = 0; r < num_rows_; r++) {
|
|
Real value = *(v_inc_data++);
|
|
for (MatrixIndexT c = 0; c < num_cols_; c++)
|
|
m_inc_data[c] = value;
|
|
m_inc_data += stride_;
|
|
}
|
|
} else {
|
|
KALDI_ERR << "Wrong size of arguments.";
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::CopyRowFromVec(const VectorBase<Real> &rv, const MatrixIndexT row) {
|
|
KALDI_ASSERT(rv.Dim() == num_cols_ &&
|
|
static_cast<UnsignedMatrixIndexT>(row) <
|
|
static_cast<UnsignedMatrixIndexT>(num_rows_));
|
|
|
|
const Real *rv_data = rv.Data();
|
|
Real *row_data = RowData(row);
|
|
|
|
std::memcpy(row_data, rv_data, num_cols_ * sizeof(Real));
|
|
}
|
|
/*
|
|
template<typename Real>
|
|
void MatrixBase<Real>::CopyDiagFromVec(const VectorBase<Real> &rv) {
|
|
KALDI_ASSERT(rv.Dim() == std::min(num_cols_, num_rows_));
|
|
const Real *rv_data = rv.Data(), *rv_end = rv_data + rv.Dim();
|
|
Real *my_data = this->Data();
|
|
for (; rv_data != rv_end; rv_data++, my_data += (this->stride_+1))
|
|
*my_data = *rv_data;
|
|
}*/
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::CopyColFromVec(const VectorBase<Real> &rv,
|
|
const MatrixIndexT col) {
|
|
KALDI_ASSERT(rv.Dim() == num_rows_ &&
|
|
static_cast<UnsignedMatrixIndexT>(col) <
|
|
static_cast<UnsignedMatrixIndexT>(num_cols_));
|
|
|
|
const Real *rv_data = rv.Data();
|
|
Real *col_data = data_ + col;
|
|
|
|
for (MatrixIndexT r = 0; r < num_rows_; r++)
|
|
col_data[r * stride_] = rv_data[r];
|
|
}
|
|
|
|
|
|
|
|
template<typename Real>
|
|
void Matrix<Real>::RemoveRow(MatrixIndexT i) {
|
|
KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(i) <
|
|
static_cast<UnsignedMatrixIndexT>(MatrixBase<Real>::num_rows_)
|
|
&& "Access out of matrix");
|
|
for (MatrixIndexT j = i + 1; j < MatrixBase<Real>::num_rows_; j++)
|
|
MatrixBase<Real>::Row(j-1).CopyFromVec( MatrixBase<Real>::Row(j));
|
|
MatrixBase<Real>::num_rows_--;
|
|
}
|
|
|
|
template<typename Real>
|
|
void Matrix<Real>::Destroy() {
|
|
// we need to free the data block if it was defined
|
|
if (NULL != MatrixBase<Real>::data_)
|
|
KALDI_MEMALIGN_FREE( MatrixBase<Real>::data_);
|
|
MatrixBase<Real>::data_ = NULL;
|
|
MatrixBase<Real>::num_rows_ = MatrixBase<Real>::num_cols_
|
|
= MatrixBase<Real>::stride_ = 0;
|
|
}
|
|
|
|
|
|
/*
|
|
template<typename Real>
|
|
void MatrixBase<Real>::MulElements(const MatrixBase<Real> &a) {
|
|
KALDI_ASSERT(a.NumRows() == num_rows_ && a.NumCols() == num_cols_);
|
|
|
|
if (num_cols_ == stride_ && num_cols_ == a.stride_) {
|
|
mul_elements(num_rows_ * num_cols_, a.data_, data_);
|
|
} else {
|
|
MatrixIndexT a_stride = a.stride_, stride = stride_;
|
|
Real *data = data_, *a_data = a.data_;
|
|
for (MatrixIndexT i = 0; i < num_rows_; i++) {
|
|
mul_elements(num_cols_, a_data, data);
|
|
a_data += a_stride;
|
|
data += stride;
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::DivElements(const MatrixBase<Real> &a) {
|
|
KALDI_ASSERT(a.NumRows() == num_rows_ && a.NumCols() == num_cols_);
|
|
MatrixIndexT i;
|
|
MatrixIndexT j;
|
|
|
|
for (i = 0; i < num_rows_; i++) {
|
|
for (j = 0; j < num_cols_; j++) {
|
|
(*this)(i, j) /= a(i, j);
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
Real MatrixBase<Real>::Sum() const {
|
|
double sum = 0.0;
|
|
|
|
for (MatrixIndexT i = 0; i < num_rows_; i++) {
|
|
for (MatrixIndexT j = 0; j < num_cols_; j++) {
|
|
sum += (*this)(i, j);
|
|
}
|
|
}
|
|
|
|
return (Real)sum;
|
|
}
|
|
|
|
template<typename Real> void MatrixBase<Real>::Max(const MatrixBase<Real> &A) {
|
|
KALDI_ASSERT(A.NumRows() == NumRows() && A.NumCols() == NumCols());
|
|
for (MatrixIndexT row = 0; row < num_rows_; row++) {
|
|
Real *row_data = RowData(row);
|
|
const Real *other_row_data = A.RowData(row);
|
|
MatrixIndexT num_cols = num_cols_;
|
|
for (MatrixIndexT col = 0; col < num_cols; col++) {
|
|
row_data[col] = std::max(row_data[col],
|
|
other_row_data[col]);
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real> void MatrixBase<Real>::Min(const MatrixBase<Real> &A) {
|
|
KALDI_ASSERT(A.NumRows() == NumRows() && A.NumCols() == NumCols());
|
|
for (MatrixIndexT row = 0; row < num_rows_; row++) {
|
|
Real *row_data = RowData(row);
|
|
const Real *other_row_data = A.RowData(row);
|
|
MatrixIndexT num_cols = num_cols_;
|
|
for (MatrixIndexT col = 0; col < num_cols; col++) {
|
|
row_data[col] = std::min(row_data[col],
|
|
other_row_data[col]);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
template<typename Real> void MatrixBase<Real>::Scale(Real alpha) {
|
|
if (alpha == 1.0) return;
|
|
if (num_rows_ == 0) return;
|
|
if (num_cols_ == stride_) {
|
|
cblas_Xscal(static_cast<size_t>(num_rows_) * static_cast<size_t>(num_cols_),
|
|
alpha, data_,1);
|
|
} else {
|
|
Real *data = data_;
|
|
for (MatrixIndexT i = 0; i < num_rows_; ++i, data += stride_) {
|
|
cblas_Xscal(num_cols_, alpha, data,1);
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real> // scales each row by scale[i].
|
|
void MatrixBase<Real>::MulRowsVec(const VectorBase<Real> &scale) {
|
|
KALDI_ASSERT(scale.Dim() == num_rows_);
|
|
MatrixIndexT M = num_rows_, N = num_cols_;
|
|
|
|
for (MatrixIndexT i = 0; i < M; i++) {
|
|
Real this_scale = scale(i);
|
|
for (MatrixIndexT j = 0; j < N; j++) {
|
|
(*this)(i, j) *= this_scale;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::MulRowsGroupMat(const MatrixBase<Real> &src) {
|
|
KALDI_ASSERT(src.NumRows() == this->NumRows() &&
|
|
this->NumCols() % src.NumCols() == 0);
|
|
int32 group_size = this->NumCols() / src.NumCols(),
|
|
num_groups = this->NumCols() / group_size,
|
|
num_rows = this->NumRows();
|
|
|
|
for (MatrixIndexT i = 0; i < num_rows; i++) {
|
|
Real *data = this->RowData(i);
|
|
for (MatrixIndexT j = 0; j < num_groups; j++, data += group_size) {
|
|
Real scale = src(i, j);
|
|
cblas_Xscal(group_size, scale, data, 1);
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::GroupPnormDeriv(const MatrixBase<Real> &input,
|
|
const MatrixBase<Real> &output,
|
|
Real power) {
|
|
KALDI_ASSERT(input.NumCols() == this->NumCols() && input.NumRows() == this->NumRows());
|
|
KALDI_ASSERT(this->NumCols() % output.NumCols() == 0 &&
|
|
this->NumRows() == output.NumRows());
|
|
|
|
int group_size = this->NumCols() / output.NumCols(),
|
|
num_rows = this->NumRows(), num_cols = this->NumCols();
|
|
|
|
if (power == 1.0) {
|
|
for (MatrixIndexT i = 0; i < num_rows; i++) {
|
|
for (MatrixIndexT j = 0; j < num_cols; j++) {
|
|
Real input_val = input(i, j);
|
|
(*this)(i, j) = (input_val == 0 ? 0 : (input_val > 0 ? 1 : -1));
|
|
}
|
|
}
|
|
} else if (power == std::numeric_limits<Real>::infinity()) {
|
|
for (MatrixIndexT i = 0; i < num_rows; i++) {
|
|
for (MatrixIndexT j = 0; j < num_cols; j++) {
|
|
Real output_val = output(i, j / group_size), input_val = input(i, j);
|
|
if (output_val == 0)
|
|
(*this)(i, j) = 0;
|
|
else
|
|
(*this)(i, j) = (std::abs(input_val) == output_val ? 1.0 : 0.0)
|
|
* (input_val >= 0 ? 1 : -1);
|
|
}
|
|
}
|
|
} else {
|
|
for (MatrixIndexT i = 0; i < num_rows; i++) {
|
|
for (MatrixIndexT j = 0; j < num_cols; j++) {
|
|
Real output_val = output(i, j / group_size),
|
|
input_val = input(i, j);
|
|
if (output_val == 0)
|
|
(*this)(i, j) = 0;
|
|
else
|
|
(*this)(i, j) = pow(std::abs(input_val), power - 1) *
|
|
pow(output_val, 1 - power) * (input_val >= 0 ? 1 : -1) ;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::GroupMaxDeriv(const MatrixBase<Real> &input,
|
|
const MatrixBase<Real> &output) {
|
|
KALDI_ASSERT(input.NumCols() == this->NumCols() &&
|
|
input.NumRows() == this->NumRows());
|
|
KALDI_ASSERT(this->NumCols() % output.NumCols() == 0 &&
|
|
this->NumRows() == output.NumRows());
|
|
|
|
int group_size = this->NumCols() / output.NumCols(),
|
|
num_rows = this->NumRows(), num_cols = this->NumCols();
|
|
|
|
for (MatrixIndexT i = 0; i < num_rows; i++) {
|
|
for (MatrixIndexT j = 0; j < num_cols; j++) {
|
|
Real input_val = input(i, j);
|
|
Real output_val = output(i, j / group_size);
|
|
(*this)(i, j) = (input_val == output_val ? 1 : 0);
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real> // scales each column by scale[i].
|
|
void MatrixBase<Real>::MulColsVec(const VectorBase<Real> &scale) {
|
|
KALDI_ASSERT(scale.Dim() == num_cols_);
|
|
for (MatrixIndexT i = 0; i < num_rows_; i++) {
|
|
for (MatrixIndexT j = 0; j < num_cols_; j++) {
|
|
Real this_scale = scale(j);
|
|
(*this)(i, j) *= this_scale;
|
|
}
|
|
}
|
|
}
|
|
*/
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::SetZero() {
|
|
if (num_cols_ == stride_)
|
|
memset(data_, 0, sizeof(Real)*num_rows_*num_cols_);
|
|
else
|
|
for (MatrixIndexT row = 0; row < num_rows_; row++)
|
|
memset(data_ + row*stride_, 0, sizeof(Real)*num_cols_);
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::Set(Real value) {
|
|
for (MatrixIndexT row = 0; row < num_rows_; row++) {
|
|
for (MatrixIndexT col = 0; col < num_cols_; col++) {
|
|
(*this)(row, col) = value;
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
template<typename Real>
|
|
void MatrixBase<Real>::SetUnit() {
|
|
SetZero();
|
|
for (MatrixIndexT row = 0; row < std::min(num_rows_, num_cols_); row++)
|
|
(*this)(row, row) = 1.0;
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::SetRandn() {
|
|
kaldi::RandomState rstate;
|
|
for (MatrixIndexT row = 0; row < num_rows_; row++) {
|
|
Real *row_data = this->RowData(row);
|
|
MatrixIndexT nc = (num_cols_ % 2 == 1) ? num_cols_ - 1 : num_cols_;
|
|
for (MatrixIndexT col = 0; col < nc; col += 2) {
|
|
kaldi::RandGauss2(row_data + col, row_data + col + 1, &rstate);
|
|
}
|
|
if (nc != num_cols_) row_data[nc] = static_cast<Real>(kaldi::RandGauss(&rstate));
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::SetRandUniform() {
|
|
kaldi::RandomState rstate;
|
|
for (MatrixIndexT row = 0; row < num_rows_; row++) {
|
|
Real *row_data = this->RowData(row);
|
|
for (MatrixIndexT col = 0; col < num_cols_; col++, row_data++) {
|
|
*row_data = static_cast<Real>(kaldi::RandUniform(&rstate));
|
|
}
|
|
}
|
|
}
|
|
*/
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::Write(std::ostream &os, bool binary) const {
|
|
if (!os.good()) {
|
|
KALDI_ERR << "Failed to write matrix to stream: stream not good";
|
|
}
|
|
if (binary) { // Use separate binary and text formats,
|
|
// since in binary mode we need to know if it's float or double.
|
|
std::string my_token = (sizeof(Real) == 4 ? "FM" : "DM");
|
|
|
|
WriteToken(os, binary, my_token);
|
|
{
|
|
int32 rows = this->num_rows_; // make the size 32-bit on disk.
|
|
int32 cols = this->num_cols_;
|
|
KALDI_ASSERT(this->num_rows_ == (MatrixIndexT) rows);
|
|
KALDI_ASSERT(this->num_cols_ == (MatrixIndexT) cols);
|
|
WriteBasicType(os, binary, rows);
|
|
WriteBasicType(os, binary, cols);
|
|
}
|
|
if (Stride() == NumCols())
|
|
os.write(reinterpret_cast<const char*> (Data()), sizeof(Real)
|
|
* static_cast<size_t>(num_rows_) * static_cast<size_t>(num_cols_));
|
|
else
|
|
for (MatrixIndexT i = 0; i < num_rows_; i++)
|
|
os.write(reinterpret_cast<const char*> (RowData(i)), sizeof(Real)
|
|
* num_cols_);
|
|
if (!os.good()) {
|
|
KALDI_ERR << "Failed to write matrix to stream";
|
|
}
|
|
} else { // text mode.
|
|
if (num_cols_ == 0) {
|
|
os << " [ ]\n";
|
|
} else {
|
|
os << " [";
|
|
for (MatrixIndexT i = 0; i < num_rows_; i++) {
|
|
os << "\n ";
|
|
for (MatrixIndexT j = 0; j < num_cols_; j++)
|
|
os << (*this)(i, j) << " ";
|
|
}
|
|
os << "]\n";
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::Read(std::istream & is, bool binary) {
|
|
// In order to avoid rewriting this, we just declare a Matrix and
|
|
// use it to read the data, then copy.
|
|
Matrix<Real> tmp;
|
|
tmp.Read(is, binary);
|
|
if (tmp.NumRows() != NumRows() || tmp.NumCols() != NumCols()) {
|
|
KALDI_ERR << "MatrixBase<Real>::Read, size mismatch "
|
|
<< NumRows() << " x " << NumCols() << " versus "
|
|
<< tmp.NumRows() << " x " << tmp.NumCols();
|
|
}
|
|
CopyFromMat(tmp);
|
|
}
|
|
|
|
|
|
template<typename Real>
|
|
void Matrix<Real>::Read(std::istream & is, bool binary) {
|
|
// now assume add == false.
|
|
MatrixIndexT pos_at_start = is.tellg();
|
|
std::ostringstream specific_error;
|
|
|
|
if (binary) { // Read in binary mode.
|
|
int peekval = Peek(is, binary);
|
|
if (peekval == 'C') {
|
|
// This code enables us to read CompressedMatrix as a regular matrix.
|
|
//CompressedMatrix compressed_mat;
|
|
//compressed_mat.Read(is, binary); // at this point, add == false.
|
|
//this->Resize(compressed_mat.NumRows(), compressed_mat.NumCols());
|
|
//compressed_mat.CopyToMat(this);
|
|
return;
|
|
}
|
|
const char *my_token = (sizeof(Real) == 4 ? "FM" : "DM");
|
|
char other_token_start = (sizeof(Real) == 4 ? 'D' : 'F');
|
|
if (peekval == other_token_start) { // need to instantiate the other type to read it.
|
|
typedef typename OtherReal<Real>::Real OtherType; // if Real == float, OtherType == double, and vice versa.
|
|
Matrix<OtherType> other(this->num_rows_, this->num_cols_);
|
|
other.Read(is, binary); // add is false at this point anyway.
|
|
this->Resize(other.NumRows(), other.NumCols());
|
|
this->CopyFromMat(other);
|
|
return;
|
|
}
|
|
std::string token;
|
|
ReadToken(is, binary, &token);
|
|
if (token != my_token) {
|
|
if (token.length() > 20) token = token.substr(0, 17) + "...";
|
|
specific_error << ": Expected token " << my_token << ", got " << token;
|
|
goto bad;
|
|
}
|
|
int32 rows, cols;
|
|
ReadBasicType(is, binary, &rows); // throws on error.
|
|
ReadBasicType(is, binary, &cols); // throws on error.
|
|
if ((MatrixIndexT)rows != this->num_rows_ || (MatrixIndexT)cols != this->num_cols_) {
|
|
this->Resize(rows, cols);
|
|
}
|
|
if (this->Stride() == this->NumCols() && rows*cols!=0) {
|
|
is.read(reinterpret_cast<char*>(this->Data()),
|
|
sizeof(Real)*rows*cols);
|
|
if (is.fail()) goto bad;
|
|
} else {
|
|
for (MatrixIndexT i = 0; i < (MatrixIndexT)rows; i++) {
|
|
is.read(reinterpret_cast<char*>(this->RowData(i)), sizeof(Real)*cols);
|
|
if (is.fail()) goto bad;
|
|
}
|
|
}
|
|
if (is.eof()) return;
|
|
if (is.fail()) goto bad;
|
|
return;
|
|
} else { // Text mode.
|
|
std::string str;
|
|
is >> str; // get a token
|
|
if (is.fail()) { specific_error << ": Expected \"[\", got EOF"; goto bad; }
|
|
// if ((str.compare("DM") == 0) || (str.compare("FM") == 0)) { // Back compatibility.
|
|
// is >> str; // get #rows
|
|
// is >> str; // get #cols
|
|
// is >> str; // get "["
|
|
// }
|
|
if (str == "[]") { Resize(0, 0); return; } // Be tolerant of variants.
|
|
else if (str != "[") {
|
|
if (str.length() > 20) str = str.substr(0, 17) + "...";
|
|
specific_error << ": Expected \"[\", got \"" << str << '"';
|
|
goto bad;
|
|
}
|
|
// At this point, we have read "[".
|
|
std::vector<std::vector<Real>* > data;
|
|
std::vector<Real> *cur_row = new std::vector<Real>;
|
|
while (1) {
|
|
int i = is.peek();
|
|
if (i == -1) { specific_error << "Got EOF while reading matrix data"; goto cleanup; }
|
|
else if (static_cast<char>(i) == ']') { // Finished reading matrix.
|
|
is.get(); // eat the "]".
|
|
i = is.peek();
|
|
if (static_cast<char>(i) == '\r') {
|
|
is.get();
|
|
is.get(); // get \r\n (must eat what we wrote)
|
|
} else if (static_cast<char>(i) == '\n') { is.get(); } // get \n (must eat what we wrote)
|
|
if (is.fail()) {
|
|
KALDI_WARN << "After end of matrix data, read error.";
|
|
// we got the data we needed, so just warn for this error.
|
|
}
|
|
// Now process the data.
|
|
if (!cur_row->empty()) data.push_back(cur_row);
|
|
else delete(cur_row);
|
|
cur_row = NULL;
|
|
if (data.empty()) { this->Resize(0, 0); return; }
|
|
else {
|
|
int32 num_rows = data.size(), num_cols = data[0]->size();
|
|
this->Resize(num_rows, num_cols);
|
|
for (int32 i = 0; i < num_rows; i++) {
|
|
if (static_cast<int32>(data[i]->size()) != num_cols) {
|
|
specific_error << "Matrix has inconsistent #cols: " << num_cols
|
|
<< " vs." << data[i]->size() << " (processing row"
|
|
<< i << ")";
|
|
goto cleanup;
|
|
}
|
|
for (int32 j = 0; j < num_cols; j++)
|
|
(*this)(i, j) = (*(data[i]))[j];
|
|
delete data[i];
|
|
data[i] = NULL;
|
|
}
|
|
}
|
|
return;
|
|
} else if (static_cast<char>(i) == '\n' || static_cast<char>(i) == ';') {
|
|
// End of matrix row.
|
|
is.get();
|
|
if (cur_row->size() != 0) {
|
|
data.push_back(cur_row);
|
|
cur_row = new std::vector<Real>;
|
|
cur_row->reserve(data.back()->size());
|
|
}
|
|
} else if ( (i >= '0' && i <= '9') || i == '-' ) { // A number...
|
|
Real r;
|
|
is >> r;
|
|
if (is.fail()) {
|
|
specific_error << "Stream failure/EOF while reading matrix data.";
|
|
goto cleanup;
|
|
}
|
|
cur_row->push_back(r);
|
|
} else if (isspace(i)) {
|
|
is.get(); // eat the space and do nothing.
|
|
} else { // NaN or inf or error.
|
|
std::string str;
|
|
is >> str;
|
|
if (!KALDI_STRCASECMP(str.c_str(), "inf") ||
|
|
!KALDI_STRCASECMP(str.c_str(), "infinity")) {
|
|
cur_row->push_back(std::numeric_limits<Real>::infinity());
|
|
KALDI_WARN << "Reading infinite value into matrix.";
|
|
} else if (!KALDI_STRCASECMP(str.c_str(), "nan")) {
|
|
cur_row->push_back(std::numeric_limits<Real>::quiet_NaN());
|
|
KALDI_WARN << "Reading NaN value into matrix.";
|
|
} else {
|
|
if (str.length() > 20) str = str.substr(0, 17) + "...";
|
|
specific_error << "Expecting numeric matrix data, got " << str;
|
|
goto cleanup;
|
|
}
|
|
}
|
|
}
|
|
// Note, we never leave the while () loop before this
|
|
// line (we return from it.)
|
|
cleanup: // We only reach here in case of error in the while loop above.
|
|
if(cur_row != NULL)
|
|
delete cur_row;
|
|
for (size_t i = 0; i < data.size(); i++)
|
|
if(data[i] != NULL)
|
|
delete data[i];
|
|
// and then go on to "bad" below, where we print error.
|
|
}
|
|
bad:
|
|
KALDI_ERR << "Failed to read matrix from stream. " << specific_error.str()
|
|
<< " File position at start is "
|
|
<< pos_at_start << ", currently " << is.tellg();
|
|
}
|
|
|
|
|
|
// Constructor... note that this is not const-safe as it would
|
|
// be quite complicated to implement a "const SubMatrix" class that
|
|
// would not allow its contents to be changed.
|
|
template<typename Real>
|
|
SubMatrix<Real>::SubMatrix(const MatrixBase<Real> &M,
|
|
const MatrixIndexT ro,
|
|
const MatrixIndexT r,
|
|
const MatrixIndexT co,
|
|
const MatrixIndexT c) {
|
|
if (r == 0 || c == 0) {
|
|
// we support the empty sub-matrix as a special case.
|
|
KALDI_ASSERT(c == 0 && r == 0);
|
|
this->data_ = NULL;
|
|
this->num_cols_ = 0;
|
|
this->num_rows_ = 0;
|
|
this->stride_ = 0;
|
|
return;
|
|
}
|
|
KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(ro) <
|
|
static_cast<UnsignedMatrixIndexT>(M.num_rows_) &&
|
|
static_cast<UnsignedMatrixIndexT>(co) <
|
|
static_cast<UnsignedMatrixIndexT>(M.num_cols_) &&
|
|
static_cast<UnsignedMatrixIndexT>(r) <=
|
|
static_cast<UnsignedMatrixIndexT>(M.num_rows_ - ro) &&
|
|
static_cast<UnsignedMatrixIndexT>(c) <=
|
|
static_cast<UnsignedMatrixIndexT>(M.num_cols_ - co));
|
|
// point to the begining of window
|
|
MatrixBase<Real>::num_rows_ = r;
|
|
MatrixBase<Real>::num_cols_ = c;
|
|
MatrixBase<Real>::stride_ = M.Stride();
|
|
MatrixBase<Real>::data_ = M.Data_workaround() +
|
|
static_cast<size_t>(co) +
|
|
static_cast<size_t>(ro) * static_cast<size_t>(M.Stride());
|
|
}
|
|
|
|
|
|
template<typename Real>
|
|
SubMatrix<Real>::SubMatrix(Real *data,
|
|
MatrixIndexT num_rows,
|
|
MatrixIndexT num_cols,
|
|
MatrixIndexT stride):
|
|
MatrixBase<Real>(data, num_cols, num_rows, stride) { // caution: reversed order!
|
|
if (data == NULL) {
|
|
KALDI_ASSERT(num_rows * num_cols == 0);
|
|
this->num_rows_ = 0;
|
|
this->num_cols_ = 0;
|
|
this->stride_ = 0;
|
|
} else {
|
|
KALDI_ASSERT(this->stride_ >= this->num_cols_);
|
|
}
|
|
}
|
|
|
|
/*
|
|
template<typename Real>
|
|
void MatrixBase<Real>::Add(const Real alpha) {
|
|
Real *data = data_;
|
|
MatrixIndexT stride = stride_;
|
|
for (MatrixIndexT r = 0; r < num_rows_; r++)
|
|
for (MatrixIndexT c = 0; c < num_cols_; c++)
|
|
data[c + stride*r] += alpha;
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::AddToDiag(const Real alpha) {
|
|
Real *data = data_;
|
|
MatrixIndexT this_stride = stride_ + 1,
|
|
num_to_add = std::min(num_rows_, num_cols_);
|
|
for (MatrixIndexT r = 0; r < num_to_add; r++)
|
|
data[r * this_stride] += alpha;
|
|
}
|
|
|
|
|
|
template<typename Real>
|
|
Real MatrixBase<Real>::Cond() const {
|
|
KALDI_ASSERT(num_rows_ > 0&&num_cols_ > 0);
|
|
Vector<Real> singular_values(std::min(num_rows_, num_cols_));
|
|
Svd(&singular_values); // Get singular values...
|
|
Real min = singular_values(0), max = singular_values(0); // both absolute values...
|
|
for (MatrixIndexT i = 1;i < singular_values.Dim();i++) {
|
|
min = std::min((Real)std::abs(singular_values(i)), min); max = std::max((Real)std::abs(singular_values(i)), max);
|
|
}
|
|
if (min > 0) return max/min;
|
|
else return std::numeric_limits<Real>::infinity();
|
|
}
|
|
|
|
template<typename Real>
|
|
Real MatrixBase<Real>::Trace(bool check_square) const {
|
|
KALDI_ASSERT(!check_square || num_rows_ == num_cols_);
|
|
Real ans = 0.0;
|
|
for (MatrixIndexT r = 0;r < std::min(num_rows_, num_cols_);r++) ans += data_ [r + stride_*r];
|
|
return ans;
|
|
}
|
|
|
|
template<typename Real>
|
|
Real MatrixBase<Real>::Max() const {
|
|
KALDI_ASSERT(num_rows_ > 0 && num_cols_ > 0);
|
|
Real ans= *data_;
|
|
for (MatrixIndexT r = 0; r < num_rows_; r++)
|
|
for (MatrixIndexT c = 0; c < num_cols_; c++)
|
|
if (data_[c + stride_*r] > ans)
|
|
ans = data_[c + stride_*r];
|
|
return ans;
|
|
}
|
|
|
|
template<typename Real>
|
|
Real MatrixBase<Real>::Min() const {
|
|
KALDI_ASSERT(num_rows_ > 0 && num_cols_ > 0);
|
|
Real ans= *data_;
|
|
for (MatrixIndexT r = 0; r < num_rows_; r++)
|
|
for (MatrixIndexT c = 0; c < num_cols_; c++)
|
|
if (data_[c + stride_*r] < ans)
|
|
ans = data_[c + stride_*r];
|
|
return ans;
|
|
}
|
|
|
|
|
|
|
|
template <typename Real>
|
|
void MatrixBase<Real>::AddMatMatMat(Real alpha,
|
|
const MatrixBase<Real> &A, MatrixTransposeType transA,
|
|
const MatrixBase<Real> &B, MatrixTransposeType transB,
|
|
const MatrixBase<Real> &C, MatrixTransposeType transC,
|
|
Real beta) {
|
|
// Note on time taken with different orders of computation. Assume not transposed in this /
|
|
// discussion. Firstly, normalize expressions using A.NumCols == B.NumRows and B.NumCols == C.NumRows, prefer
|
|
// rows where there is a choice.
|
|
// time taken for (AB) is: A.NumRows*B.NumRows*C.Rows
|
|
// time taken for (AB)C is A.NumRows*C.NumRows*C.Cols
|
|
// so this order is A.NumRows*B.NumRows*C.NumRows + A.NumRows*C.NumRows*C.NumCols.
|
|
|
|
// time taken for (BC) is: B.NumRows*C.NumRows*C.Cols
|
|
// time taken for A(BC) is: A.NumRows*B.NumRows*C.Cols
|
|
// so this order is B.NumRows*C.NumRows*C.NumCols + A.NumRows*B.NumRows*C.Cols
|
|
|
|
MatrixIndexT ARows = A.num_rows_, ACols = A.num_cols_, BRows = B.num_rows_, BCols = B.num_cols_,
|
|
CRows = C.num_rows_, CCols = C.num_cols_;
|
|
if (transA == kTrans) std::swap(ARows, ACols);
|
|
if (transB == kTrans) std::swap(BRows, BCols);
|
|
if (transC == kTrans) std::swap(CRows, CCols);
|
|
|
|
MatrixIndexT AB_C_time = ARows*BRows*CRows + ARows*CRows*CCols;
|
|
MatrixIndexT A_BC_time = BRows*CRows*CCols + ARows*BRows*CCols;
|
|
|
|
if (AB_C_time < A_BC_time) {
|
|
Matrix<Real> AB(ARows, BCols);
|
|
AB.AddMatMat(1.0, A, transA, B, transB, 0.0); // AB = A * B.
|
|
(*this).AddMatMat(alpha, AB, kNoTrans, C, transC, beta);
|
|
} else {
|
|
Matrix<Real> BC(BRows, CCols);
|
|
BC.AddMatMat(1.0, B, transB, C, transC, 0.0); // BC = B * C.
|
|
(*this).AddMatMat(alpha, A, transA, BC, kNoTrans, beta);
|
|
}
|
|
}
|
|
|
|
|
|
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::DestructiveSvd(VectorBase<Real> *s, MatrixBase<Real> *U, MatrixBase<Real> *Vt) {
|
|
// Svd, *this = U*diag(s)*Vt.
|
|
// With (*this).num_rows_ == m, (*this).num_cols_ == n,
|
|
// Support only skinny Svd with m>=n (NumRows>=NumCols), and zero sizes for U and Vt mean
|
|
// we do not want that output. We expect that s.Dim() == m,
|
|
// U is either 0 by 0 or m by n, and rv is either 0 by 0 or n by n.
|
|
// Throws exception on error.
|
|
|
|
KALDI_ASSERT(num_rows_>=num_cols_ && "Svd requires that #rows by >= #cols."); // For compatibility with JAMA code.
|
|
KALDI_ASSERT(s->Dim() == num_cols_); // s should be the smaller dim.
|
|
KALDI_ASSERT(U == NULL || (U->num_rows_ == num_rows_&&U->num_cols_ == num_cols_));
|
|
KALDI_ASSERT(Vt == NULL || (Vt->num_rows_ == num_cols_&&Vt->num_cols_ == num_cols_));
|
|
|
|
Real prescale = 1.0;
|
|
if ( std::abs((*this)(0, 0) ) < 1.0e-30) { // Very tiny value... can cause problems in Svd.
|
|
Real max_elem = LargestAbsElem();
|
|
if (max_elem != 0) {
|
|
prescale = 1.0 / max_elem;
|
|
if (std::abs(prescale) == std::numeric_limits<Real>::infinity()) { prescale = 1.0e+40; }
|
|
(*this).Scale(prescale);
|
|
}
|
|
}
|
|
|
|
#if !defined(HAVE_ATLAS) && !defined(USE_KALDI_SVD)
|
|
// "S" == skinny Svd (only one we support because of compatibility with Jama one which is only skinny),
|
|
// "N"== no eigenvectors wanted.
|
|
LapackGesvd(s, U, Vt);
|
|
#else
|
|
/* if (num_rows_ > 1 && num_cols_ > 1 && (*this)(0, 0) == (*this)(1, 1)
|
|
&& Max() == Min() && (*this)(0, 0) != 0.0) { // special case that JamaSvd sometimes crashes on.
|
|
KALDI_WARN << "Jama SVD crashes on this type of matrix, perturbing it to prevent crash.";
|
|
for(int32 i = 0; i < NumRows(); i++)
|
|
(*this)(i, i) *= 1.00001;
|
|
}*/
|
|
// bool ans = JamaSvd(s, U, Vt);
|
|
//if (Vt != NULL) Vt->Transpose(); // possibly to do: change this and also the transpose inside the JamaSvd routine. note, Vt is square.
|
|
//if (!ans) {
|
|
//KALDI_ERR << "Error doing Svd"; // This one will be caught.
|
|
//}
|
|
//#endif
|
|
//if (prescale != 1.0) s->Scale(1.0/prescale);
|
|
//}
|
|
/*
|
|
template<typename Real>
|
|
void MatrixBase<Real>::Svd(VectorBase<Real> *s, MatrixBase<Real> *U, MatrixBase<Real> *Vt) const {
|
|
try {
|
|
if (num_rows_ >= num_cols_) {
|
|
Matrix<Real> tmp(*this);
|
|
tmp.DestructiveSvd(s, U, Vt);
|
|
} else {
|
|
Matrix<Real> tmp(*this, kTrans); // transpose of *this.
|
|
// rVt will have different dim so cannot transpose in-place --> use a temp matrix.
|
|
Matrix<Real> Vt_Trans(Vt ? Vt->num_cols_ : 0, Vt ? Vt->num_rows_ : 0);
|
|
// U will be transpose
|
|
tmp.DestructiveSvd(s, Vt ? &Vt_Trans : NULL, U);
|
|
if (U) U->Transpose();
|
|
if (Vt) Vt->CopyFromMat(Vt_Trans, kTrans); // copy with transpose.
|
|
}
|
|
} catch (...) {
|
|
KALDI_ERR << "Error doing Svd (did not converge), first part of matrix is\n"
|
|
<< SubMatrix<Real>(*this, 0, std::min((MatrixIndexT)10, num_rows_),
|
|
0, std::min((MatrixIndexT)10, num_cols_))
|
|
<< ", min and max are: " << Min() << ", " << Max();
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
bool MatrixBase<Real>::IsSymmetric(Real cutoff) const {
|
|
MatrixIndexT R = num_rows_, C = num_cols_;
|
|
if (R != C) return false;
|
|
Real bad_sum = 0.0, good_sum = 0.0;
|
|
for (MatrixIndexT i = 0;i < R;i++) {
|
|
for (MatrixIndexT j = 0;j < i;j++) {
|
|
Real a = (*this)(i, j), b = (*this)(j, i), avg = 0.5*(a+b), diff = 0.5*(a-b);
|
|
good_sum += std::abs(avg); bad_sum += std::abs(diff);
|
|
}
|
|
good_sum += std::abs((*this)(i, i));
|
|
}
|
|
if (bad_sum > cutoff*good_sum) return false;
|
|
return true;
|
|
}
|
|
|
|
template<typename Real>
|
|
bool MatrixBase<Real>::IsDiagonal(Real cutoff) const{
|
|
MatrixIndexT R = num_rows_, C = num_cols_;
|
|
Real bad_sum = 0.0, good_sum = 0.0;
|
|
for (MatrixIndexT i = 0;i < R;i++) {
|
|
for (MatrixIndexT j = 0;j < C;j++) {
|
|
if (i == j) good_sum += std::abs((*this)(i, j));
|
|
else bad_sum += std::abs((*this)(i, j));
|
|
}
|
|
}
|
|
return (!(bad_sum > good_sum * cutoff));
|
|
}
|
|
|
|
// This does nothing, it's designed to trigger Valgrind errors
|
|
// if any memory is uninitialized.
|
|
template<typename Real>
|
|
void MatrixBase<Real>::TestUninitialized() const {
|
|
MatrixIndexT R = num_rows_, C = num_cols_, positive = 0;
|
|
for (MatrixIndexT i = 0; i < R; i++)
|
|
for (MatrixIndexT j = 0; j < C; j++)
|
|
if ((*this)(i, j) > 0.0) positive++;
|
|
if (positive > R * C)
|
|
KALDI_ERR << "Error....";
|
|
}
|
|
|
|
|
|
template<typename Real>
|
|
bool MatrixBase<Real>::IsUnit(Real cutoff) const {
|
|
MatrixIndexT R = num_rows_, C = num_cols_;
|
|
Real bad_max = 0.0;
|
|
for (MatrixIndexT i = 0; i < R;i++)
|
|
for (MatrixIndexT j = 0; j < C;j++)
|
|
bad_max = std::max(bad_max, static_cast<Real>(std::abs( (*this)(i, j) - (i == j?1.0:0.0))));
|
|
return (bad_max <= cutoff);
|
|
}
|
|
|
|
template<typename Real>
|
|
bool MatrixBase<Real>::IsZero(Real cutoff)const {
|
|
MatrixIndexT R = num_rows_, C = num_cols_;
|
|
Real bad_max = 0.0;
|
|
for (MatrixIndexT i = 0;i < R;i++)
|
|
for (MatrixIndexT j = 0;j < C;j++)
|
|
bad_max = std::max(bad_max, static_cast<Real>(std::abs( (*this)(i, j) )));
|
|
return (bad_max <= cutoff);
|
|
}
|
|
|
|
template<typename Real>
|
|
Real MatrixBase<Real>::FrobeniusNorm() const{
|
|
return std::sqrt(TraceMatMat(*this, *this, kTrans));
|
|
}
|
|
|
|
template<typename Real>
|
|
bool MatrixBase<Real>::ApproxEqual(const MatrixBase<Real> &other, float tol) const {
|
|
if (num_rows_ != other.num_rows_ || num_cols_ != other.num_cols_)
|
|
KALDI_ERR << "ApproxEqual: size mismatch.";
|
|
Matrix<Real> tmp(*this);
|
|
tmp.AddMat(-1.0, other);
|
|
return (tmp.FrobeniusNorm() <= static_cast<Real>(tol) *
|
|
this->FrobeniusNorm());
|
|
}
|
|
|
|
template<typename Real>
|
|
bool MatrixBase<Real>::Equal(const MatrixBase<Real> &other) const {
|
|
if (num_rows_ != other.num_rows_ || num_cols_ != other.num_cols_)
|
|
KALDI_ERR << "Equal: size mismatch.";
|
|
for (MatrixIndexT i = 0; i < num_rows_; i++)
|
|
for (MatrixIndexT j = 0; j < num_cols_; j++)
|
|
if ( (*this)(i, j) != other(i, j))
|
|
return false;
|
|
return true;
|
|
}
|
|
|
|
|
|
template<typename Real>
|
|
Real MatrixBase<Real>::LargestAbsElem() const{
|
|
MatrixIndexT R = num_rows_, C = num_cols_;
|
|
Real largest = 0.0;
|
|
for (MatrixIndexT i = 0;i < R;i++)
|
|
for (MatrixIndexT j = 0;j < C;j++)
|
|
largest = std::max(largest, (Real)std::abs((*this)(i, j)));
|
|
return largest;
|
|
}
|
|
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::OrthogonalizeRows() {
|
|
KALDI_ASSERT(NumRows() <= NumCols());
|
|
MatrixIndexT num_rows = num_rows_;
|
|
for (MatrixIndexT i = 0; i < num_rows; i++) {
|
|
int32 counter = 0;
|
|
while (1) {
|
|
Real start_prod = VecVec(this->Row(i), this->Row(i));
|
|
if (start_prod - start_prod != 0.0 || start_prod == 0.0) {
|
|
KALDI_WARN << "Self-product of row " << i << " of matrix is "
|
|
<< start_prod << ", randomizing.";
|
|
this->Row(i).SetRandn();
|
|
counter++;
|
|
continue; // loop again.
|
|
}
|
|
for (MatrixIndexT j = 0; j < i; j++) {
|
|
Real prod = VecVec(this->Row(i), this->Row(j));
|
|
this->Row(i).AddVec(-prod, this->Row(j));
|
|
}
|
|
Real end_prod = VecVec(this->Row(i), this->Row(i));
|
|
if (end_prod <= 0.01 * start_prod) { // We removed
|
|
// almost all of the vector during orthogonalization,
|
|
// so we have reason to doubt (for roundoff reasons)
|
|
// that it's still orthogonal to the other vectors.
|
|
// We need to orthogonalize again.
|
|
if (end_prod == 0.0) { // Row is exactly zero:
|
|
// generate random direction.
|
|
this->Row(i).SetRandn();
|
|
}
|
|
counter++;
|
|
if (counter > 100)
|
|
KALDI_ERR << "Loop detected while orthogalizing matrix.";
|
|
} else {
|
|
this->Row(i).Scale(1.0 / std::sqrt(end_prod));
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Uses Svd to compute the eigenvalue decomposition of a symmetric positive semidefinite
|
|
// matrix:
|
|
// (*this) = rU * diag(rs) * rU^T, with rU an orthogonal matrix so rU^{-1} = rU^T.
|
|
// Does this by computing svd (*this) = U diag(rs) V^T ... answer is just U diag(rs) U^T.
|
|
// Throws exception if this failed to within supplied precision (typically because *this was not
|
|
// symmetric positive definite).
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::SymPosSemiDefEig(VectorBase<Real> *rs, MatrixBase<Real> *rU, Real check_thresh) // e.g. check_thresh = 0.001
|
|
{
|
|
const MatrixIndexT D = num_rows_;
|
|
|
|
KALDI_ASSERT(num_rows_ == num_cols_);
|
|
KALDI_ASSERT(IsSymmetric() && "SymPosSemiDefEig: expecting input to be symmetrical.");
|
|
KALDI_ASSERT(rU->num_rows_ == D && rU->num_cols_ == D && rs->Dim() == D);
|
|
|
|
Matrix<Real> Vt(D, D);
|
|
Svd(rs, rU, &Vt);
|
|
|
|
// First just zero any singular values if the column of U and V do not have +ve dot product--
|
|
// this may mean we have small negative eigenvalues, and if we zero them the result will be closer to correct.
|
|
for (MatrixIndexT i = 0;i < D;i++) {
|
|
Real sum = 0.0;
|
|
for (MatrixIndexT j = 0;j < D;j++) sum += (*rU)(j, i) * Vt(i, j);
|
|
if (sum < 0.0) (*rs)(i) = 0.0;
|
|
}
|
|
|
|
{
|
|
Matrix<Real> tmpU(*rU); Vector<Real> tmps(*rs); tmps.ApplyPow(0.5);
|
|
tmpU.MulColsVec(tmps);
|
|
SpMatrix<Real> tmpThis(D);
|
|
tmpThis.AddMat2(1.0, tmpU, kNoTrans, 0.0);
|
|
Matrix<Real> tmpThisFull(tmpThis);
|
|
float new_norm = tmpThisFull.FrobeniusNorm();
|
|
float old_norm = (*this).FrobeniusNorm();
|
|
tmpThisFull.AddMat(-1.0, (*this));
|
|
|
|
if (!(old_norm == 0 && new_norm == 0)) {
|
|
float diff_norm = tmpThisFull.FrobeniusNorm();
|
|
if (std::abs(new_norm-old_norm) > old_norm*check_thresh || diff_norm > old_norm*check_thresh) {
|
|
KALDI_WARN << "SymPosSemiDefEig seems to have failed " << diff_norm << " !<< "
|
|
<< check_thresh << "*" << old_norm << ", maybe matrix was not "
|
|
<< "positive semi definite. Continuing anyway.";
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
template<typename Real>
|
|
Real MatrixBase<Real>::LogDet(Real *det_sign) const {
|
|
Real log_det;
|
|
Matrix<Real> tmp(*this);
|
|
tmp.Invert(&log_det, det_sign, false); // false== output not needed (saves some computation).
|
|
return log_det;
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::InvertDouble(Real *log_det, Real *det_sign,
|
|
bool inverse_needed) {
|
|
double log_det_tmp, det_sign_tmp;
|
|
Matrix<double> dmat(*this);
|
|
dmat.Invert(&log_det_tmp, &det_sign_tmp, inverse_needed);
|
|
if (inverse_needed) (*this).CopyFromMat(dmat);
|
|
if (log_det) *log_det = log_det_tmp;
|
|
if (det_sign) *det_sign = det_sign_tmp;
|
|
}
|
|
*/
|
|
|
|
//template<class Real>
|
|
//void MatrixBase<Real>::CopyFromMat(const CompressedMatrix &mat) {
|
|
//mat.CopyToMat(this);
|
|
//}
|
|
|
|
//template<class Real>
|
|
//Matrix<Real>::Matrix(const CompressedMatrix &M): MatrixBase<Real>() {
|
|
//Resize(M.NumRows(), M.NumCols(), kUndefined);
|
|
//M.CopyToMat(this);
|
|
//}
|
|
|
|
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::InvertElements() {
|
|
for (MatrixIndexT r = 0; r < num_rows_; r++) {
|
|
for (MatrixIndexT c = 0; c < num_cols_; c++) {
|
|
(*this)(r, c) = static_cast<Real>(1.0 / (*this)(r, c));
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
template<typename Real>
|
|
void MatrixBase<Real>::Transpose() {
|
|
KALDI_ASSERT(num_rows_ == num_cols_);
|
|
MatrixIndexT M = num_rows_;
|
|
for (MatrixIndexT i = 0;i < M;i++)
|
|
for (MatrixIndexT j = 0;j < i;j++) {
|
|
Real &a = (*this)(i, j), &b = (*this)(j, i);
|
|
std::swap(a, b);
|
|
}
|
|
}
|
|
|
|
|
|
template<typename Real>
|
|
void Matrix<Real>::Transpose() {
|
|
if (this->num_rows_ != this->num_cols_) {
|
|
Matrix<Real> tmp(*this, kTrans);
|
|
Resize(this->num_cols_, this->num_rows_);
|
|
this->CopyFromMat(tmp);
|
|
} else {
|
|
(static_cast<MatrixBase<Real>&>(*this)).Transpose();
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::Heaviside(const MatrixBase<Real> &src) {
|
|
KALDI_ASSERT(SameDim(*this, src));
|
|
MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
|
|
Real *row_data = data_;
|
|
const Real *src_row_data = src.Data();
|
|
for (MatrixIndexT row = 0; row < num_rows;
|
|
row++,row_data += stride_, src_row_data += src.stride_) {
|
|
for (MatrixIndexT col = 0; col < num_cols; col++)
|
|
row_data[col] = (src_row_data[col] > 0 ? 1.0 : 0.0);
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::Exp(const MatrixBase<Real> &src) {
|
|
KALDI_ASSERT(SameDim(*this, src));
|
|
MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
|
|
Real *row_data = data_;
|
|
const Real *src_row_data = src.Data();
|
|
for (MatrixIndexT row = 0; row < num_rows;
|
|
row++,row_data += stride_, src_row_data += src.stride_) {
|
|
for (MatrixIndexT col = 0; col < num_cols; col++)
|
|
row_data[col] = kaldi::Exp(src_row_data[col]);
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::Pow(const MatrixBase<Real> &src, Real power) {
|
|
KALDI_ASSERT(SameDim(*this, src));
|
|
MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
|
|
Real *row_data = data_;
|
|
const Real *src_row_data = src.Data();
|
|
for (MatrixIndexT row = 0; row < num_rows;
|
|
row++,row_data += stride_, src_row_data += src.stride_) {
|
|
for (MatrixIndexT col = 0; col < num_cols; col++) {
|
|
row_data[col] = pow(src_row_data[col], power);
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::PowAbs(const MatrixBase<Real> &src, Real power, bool include_sign) {
|
|
KALDI_ASSERT(SameDim(*this, src));
|
|
MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
|
|
Real *row_data = data_;
|
|
const Real *src_row_data = src.Data();
|
|
for (MatrixIndexT row = 0; row < num_rows;
|
|
row++,row_data += stride_, src_row_data += src.stride_) {
|
|
for (MatrixIndexT col = 0; col < num_cols; col ++) {
|
|
if (include_sign == true && src_row_data[col] < 0) {
|
|
row_data[col] = -pow(std::abs(src_row_data[col]), power);
|
|
} else {
|
|
row_data[col] = pow(std::abs(src_row_data[col]), power);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::Floor(const MatrixBase<Real> &src, Real floor_val) {
|
|
KALDI_ASSERT(SameDim(*this, src));
|
|
MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
|
|
Real *row_data = data_;
|
|
const Real *src_row_data = src.Data();
|
|
for (MatrixIndexT row = 0; row < num_rows;
|
|
row++,row_data += stride_, src_row_data += src.stride_) {
|
|
for (MatrixIndexT col = 0; col < num_cols; col++)
|
|
row_data[col] = (src_row_data[col] < floor_val ? floor_val : src_row_data[col]);
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::Ceiling(const MatrixBase<Real> &src, Real ceiling_val) {
|
|
KALDI_ASSERT(SameDim(*this, src));
|
|
MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
|
|
Real *row_data = data_;
|
|
const Real *src_row_data = src.Data();
|
|
for (MatrixIndexT row = 0; row < num_rows;
|
|
row++,row_data += stride_, src_row_data += src.stride_) {
|
|
for (MatrixIndexT col = 0; col < num_cols; col++)
|
|
row_data[col] = (src_row_data[col] > ceiling_val ? ceiling_val : src_row_data[col]);
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::Log(const MatrixBase<Real> &src) {
|
|
KALDI_ASSERT(SameDim(*this, src));
|
|
MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
|
|
Real *row_data = data_;
|
|
const Real *src_row_data = src.Data();
|
|
for (MatrixIndexT row = 0; row < num_rows;
|
|
row++,row_data += stride_, src_row_data += src.stride_) {
|
|
for (MatrixIndexT col = 0; col < num_cols; col++)
|
|
row_data[col] = kaldi::Log(src_row_data[col]);
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::ExpSpecial(const MatrixBase<Real> &src) {
|
|
KALDI_ASSERT(SameDim(*this, src));
|
|
MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
|
|
Real *row_data = data_;
|
|
const Real *src_row_data = src.Data();
|
|
for (MatrixIndexT row = 0; row < num_rows;
|
|
row++,row_data += stride_, src_row_data += src.stride_) {
|
|
for (MatrixIndexT col = 0; col < num_cols; col++)
|
|
row_data[col] = (src_row_data[col] < Real(0) ? kaldi::Exp(src_row_data[col]) : (src_row_data[col] + Real(1)));
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::ExpLimited(const MatrixBase<Real> &src, Real lower_limit, Real upper_limit) {
|
|
KALDI_ASSERT(SameDim(*this, src));
|
|
MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
|
|
Real *row_data = data_;
|
|
const Real *src_row_data = src.Data();
|
|
for (MatrixIndexT row = 0; row < num_rows;
|
|
row++,row_data += stride_, src_row_data += src.stride_) {
|
|
for (MatrixIndexT col = 0; col < num_cols; col++) {
|
|
const Real x = src_row_data[col];
|
|
if (!(x >= lower_limit))
|
|
row_data[col] = kaldi::Exp(lower_limit);
|
|
else if (x > upper_limit)
|
|
row_data[col] = kaldi::Exp(upper_limit);
|
|
else
|
|
row_data[col] = kaldi::Exp(x);
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
bool MatrixBase<Real>::Power(Real power) {
|
|
KALDI_ASSERT(num_rows_ > 0 && num_rows_ == num_cols_);
|
|
MatrixIndexT n = num_rows_;
|
|
Matrix<Real> P(n, n);
|
|
Vector<Real> re(n), im(n);
|
|
this->Eig(&P, &re, &im);
|
|
// Now attempt to take the complex eigenvalues to this power.
|
|
for (MatrixIndexT i = 0; i < n; i++)
|
|
if (!AttemptComplexPower(&(re(i)), &(im(i)), power))
|
|
return false; // e.g. real and negative, or zero, eigenvalues.
|
|
|
|
Matrix<Real> D(n, n); // D to the power.
|
|
CreateEigenvalueMatrix(re, im, &D);
|
|
|
|
Matrix<Real> tmp(n, n); // P times D
|
|
tmp.AddMatMat(1.0, P, kNoTrans, D, kNoTrans, 0.0); // tmp := P*D
|
|
P.Invert();
|
|
// next line is: *this = tmp * P^{-1} = P * D * P^{-1}
|
|
(*this).AddMatMat(1.0, tmp, kNoTrans, P, kNoTrans, 0.0);
|
|
return true;
|
|
}
|
|
*/
|
|
template<typename Real>
|
|
void Matrix<Real>::Swap(Matrix<Real> *other) {
|
|
std::swap(this->data_, other->data_);
|
|
std::swap(this->num_cols_, other->num_cols_);
|
|
std::swap(this->num_rows_, other->num_rows_);
|
|
std::swap(this->stride_, other->stride_);
|
|
}
|
|
/*
|
|
// Repeating this comment that appeared in the header:
|
|
// 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.
|
|
//
|
|
// By making the pointer arguments non-NULL or NULL, the user can choose to take
|
|
// not to take the eigenvalues directly, and/or the matrix D which is block-diagonal
|
|
// with 2x2 blocks.
|
|
template<typename Real>
|
|
void MatrixBase<Real>::Eig(MatrixBase<Real> *P,
|
|
VectorBase<Real> *r,
|
|
VectorBase<Real> *i) const {
|
|
EigenvalueDecomposition<Real> eig(*this);
|
|
if (P) eig.GetV(P);
|
|
if (r) eig.GetRealEigenvalues(r);
|
|
if (i) eig.GetImagEigenvalues(i);
|
|
}
|
|
|
|
|
|
// Begin non-member function definitions.
|
|
|
|
// /**
|
|
// * @brief Extension of the HTK header
|
|
// */
|
|
// struct HtkHeaderExt
|
|
// {
|
|
// INT_32 mHeaderSize;
|
|
// INT_32 mVersion;
|
|
// INT_32 mSampSize;
|
|
// };
|
|
/*
|
|
template<typename Real>
|
|
bool ReadHtk(std::istream &is, Matrix<Real> *M_ptr, HtkHeader *header_ptr)
|
|
{
|
|
// check instantiated with double or float.
|
|
KALDI_ASSERT_IS_FLOATING_TYPE(Real);
|
|
Matrix<Real> &M = *M_ptr;
|
|
HtkHeader htk_hdr;
|
|
|
|
// TODO(arnab): this fails if the HTK file has CRC cheksum or is compressed.
|
|
is.read((char*)&htk_hdr, sizeof(htk_hdr)); // we're being really POSIX here!
|
|
if (is.fail()) {
|
|
KALDI_WARN << "Could not read header from HTK feature file ";
|
|
return false;
|
|
}
|
|
|
|
KALDI_SWAP4(htk_hdr.mNSamples);
|
|
KALDI_SWAP4(htk_hdr.mSamplePeriod);
|
|
KALDI_SWAP2(htk_hdr.mSampleSize);
|
|
KALDI_SWAP2(htk_hdr.mSampleKind);
|
|
|
|
bool has_checksum = false;
|
|
{
|
|
// See HParm.h in HTK code for sources of these things.
|
|
enum BaseParmKind{
|
|
Waveform, Lpc, Lprefc, Lpcepstra, Lpdelcep,
|
|
Irefc, Mfcc, Fbank, Melspec, User, Discrete, Plp, Anon };
|
|
|
|
const int32 IsCompressed = 02000, HasChecksum = 010000, HasVq = 040000,
|
|
Problem = IsCompressed | HasVq;
|
|
int32 base_parm = htk_hdr.mSampleKind & (077);
|
|
has_checksum = (base_parm & HasChecksum) != 0;
|
|
htk_hdr.mSampleKind &= ~HasChecksum; // We don't support writing with
|
|
// checksum so turn it off.
|
|
if (htk_hdr.mSampleKind & Problem)
|
|
KALDI_ERR << "Code to read HTK features does not support compressed "
|
|
"features, or features with VQ.";
|
|
if (base_parm == Waveform || base_parm == Irefc || base_parm == Discrete)
|
|
KALDI_ERR << "Attempting to read HTK features from unsupported type "
|
|
"(e.g. waveform or discrete features.";
|
|
}
|
|
|
|
KALDI_VLOG(3) << "HTK header: Num Samples: " << htk_hdr.mNSamples
|
|
<< "; Sample period: " << htk_hdr.mSamplePeriod
|
|
<< "; Sample size: " << htk_hdr.mSampleSize
|
|
<< "; Sample kind: " << htk_hdr.mSampleKind;
|
|
|
|
M.Resize(htk_hdr.mNSamples, htk_hdr.mSampleSize / sizeof(float));
|
|
|
|
MatrixIndexT i;
|
|
MatrixIndexT j;
|
|
if (sizeof(Real) == sizeof(float)) {
|
|
for (i = 0; i< M.NumRows(); i++) {
|
|
is.read((char*)M.RowData(i), sizeof(float)*M.NumCols());
|
|
if (is.fail()) {
|
|
KALDI_WARN << "Could not read data from HTK feature file ";
|
|
return false;
|
|
}
|
|
if (MachineIsLittleEndian()) {
|
|
MatrixIndexT C = M.NumCols();
|
|
for (j = 0; j < C; j++) {
|
|
KALDI_SWAP4((M(i, j))); // The HTK standard is big-endian!
|
|
}
|
|
}
|
|
}
|
|
} else {
|
|
float *pmem = new float[M.NumCols()];
|
|
for (i = 0; i < M.NumRows(); i++) {
|
|
is.read((char*)pmem, sizeof(float)*M.NumCols());
|
|
if (is.fail()) {
|
|
KALDI_WARN << "Could not read data from HTK feature file ";
|
|
delete [] pmem;
|
|
return false;
|
|
}
|
|
MatrixIndexT C = M.NumCols();
|
|
for (j = 0; j < C; j++) {
|
|
if (MachineIsLittleEndian()) // HTK standard is big-endian!
|
|
KALDI_SWAP4(pmem[j]);
|
|
M(i, j) = static_cast<Real>(pmem[j]);
|
|
}
|
|
}
|
|
delete [] pmem;
|
|
}
|
|
if (header_ptr) *header_ptr = htk_hdr;
|
|
if (has_checksum) {
|
|
int16 checksum;
|
|
is.read((char*)&checksum, sizeof(checksum));
|
|
if (is.fail())
|
|
KALDI_WARN << "Could not read checksum from HTK feature file ";
|
|
// We ignore the checksum.
|
|
}
|
|
return true;
|
|
}
|
|
|
|
|
|
template
|
|
bool ReadHtk(std::istream &is, Matrix<float> *M, HtkHeader *header_ptr);
|
|
|
|
template
|
|
bool ReadHtk(std::istream &is, Matrix<double> *M, HtkHeader *header_ptr);
|
|
|
|
template<typename Real>
|
|
bool WriteHtk(std::ostream &os, const MatrixBase<Real> &M, HtkHeader htk_hdr) // header may be derived from a previous call to ReadHtk. Must be in binary mode.
|
|
{
|
|
KALDI_ASSERT(M.NumRows() == static_cast<MatrixIndexT>(htk_hdr.mNSamples));
|
|
KALDI_ASSERT(M.NumCols() == static_cast<MatrixIndexT>(htk_hdr.mSampleSize) /
|
|
static_cast<MatrixIndexT>(sizeof(float)));
|
|
|
|
KALDI_SWAP4(htk_hdr.mNSamples);
|
|
KALDI_SWAP4(htk_hdr.mSamplePeriod);
|
|
KALDI_SWAP2(htk_hdr.mSampleSize);
|
|
KALDI_SWAP2(htk_hdr.mSampleKind);
|
|
|
|
os.write((char*)&htk_hdr, sizeof(htk_hdr));
|
|
if (os.fail()) goto bad;
|
|
|
|
MatrixIndexT i;
|
|
MatrixIndexT j;
|
|
if (sizeof(Real) == sizeof(float) && !MachineIsLittleEndian()) {
|
|
for (i = 0; i< M.NumRows(); i++) { // Unlikely to reach here ever!
|
|
os.write((char*)M.RowData(i), sizeof(float)*M.NumCols());
|
|
if (os.fail()) goto bad;
|
|
}
|
|
} else {
|
|
float *pmem = new float[M.NumCols()];
|
|
|
|
for (i = 0; i < M.NumRows(); i++) {
|
|
const Real *rowData = M.RowData(i);
|
|
for (j = 0;j < M.NumCols();j++)
|
|
pmem[j] = static_cast<float> ( rowData[j] );
|
|
if (MachineIsLittleEndian())
|
|
for (j = 0;j < M.NumCols();j++)
|
|
KALDI_SWAP4(pmem[j]);
|
|
os.write((char*)pmem, sizeof(float)*M.NumCols());
|
|
if (os.fail()) {
|
|
delete [] pmem;
|
|
goto bad;
|
|
}
|
|
}
|
|
delete [] pmem;
|
|
}
|
|
return true;
|
|
bad:
|
|
KALDI_WARN << "Could not write to HTK feature file ";
|
|
return false;
|
|
}
|
|
|
|
template
|
|
bool WriteHtk(std::ostream &os, const MatrixBase<float> &M, HtkHeader htk_hdr);
|
|
|
|
template
|
|
bool WriteHtk(std::ostream &os, const MatrixBase<double> &M, HtkHeader htk_hdr);
|
|
|
|
template<class Real>
|
|
bool WriteSphinx(std::ostream &os, const MatrixBase<Real> &M)
|
|
{
|
|
// CMUSphinx mfc file header contains count of the floats, followed
|
|
// by the data in float little endian format.
|
|
|
|
int size = M.NumRows() * M.NumCols();
|
|
os.write((char*)&size, sizeof(int));
|
|
if (os.fail()) goto bad;
|
|
|
|
MatrixIndexT i;
|
|
MatrixIndexT j;
|
|
if (sizeof(Real) == sizeof(float) && MachineIsLittleEndian()) {
|
|
for (i = 0; i< M.NumRows(); i++) { // Unlikely to reach here ever!
|
|
os.write((char*)M.RowData(i), sizeof(float)*M.NumCols());
|
|
if (os.fail()) goto bad;
|
|
}
|
|
} else {
|
|
float *pmem = new float[M.NumCols()];
|
|
|
|
for (i = 0; i < M.NumRows(); i++) {
|
|
const Real *rowData = M.RowData(i);
|
|
for (j = 0;j < M.NumCols();j++)
|
|
pmem[j] = static_cast<float> ( rowData[j] );
|
|
if (!MachineIsLittleEndian())
|
|
for (j = 0;j < M.NumCols();j++)
|
|
KALDI_SWAP4(pmem[j]);
|
|
os.write((char*)pmem, sizeof(float)*M.NumCols());
|
|
if (os.fail()) {
|
|
delete [] pmem;
|
|
goto bad;
|
|
}
|
|
}
|
|
delete [] pmem;
|
|
}
|
|
return true;
|
|
bad:
|
|
KALDI_WARN << "Could not write to Sphinx feature file";
|
|
return false;
|
|
}
|
|
|
|
template
|
|
bool WriteSphinx(std::ostream &os, const MatrixBase<float> &M);
|
|
|
|
template
|
|
bool WriteSphinx(std::ostream &os, const MatrixBase<double> &M);
|
|
|
|
template <typename Real>
|
|
Real TraceMatMatMat(const MatrixBase<Real> &A, MatrixTransposeType transA,
|
|
const MatrixBase<Real> &B, MatrixTransposeType transB,
|
|
const MatrixBase<Real> &C, MatrixTransposeType transC) {
|
|
MatrixIndexT ARows = A.NumRows(), ACols = A.NumCols(), BRows = B.NumRows(), BCols = B.NumCols(),
|
|
CRows = C.NumRows(), CCols = C.NumCols();
|
|
if (transA == kTrans) std::swap(ARows, ACols);
|
|
if (transB == kTrans) std::swap(BRows, BCols);
|
|
if (transC == kTrans) std::swap(CRows, CCols);
|
|
KALDI_ASSERT( CCols == ARows && ACols == BRows && BCols == CRows && "TraceMatMatMat: args have mismatched dimensions.");
|
|
if (ARows*BCols < std::min(BRows*CCols, CRows*ACols)) {
|
|
Matrix<Real> AB(ARows, BCols);
|
|
AB.AddMatMat(1.0, A, transA, B, transB, 0.0); // AB = A * B.
|
|
return TraceMatMat(AB, C, transC);
|
|
} else if ( BRows*CCols < CRows*ACols) {
|
|
Matrix<Real> BC(BRows, CCols);
|
|
BC.AddMatMat(1.0, B, transB, C, transC, 0.0); // BC = B * C.
|
|
return TraceMatMat(BC, A, transA);
|
|
} else {
|
|
Matrix<Real> CA(CRows, ACols);
|
|
CA.AddMatMat(1.0, C, transC, A, transA, 0.0); // CA = C * A
|
|
return TraceMatMat(CA, B, transB);
|
|
}
|
|
}
|
|
|
|
template
|
|
float TraceMatMatMat(const MatrixBase<float> &A, MatrixTransposeType transA,
|
|
const MatrixBase<float> &B, MatrixTransposeType transB,
|
|
const MatrixBase<float> &C, MatrixTransposeType transC);
|
|
|
|
template
|
|
double TraceMatMatMat(const MatrixBase<double> &A, MatrixTransposeType transA,
|
|
const MatrixBase<double> &B, MatrixTransposeType transB,
|
|
const MatrixBase<double> &C, MatrixTransposeType transC);
|
|
|
|
|
|
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) {
|
|
MatrixIndexT ARows = A.NumRows(), ACols = A.NumCols(), BRows = B.NumRows(), BCols = B.NumCols(),
|
|
CRows = C.NumRows(), CCols = C.NumCols(), DRows = D.NumRows(), DCols = D.NumCols();
|
|
if (transA == kTrans) std::swap(ARows, ACols);
|
|
if (transB == kTrans) std::swap(BRows, BCols);
|
|
if (transC == kTrans) std::swap(CRows, CCols);
|
|
if (transD == kTrans) std::swap(DRows, DCols);
|
|
KALDI_ASSERT( DCols == ARows && ACols == BRows && BCols == CRows && CCols == DRows && "TraceMatMatMat: args have mismatched dimensions.");
|
|
if (ARows*BCols < std::min(BRows*CCols, std::min(CRows*DCols, DRows*ACols))) {
|
|
Matrix<Real> AB(ARows, BCols);
|
|
AB.AddMatMat(1.0, A, transA, B, transB, 0.0); // AB = A * B.
|
|
return TraceMatMatMat(AB, kNoTrans, C, transC, D, transD);
|
|
} else if ((BRows*CCols) < std::min(CRows*DCols, DRows*ACols)) {
|
|
Matrix<Real> BC(BRows, CCols);
|
|
BC.AddMatMat(1.0, B, transB, C, transC, 0.0); // BC = B * C.
|
|
return TraceMatMatMat(BC, kNoTrans, D, transD, A, transA);
|
|
} else if (CRows*DCols < DRows*ACols) {
|
|
Matrix<Real> CD(CRows, DCols);
|
|
CD.AddMatMat(1.0, C, transC, D, transD, 0.0); // CD = C * D
|
|
return TraceMatMatMat(CD, kNoTrans, A, transA, B, transB);
|
|
} else {
|
|
Matrix<Real> DA(DRows, ACols);
|
|
DA.AddMatMat(1.0, D, transD, A, transA, 0.0); // DA = D * A
|
|
return TraceMatMatMat(DA, kNoTrans, B, transB, C, transC);
|
|
}
|
|
}
|
|
|
|
template
|
|
float TraceMatMatMatMat(const MatrixBase<float> &A, MatrixTransposeType transA,
|
|
const MatrixBase<float> &B, MatrixTransposeType transB,
|
|
const MatrixBase<float> &C, MatrixTransposeType transC,
|
|
const MatrixBase<float> &D, MatrixTransposeType transD);
|
|
|
|
template
|
|
double TraceMatMatMatMat(const MatrixBase<double> &A, MatrixTransposeType transA,
|
|
const MatrixBase<double> &B, MatrixTransposeType transB,
|
|
const MatrixBase<double> &C, MatrixTransposeType transC,
|
|
const MatrixBase<double> &D, MatrixTransposeType transD);
|
|
|
|
template<typename Real> void SortSvd(VectorBase<Real> *s, MatrixBase<Real> *U,
|
|
MatrixBase<Real> *Vt, bool sort_on_absolute_value) {
|
|
/// Makes sure the Svd is sorted (from greatest to least absolute value).
|
|
MatrixIndexT num_singval = s->Dim();
|
|
KALDI_ASSERT(U == NULL || U->NumCols() == num_singval);
|
|
KALDI_ASSERT(Vt == NULL || Vt->NumRows() == num_singval);
|
|
|
|
std::vector<std::pair<Real, MatrixIndexT> > vec(num_singval);
|
|
// negative because we want revese order.
|
|
for (MatrixIndexT d = 0; d < num_singval; d++) {
|
|
Real val = (*s)(d),
|
|
sort_val = -(sort_on_absolute_value ? std::abs(val) : val);
|
|
vec[d] = std::pair<Real, MatrixIndexT>(sort_val, d);
|
|
}
|
|
std::sort(vec.begin(), vec.end());
|
|
Vector<Real> s_copy(*s);
|
|
for (MatrixIndexT d = 0; d < num_singval; d++)
|
|
(*s)(d) = s_copy(vec[d].second);
|
|
if (U != NULL) {
|
|
Matrix<Real> Utmp(*U);
|
|
MatrixIndexT dim = Utmp.NumRows();
|
|
for (MatrixIndexT d = 0; d < num_singval; d++) {
|
|
MatrixIndexT oldidx = vec[d].second;
|
|
for (MatrixIndexT e = 0; e < dim; e++)
|
|
(*U)(e, d) = Utmp(e, oldidx);
|
|
}
|
|
}
|
|
if (Vt != NULL) {
|
|
Matrix<Real> Vttmp(*Vt);
|
|
for (MatrixIndexT d = 0; d < num_singval; d++)
|
|
(*Vt).Row(d).CopyFromVec(Vttmp.Row(vec[d].second));
|
|
}
|
|
}
|
|
|
|
template
|
|
void SortSvd(VectorBase<float> *s, MatrixBase<float> *U,
|
|
MatrixBase<float> *Vt, bool);
|
|
|
|
template
|
|
void SortSvd(VectorBase<double> *s, MatrixBase<double> *U,
|
|
MatrixBase<double> *Vt, bool);
|
|
|
|
template<typename Real>
|
|
void CreateEigenvalueMatrix(const VectorBase<Real> &re, const VectorBase<Real> &im,
|
|
MatrixBase<Real> *D) {
|
|
MatrixIndexT n = re.Dim();
|
|
KALDI_ASSERT(im.Dim() == n && D->NumRows() == n && D->NumCols() == n);
|
|
|
|
MatrixIndexT j = 0;
|
|
D->SetZero();
|
|
while (j < n) {
|
|
if (im(j) == 0) { // Real eigenvalue
|
|
(*D)(j, j) = re(j);
|
|
j++;
|
|
} else { // First of a complex pair
|
|
KALDI_ASSERT(j+1 < n && ApproxEqual(im(j+1), -im(j))
|
|
&& ApproxEqual(re(j+1), re(j)));
|
|
/// if (im(j) < 0.0) KALDI_WARN << "Negative first im part of pair"; // TEMP
|
|
Real lambda = re(j), mu = im(j);
|
|
// create 2x2 block [lambda, mu; -mu, lambda]
|
|
(*D)(j, j) = lambda;
|
|
(*D)(j, j+1) = mu;
|
|
(*D)(j+1, j) = -mu;
|
|
(*D)(j+1, j+1) = lambda;
|
|
j += 2;
|
|
}
|
|
}
|
|
}
|
|
|
|
template
|
|
void CreateEigenvalueMatrix(const VectorBase<float> &re, const VectorBase<float> &im,
|
|
MatrixBase<float> *D);
|
|
template
|
|
void CreateEigenvalueMatrix(const VectorBase<double> &re, const VectorBase<double> &im,
|
|
MatrixBase<double> *D);
|
|
|
|
|
|
|
|
template<typename Real>
|
|
bool AttemptComplexPower(Real *x_re, Real *x_im, Real power) {
|
|
// Used in Matrix<Real>::Power().
|
|
// Attempts to take the complex value x to the power "power",
|
|
// assuming that power is fractional (i.e. we don't treat integers as a
|
|
// special case). Returns false if this is not possible, either
|
|
// because x is negative and real (hence there is no obvious answer
|
|
// that is "closest to 1", and anyway this case does not make sense
|
|
// in the Matrix<Real>::Power() routine);
|
|
// or because power is negative, and x is zero.
|
|
|
|
// First solve for r and theta in
|
|
// x_re = r*cos(theta), x_im = r*sin(theta)
|
|
if (*x_re < 0.0 && *x_im == 0.0) return false; // can't do
|
|
// it for negative real values.
|
|
Real r = std::sqrt((*x_re * *x_re) + (*x_im * *x_im)); // r == radius.
|
|
if (power < 0.0 && r == 0.0) return false;
|
|
Real theta = std::atan2(*x_im, *x_re);
|
|
// Take the power.
|
|
r = std::pow(r, power);
|
|
theta *= power;
|
|
*x_re = r * std::cos(theta);
|
|
*x_im = r * std::sin(theta);
|
|
return true;
|
|
}
|
|
|
|
template
|
|
bool AttemptComplexPower(float *x_re, float *x_im, float power);
|
|
template
|
|
bool AttemptComplexPower(double *x_re, double *x_im, double power);
|
|
|
|
|
|
|
|
template <typename Real>
|
|
Real TraceMatMat(const MatrixBase<Real> &A,
|
|
const MatrixBase<Real> &B,
|
|
MatrixTransposeType trans) { // tr(A B), equivalent to sum of each element of A times same element in B'
|
|
MatrixIndexT aStride = A.stride_, bStride = B.stride_;
|
|
if (trans == kNoTrans) {
|
|
KALDI_ASSERT(A.NumRows() == B.NumCols() && A.NumCols() == B.NumRows());
|
|
Real ans = 0.0;
|
|
Real *adata = A.data_, *bdata = B.data_;
|
|
MatrixIndexT arows = A.NumRows(), acols = A.NumCols();
|
|
for (MatrixIndexT row = 0;row < arows;row++, adata+=aStride, bdata++)
|
|
ans += cblas_Xdot(acols, adata, 1, bdata, bStride);
|
|
return ans;
|
|
} else {
|
|
KALDI_ASSERT(A.NumRows() == B.NumRows() && A.NumCols() == B.NumCols());
|
|
Real ans = 0.0;
|
|
Real *adata = A.data_, *bdata = B.data_;
|
|
MatrixIndexT arows = A.NumRows(), acols = A.NumCols();
|
|
for (MatrixIndexT row = 0;row < arows;row++, adata+=aStride, bdata+=bStride)
|
|
ans += cblas_Xdot(acols, adata, 1, bdata, 1);
|
|
return ans;
|
|
}
|
|
}
|
|
|
|
|
|
// Instantiate the template above for float and double.
|
|
template
|
|
float TraceMatMat(const MatrixBase<float> &A,
|
|
const MatrixBase<float> &B,
|
|
MatrixTransposeType trans);
|
|
template
|
|
double TraceMatMat(const MatrixBase<double> &A,
|
|
const MatrixBase<double> &B,
|
|
MatrixTransposeType trans);
|
|
|
|
|
|
template<typename Real>
|
|
Real MatrixBase<Real>::LogSumExp(Real prune) const {
|
|
Real sum;
|
|
if (sizeof(sum) == 8) sum = kLogZeroDouble;
|
|
else sum = kLogZeroFloat;
|
|
Real max_elem = Max(), cutoff;
|
|
if (sizeof(Real) == 4) cutoff = max_elem + kMinLogDiffFloat;
|
|
else cutoff = max_elem + kMinLogDiffDouble;
|
|
if (prune > 0.0 && max_elem - prune > cutoff) // explicit pruning...
|
|
cutoff = max_elem - prune;
|
|
|
|
double sum_relto_max_elem = 0.0;
|
|
|
|
for (MatrixIndexT i = 0; i < num_rows_; i++) {
|
|
for (MatrixIndexT j = 0; j < num_cols_; j++) {
|
|
BaseFloat f = (*this)(i, j);
|
|
if (f >= cutoff)
|
|
sum_relto_max_elem += kaldi::Exp(f - max_elem);
|
|
}
|
|
}
|
|
return max_elem + kaldi::Log(sum_relto_max_elem);
|
|
}
|
|
|
|
template<typename Real>
|
|
Real MatrixBase<Real>::ApplySoftMax() {
|
|
Real max = this->Max(), sum = 0.0;
|
|
// the 'max' helps to get in good numeric range.
|
|
for (MatrixIndexT i = 0; i < num_rows_; i++)
|
|
for (MatrixIndexT j = 0; j < num_cols_; j++)
|
|
sum += ((*this)(i, j) = kaldi::Exp((*this)(i, j) - max));
|
|
this->Scale(1.0 / sum);
|
|
return max + kaldi::Log(sum);
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::Tanh(const MatrixBase<Real> &src) {
|
|
KALDI_ASSERT(SameDim(*this, src));
|
|
|
|
if (num_cols_ == stride_ && src.num_cols_ == src.stride_) {
|
|
SubVector<Real> src_vec(src.data_, num_rows_ * num_cols_),
|
|
dst_vec(this->data_, num_rows_ * num_cols_);
|
|
dst_vec.Tanh(src_vec);
|
|
} else {
|
|
for (MatrixIndexT r = 0; r < num_rows_; r++) {
|
|
SubVector<Real> src_vec(src, r), dest_vec(*this, r);
|
|
dest_vec.Tanh(src_vec);
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::SoftHinge(const MatrixBase<Real> &src) {
|
|
KALDI_ASSERT(SameDim(*this, src));
|
|
int32 num_rows = num_rows_, num_cols = num_cols_;
|
|
for (MatrixIndexT r = 0; r < num_rows; r++) {
|
|
Real *row_data = this->RowData(r);
|
|
const Real *src_row_data = src.RowData(r);
|
|
for (MatrixIndexT c = 0; c < num_cols; c++) {
|
|
Real x = src_row_data[c], y;
|
|
if (x > 10.0) y = x; // avoid exponentiating large numbers; function
|
|
// approaches y=x.
|
|
else y = Log1p(kaldi::Exp(x)); // these defined in kaldi-math.h
|
|
row_data[c] = y;
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::GroupPnorm(const MatrixBase<Real> &src, Real power) {
|
|
KALDI_ASSERT(src.NumCols() % this->NumCols() == 0 &&
|
|
src.NumRows() == this->NumRows());
|
|
int group_size = src.NumCols() / this->NumCols(),
|
|
num_rows = this->NumRows(), num_cols = this->NumCols();
|
|
for (MatrixIndexT i = 0; i < num_rows; i++)
|
|
for (MatrixIndexT j = 0; j < num_cols; j++)
|
|
(*this)(i, j) = src.Row(i).Range(j * group_size, group_size).Norm(power);
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::GroupMax(const MatrixBase<Real> &src) {
|
|
KALDI_ASSERT(src.NumCols() % this->NumCols() == 0 &&
|
|
src.NumRows() == this->NumRows());
|
|
int group_size = src.NumCols() / this->NumCols(),
|
|
num_rows = this->NumRows(), num_cols = this->NumCols();
|
|
for (MatrixIndexT i = 0; i < num_rows; i++) {
|
|
const Real *src_row_data = src.RowData(i);
|
|
for (MatrixIndexT j = 0; j < num_cols; j++) {
|
|
Real max_val = -1e20;
|
|
for (MatrixIndexT k = 0; k < group_size; k++) {
|
|
Real src_data = src_row_data[j * group_size + k];
|
|
if (src_data > max_val)
|
|
max_val = src_data;
|
|
}
|
|
(*this)(i, j) = max_val;
|
|
}
|
|
}
|
|
}
|
|
*/
|
|
template<typename Real>
|
|
void MatrixBase<Real>::CopyCols(const MatrixBase<Real> &src,
|
|
const MatrixIndexT *indices) {
|
|
KALDI_ASSERT(NumRows() == src.NumRows());
|
|
MatrixIndexT num_rows = num_rows_, num_cols = num_cols_,
|
|
this_stride = stride_, src_stride = src.stride_;
|
|
Real *this_data = this->data_;
|
|
const Real *src_data = src.data_;
|
|
#ifdef KALDI_PARANOID
|
|
MatrixIndexT src_cols = src.NumCols();
|
|
for (MatrixIndexT i = 0; i < num_cols; i++)
|
|
KALDI_ASSERT(indices[i] >= -1 && indices[i] < src_cols);
|
|
#endif
|
|
|
|
// For the sake of memory locality we do this row by row, rather
|
|
// than doing it column-wise using cublas_Xcopy
|
|
for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride, src_data += src_stride) {
|
|
const MatrixIndexT *index_ptr = &(indices[0]);
|
|
for (MatrixIndexT c = 0; c < num_cols; c++, index_ptr++) {
|
|
if (*index_ptr < 0) this_data[c] = 0;
|
|
else this_data[c] = src_data[*index_ptr];
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
template<typename Real>
|
|
void MatrixBase<Real>::AddCols(const MatrixBase<Real> &src,
|
|
const MatrixIndexT *indices) {
|
|
KALDI_ASSERT(NumRows() == src.NumRows());
|
|
MatrixIndexT num_rows = num_rows_, num_cols = num_cols_,
|
|
this_stride = stride_, src_stride = src.stride_;
|
|
Real *this_data = this->data_;
|
|
const Real *src_data = src.data_;
|
|
#ifdef KALDI_PARANOID
|
|
MatrixIndexT src_cols = src.NumCols();
|
|
for (MatrixIndexT i = 0; i < num_cols; i++)
|
|
KALDI_ASSERT(indices[i] >= -1 && indices[i] < src_cols);
|
|
#endif
|
|
|
|
// For the sake of memory locality we do this row by row, rather
|
|
// than doing it column-wise using cublas_Xcopy
|
|
for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride, src_data += src_stride) {
|
|
const MatrixIndexT *index_ptr = &(indices[0]);
|
|
for (MatrixIndexT c = 0; c < num_cols; c++, index_ptr++) {
|
|
if (*index_ptr >= 0)
|
|
this_data[c] += src_data[*index_ptr];
|
|
}
|
|
}
|
|
}*/
|
|
|
|
/*
|
|
template<typename Real>
|
|
void MatrixBase<Real>::CopyRows(const MatrixBase<Real> &src,
|
|
const MatrixIndexT *indices) {
|
|
KALDI_ASSERT(NumCols() == src.NumCols());
|
|
MatrixIndexT num_rows = num_rows_, num_cols = num_cols_,
|
|
this_stride = stride_;
|
|
Real *this_data = this->data_;
|
|
|
|
for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride) {
|
|
MatrixIndexT index = indices[r];
|
|
if (index < 0) memset(this_data, 0, sizeof(Real) * num_cols_);
|
|
else cblas_Xcopy(num_cols, src.RowData(index), 1, this_data, 1);
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::CopyRows(const Real *const *src) {
|
|
MatrixIndexT num_rows = num_rows_,
|
|
num_cols = num_cols_, this_stride = stride_;
|
|
Real *this_data = this->data_;
|
|
|
|
for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride) {
|
|
const Real *const src_data = src[r];
|
|
if (src_data == NULL) memset(this_data, 0, sizeof(Real) * num_cols);
|
|
else cblas_Xcopy(num_cols, src_data, 1, this_data, 1);
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::CopyToRows(Real *const *dst) const {
|
|
MatrixIndexT num_rows = num_rows_,
|
|
num_cols = num_cols_, this_stride = stride_;
|
|
const Real *this_data = this->data_;
|
|
|
|
for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride) {
|
|
Real *const dst_data = dst[r];
|
|
if (dst_data != NULL)
|
|
cblas_Xcopy(num_cols, this_data, 1, dst_data, 1);
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::AddRows(Real alpha,
|
|
const MatrixBase<Real> &src,
|
|
const MatrixIndexT *indexes) {
|
|
KALDI_ASSERT(NumCols() == src.NumCols());
|
|
MatrixIndexT num_rows = num_rows_,
|
|
num_cols = num_cols_, this_stride = stride_;
|
|
Real *this_data = this->data_;
|
|
|
|
for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride) {
|
|
MatrixIndexT index = indexes[r];
|
|
KALDI_ASSERT(index >= -1 && index < src.NumRows());
|
|
if (index != -1)
|
|
cblas_Xaxpy(num_cols, alpha, src.RowData(index), 1, this_data, 1);
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::AddRows(Real alpha, const Real *const *src) {
|
|
MatrixIndexT num_rows = num_rows_,
|
|
num_cols = num_cols_, this_stride = stride_;
|
|
Real *this_data = this->data_;
|
|
|
|
for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride) {
|
|
const Real *const src_data = src[r];
|
|
if (src_data != NULL)
|
|
cblas_Xaxpy(num_cols, alpha, src_data, 1, this_data, 1);
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::AddToRows(Real alpha,
|
|
const MatrixIndexT *indexes,
|
|
MatrixBase<Real> *dst) const {
|
|
KALDI_ASSERT(NumCols() == dst->NumCols());
|
|
MatrixIndexT num_rows = num_rows_,
|
|
num_cols = num_cols_, this_stride = stride_;
|
|
Real *this_data = this->data_;
|
|
|
|
for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride) {
|
|
MatrixIndexT index = indexes[r];
|
|
KALDI_ASSERT(index >= -1 && index < dst->NumRows());
|
|
if (index != -1)
|
|
cblas_Xaxpy(num_cols, alpha, this_data, 1, dst->RowData(index), 1);
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::AddToRows(Real alpha, Real *const *dst) const {
|
|
MatrixIndexT num_rows = num_rows_,
|
|
num_cols = num_cols_, this_stride = stride_;
|
|
const Real *this_data = this->data_;
|
|
|
|
for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride) {
|
|
Real *const dst_data = dst[r];
|
|
if (dst_data != NULL)
|
|
cblas_Xaxpy(num_cols, alpha, this_data, 1, dst_data, 1);
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::Sigmoid(const MatrixBase<Real> &src) {
|
|
KALDI_ASSERT(SameDim(*this, src));
|
|
|
|
if (num_cols_ == stride_ && src.num_cols_ == src.stride_) {
|
|
SubVector<Real> src_vec(src.data_, num_rows_ * num_cols_),
|
|
dst_vec(this->data_, num_rows_ * num_cols_);
|
|
dst_vec.Sigmoid(src_vec);
|
|
} else {
|
|
for (MatrixIndexT r = 0; r < num_rows_; r++) {
|
|
SubVector<Real> src_vec(src, r), dest_vec(*this, r);
|
|
dest_vec.Sigmoid(src_vec);
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::DiffSigmoid(const MatrixBase<Real> &value,
|
|
const MatrixBase<Real> &diff) {
|
|
KALDI_ASSERT(SameDim(*this, value) && SameDim(*this, diff));
|
|
MatrixIndexT num_rows = num_rows_, num_cols = num_cols_,
|
|
stride = stride_, value_stride = value.stride_, diff_stride = diff.stride_;
|
|
Real *data = data_;
|
|
const Real *value_data = value.data_, *diff_data = diff.data_;
|
|
for (MatrixIndexT r = 0; r < num_rows; r++) {
|
|
for (MatrixIndexT c = 0; c < num_cols; c++)
|
|
data[c] = diff_data[c] * value_data[c] * (1.0 - value_data[c]);
|
|
data += stride;
|
|
value_data += value_stride;
|
|
diff_data += diff_stride;
|
|
}
|
|
}
|
|
|
|
template<typename Real>
|
|
void MatrixBase<Real>::DiffTanh(const MatrixBase<Real> &value,
|
|
const MatrixBase<Real> &diff) {
|
|
KALDI_ASSERT(SameDim(*this, value) && SameDim(*this, diff));
|
|
MatrixIndexT num_rows = num_rows_, num_cols = num_cols_,
|
|
stride = stride_, value_stride = value.stride_, diff_stride = diff.stride_;
|
|
Real *data = data_;
|
|
const Real *value_data = value.data_, *diff_data = diff.data_;
|
|
for (MatrixIndexT r = 0; r < num_rows; r++) {
|
|
for (MatrixIndexT c = 0; c < num_cols; c++)
|
|
data[c] = diff_data[c] * (1.0 - (value_data[c] * value_data[c]));
|
|
data += stride;
|
|
value_data += value_stride;
|
|
diff_data += diff_stride;
|
|
}
|
|
}*/
|
|
|
|
/*
|
|
template<typename Real>
|
|
template<typename OtherReal>
|
|
void MatrixBase<Real>::AddVecToRows(const Real alpha, const VectorBase<OtherReal> &v) {
|
|
const MatrixIndexT num_rows = num_rows_, num_cols = num_cols_,
|
|
stride = stride_;
|
|
KALDI_ASSERT(v.Dim() == num_cols);
|
|
if(num_cols <= 64) {
|
|
Real *data = data_;
|
|
const OtherReal *vdata = v.Data();
|
|
for (MatrixIndexT i = 0; i < num_rows; i++, data += stride) {
|
|
for (MatrixIndexT j = 0; j < num_cols; j++)
|
|
data[j] += alpha * vdata[j];
|
|
}
|
|
|
|
} else {
|
|
Vector<OtherReal> ones(num_rows);
|
|
ones.Set(1.0);
|
|
this->AddVecVec(alpha, ones, v);
|
|
}
|
|
}
|
|
|
|
template void MatrixBase<float>::AddVecToRows(const float alpha,
|
|
const VectorBase<float> &v);
|
|
template void MatrixBase<float>::AddVecToRows(const float alpha,
|
|
const VectorBase<double> &v);
|
|
template void MatrixBase<double>::AddVecToRows(const double alpha,
|
|
const VectorBase<float> &v);
|
|
template void MatrixBase<double>::AddVecToRows(const double alpha,
|
|
const VectorBase<double> &v);
|
|
|
|
|
|
template<typename Real>
|
|
template<typename OtherReal>
|
|
void MatrixBase<Real>::AddVecToCols(const Real alpha, const VectorBase<OtherReal> &v) {
|
|
const MatrixIndexT num_rows = num_rows_, num_cols = num_cols_,
|
|
stride = stride_;
|
|
KALDI_ASSERT(v.Dim() == num_rows);
|
|
|
|
if (num_rows <= 64) {
|
|
Real *data = data_;
|
|
const OtherReal *vdata = v.Data();
|
|
for (MatrixIndexT i = 0; i < num_rows; i++, data += stride) {
|
|
Real to_add = alpha * vdata[i];
|
|
for (MatrixIndexT j = 0; j < num_cols; j++)
|
|
data[j] += to_add;
|
|
}
|
|
|
|
} else {
|
|
Vector<OtherReal> ones(num_cols);
|
|
ones.Set(1.0);
|
|
this->AddVecVec(alpha, v, ones);
|
|
}
|
|
}
|
|
|
|
template void MatrixBase<float>::AddVecToCols(const float alpha,
|
|
const VectorBase<float> &v);
|
|
template void MatrixBase<float>::AddVecToCols(const float alpha,
|
|
const VectorBase<double> &v);
|
|
template void MatrixBase<double>::AddVecToCols(const double alpha,
|
|
const VectorBase<float> &v);
|
|
template void MatrixBase<double>::AddVecToCols(const double alpha,
|
|
const VectorBase<double> &v);
|
|
*/
|
|
//Explicit instantiation of the classes
|
|
//Apparently, it seems to be necessary that the instantiation
|
|
//happens at the end of the file. Otherwise, not all the member
|
|
//functions will get instantiated.
|
|
|
|
template class Matrix<float>;
|
|
template class Matrix<double>;
|
|
template class MatrixBase<float>;
|
|
template class MatrixBase<double>;
|
|
template class SubMatrix<float>;
|
|
template class SubMatrix<double>;
|
|
|
|
} // namespace kaldi
|