kaldi-matrix.cc
Go to the documentation of this file.
1 // matrix/kaldi-matrix.cc
2 
3 // Copyright 2009-2011 Lukas Burget; Ondrej Glembek; Go Vivace Inc.;
4 // Microsoft Corporation; Saarland University;
5 // Yanmin Qian; Petr Schwarz; Jan Silovsky;
6 // Haihua Xu
7 // 2017 Shiyin Kang
8 // 2019 Yiwen Shao
9 
10 // See ../../COPYING for clarification regarding multiple authors
11 //
12 // Licensed under the Apache License, Version 2.0 (the "License");
13 // you may not use this file except in compliance with the License.
14 // You may obtain a copy of the License at
15 //
16 // http://www.apache.org/licenses/LICENSE-2.0
17 //
18 // THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
19 // KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
20 // WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
21 // MERCHANTABLITY OR NON-INFRINGEMENT.
22 // See the Apache 2 License for the specific language governing permissions and
23 // limitations under the License.
24 
25 #include "matrix/kaldi-matrix.h"
26 #include "matrix/sp-matrix.h"
27 #include "matrix/jama-svd.h"
28 #include "matrix/jama-eig.h"
30 #include "matrix/sparse-matrix.h"
31 
32 static_assert(int(kaldi::kNoTrans) == int(CblasNoTrans) && int(kaldi::kTrans) == int(CblasTrans),
33  "kaldi::kNoTrans and kaldi::kTrans must be equal to the appropriate CBLAS library constants!");
34 
35 namespace kaldi {
36 
37 template<typename Real>
38 void MatrixBase<Real>::Invert(Real *log_det, Real *det_sign,
39  bool inverse_needed) {
40  KALDI_ASSERT(num_rows_ == num_cols_);
41  if (num_rows_ == 0) {
42  if (det_sign) *det_sign = 1;
43  if (log_det) *log_det = 0.0;
44  return;
45  }
46 #ifndef HAVE_ATLAS
47  KaldiBlasInt *pivot = new KaldiBlasInt[num_rows_];
48  KaldiBlasInt M = num_rows_;
49  KaldiBlasInt N = num_cols_;
50  KaldiBlasInt LDA = stride_;
51  KaldiBlasInt result = -1;
52  KaldiBlasInt l_work = std::max<KaldiBlasInt>(1, N);
53  Real *p_work;
54  void *temp;
55  if ((p_work = static_cast<Real*>(
56  KALDI_MEMALIGN(16, sizeof(Real)*l_work, &temp))) == NULL) {
57  delete[] pivot;
58  throw std::bad_alloc();
59  }
60 
61  clapack_Xgetrf2(&M, &N, data_, &LDA, pivot, &result);
62  const int pivot_offset = 1;
63 #else
64  int *pivot = new int[num_rows_];
65  int result;
66  clapack_Xgetrf(num_rows_, num_cols_, data_, stride_, pivot, &result);
67  const int pivot_offset = 0;
68 #endif
69  KALDI_ASSERT(result >= 0 && "Call to CLAPACK sgetrf_ or ATLAS clapack_sgetrf "
70  "called with wrong arguments");
71  if (result > 0) {
72  if (inverse_needed) {
73  KALDI_ERR << "Cannot invert: matrix is singular";
74  } else {
75  if (log_det) *log_det = -std::numeric_limits<Real>::infinity();
76  if (det_sign) *det_sign = 0;
77  delete[] pivot;
78 #ifndef HAVE_ATLAS
79  KALDI_MEMALIGN_FREE(p_work);
80 #endif
81  return;
82  }
83  }
84  if (det_sign != NULL) {
85  int sign = 1;
86  for (MatrixIndexT i = 0; i < num_rows_; i++)
87  if (pivot[i] != static_cast<int>(i) + pivot_offset) sign *= -1;
88  *det_sign = sign;
89  }
90  if (log_det != NULL || det_sign != NULL) { // Compute log determinant.
91  if (log_det != NULL) *log_det = 0.0;
92  Real prod = 1.0;
93  for (MatrixIndexT i = 0; i < num_rows_; i++) {
94  prod *= (*this)(i, i);
95  if (i == num_rows_ - 1 || std::fabs(prod) < 1.0e-10 ||
96  std::fabs(prod) > 1.0e+10) {
97  if (log_det != NULL) *log_det += kaldi::Log(std::fabs(prod));
98  if (det_sign != NULL) *det_sign *= (prod > 0 ? 1.0 : -1.0);
99  prod = 1.0;
100  }
101  }
102  }
103 #ifndef HAVE_ATLAS
104  if (inverse_needed) clapack_Xgetri2(&M, data_, &LDA, pivot, p_work, &l_work,
105  &result);
106  delete[] pivot;
107  KALDI_MEMALIGN_FREE(p_work);
108 #else
109  if (inverse_needed)
110  clapack_Xgetri(num_rows_, data_, stride_, pivot, &result);
111  delete [] pivot;
112 #endif
113  KALDI_ASSERT(result == 0 && "Call to CLAPACK sgetri_ or ATLAS clapack_sgetri "
114  "called with wrong arguments");
115 }
116 
117 template<>
118 template<>
119 void MatrixBase<float>::AddVecVec(const float alpha,
120  const VectorBase<float> &a,
121  const VectorBase<float> &rb) {
122  KALDI_ASSERT(a.Dim() == num_rows_ && rb.Dim() == num_cols_);
123  cblas_Xger(a.Dim(), rb.Dim(), alpha, a.Data(), 1, rb.Data(),
124  1, data_, stride_);
125 }
126 
127 template<typename Real>
128 template<typename OtherReal>
129 void MatrixBase<Real>::AddVecVec(const Real alpha,
130  const VectorBase<OtherReal> &a,
131  const VectorBase<OtherReal> &b) {
132  KALDI_ASSERT(a.Dim() == num_rows_ && b.Dim() == num_cols_);
133  if (num_rows_ * num_cols_ > 100) { // It's probably worth it to allocate
134  // temporary vectors of the right type and use BLAS.
135  Vector<Real> temp_a(a), temp_b(b);
136  cblas_Xger(num_rows_, num_cols_, alpha, temp_a.Data(), 1,
137  temp_b.Data(), 1, data_, stride_);
138  } else {
139  const OtherReal *a_data = a.Data(), *b_data = b.Data();
140  Real *row_data = data_;
141  for (MatrixIndexT i = 0; i < num_rows_; i++, row_data += stride_) {
142  BaseFloat alpha_ai = alpha * a_data[i];
143  for (MatrixIndexT j = 0; j < num_cols_; j++)
144  row_data[j] += alpha_ai * b_data[j];
145  }
146  }
147 }
148 
149 // instantiate the template above.
150 template
151 void MatrixBase<float>::AddVecVec(const float alpha,
152  const VectorBase<double> &a,
153  const VectorBase<double> &b);
154 template
155 void MatrixBase<double>::AddVecVec(const double alpha,
156  const VectorBase<float> &a,
157  const VectorBase<float> &b);
158 
159 template<>
160 template<>
161 void MatrixBase<double>::AddVecVec(const double alpha,
162  const VectorBase<double> &a,
163  const VectorBase<double> &rb) {
164  KALDI_ASSERT(a.Dim() == num_rows_ && rb.Dim() == num_cols_);
165  if (num_rows_ == 0) return;
166  cblas_Xger(a.Dim(), rb.Dim(), alpha, a.Data(), 1, rb.Data(),
167  1, data_, stride_);
168 }
169 
170 template<typename Real>
171 void MatrixBase<Real>::AddMatMat(const Real alpha,
172  const MatrixBase<Real>& A,
173  MatrixTransposeType transA,
174  const MatrixBase<Real>& B,
175  MatrixTransposeType transB,
176  const Real beta) {
177  KALDI_ASSERT((transA == kNoTrans && transB == kNoTrans && A.num_cols_ == B.num_rows_ && A.num_rows_ == num_rows_ && B.num_cols_ == num_cols_)
178  || (transA == kTrans && transB == kNoTrans && A.num_rows_ == B.num_rows_ && A.num_cols_ == num_rows_ && B.num_cols_ == num_cols_)
179  || (transA == kNoTrans && transB == kTrans && A.num_cols_ == B.num_cols_ && A.num_rows_ == num_rows_ && B.num_rows_ == num_cols_)
180  || (transA == kTrans && transB == kTrans && A.num_rows_ == B.num_cols_ && A.num_cols_ == num_rows_ && B.num_rows_ == num_cols_));
181  KALDI_ASSERT(&A != this && &B != this);
182  if (num_rows_ == 0) return;
183  cblas_Xgemm(alpha, transA, A.data_, A.num_rows_, A.num_cols_, A.stride_,
184  transB, B.data_, B.stride_, beta, data_, num_rows_, num_cols_, stride_);
185 
186 }
187 
188 template<typename Real>
190  const MatrixBase<Real>& B,
191  const MatrixBase<Real>& C) {
192  KALDI_ASSERT(A.NumRows() == B.NumRows() && A.NumCols() == B.NumCols());
193  KALDI_ASSERT(A.NumRows() == C.NumRows() && A.NumCols() == C.NumCols());
194  for (int32 r = 0; r < A.NumRows(); r++) { // each frame...
195  for (int32 c = 0; c < A.NumCols(); c++) {
196  BaseFloat i = C(r, c), o = B(r, c), od = A(r, c),
197  id;
198  if (i != 0.0) {
199  id = od * (o / i);
200  } else {
201  id = od;
202  }
205  (*this)(r, c) = id;
206  }
207  }
208 }
209 
210 
211 template<typename Real>
213  KALDI_ASSERT(num_rows_ == num_cols_);
214  Real *data = data_;
215  MatrixIndexT num_rows = num_rows_, stride = stride_;
216  for (int32 i = 0; i < num_rows; i++)
217  for (int32 j = 0; j < i; j++)
218  data[j * stride + i ] = data[i * stride + j];
219 }
220 
221 
222 template<typename Real>
224  KALDI_ASSERT(num_rows_ == num_cols_);
225  Real *data = data_;
226  MatrixIndexT num_rows = num_rows_, stride = stride_;
227  for (int32 i = 0; i < num_rows; i++)
228  for (int32 j = 0; j < i; j++)
229  data[i * stride + j] = data[j * stride + i];
230 }
231 
232 template<typename Real>
233 void MatrixBase<Real>::SymAddMat2(const Real alpha,
234  const MatrixBase<Real> &A,
235  MatrixTransposeType transA,
236  Real beta) {
237  KALDI_ASSERT(num_rows_ == num_cols_ &&
238  ((transA == kNoTrans && A.num_rows_ == num_rows_) ||
239  (transA == kTrans && A.num_cols_ == num_cols_)));
240  KALDI_ASSERT(A.data_ != data_);
241  if (num_rows_ == 0) return;
242 
251 #ifdef HAVE_ATLAS
252  if (transA == kTrans && num_rows_ >= 56) {
253  this->AddMatMat(alpha, A, kTrans, A, kNoTrans, beta);
254  return;
255  }
256 #endif // HAVE_ATLAS
257 
258  MatrixIndexT A_other_dim = (transA == kNoTrans ? A.num_cols_ : A.num_rows_);
259 
260  // This function call is hard-coded to update the lower triangle.
261  cblas_Xsyrk(transA, num_rows_, A_other_dim, alpha, A.Data(),
262  A.Stride(), beta, this->data_, this->stride_);
263 }
264 
265 
266 template<typename Real>
267 void MatrixBase<Real>::AddMatSmat(const Real alpha,
268  const MatrixBase<Real> &A,
269  MatrixTransposeType transA,
270  const MatrixBase<Real> &B,
271  MatrixTransposeType transB,
272  const Real beta) {
273  KALDI_ASSERT((transA == kNoTrans && transB == kNoTrans && A.num_cols_ == B.num_rows_ && A.num_rows_ == num_rows_ && B.num_cols_ == num_cols_)
274  || (transA == kTrans && transB == kNoTrans && A.num_rows_ == B.num_rows_ && A.num_cols_ == num_rows_ && B.num_cols_ == num_cols_)
275  || (transA == kNoTrans && transB == kTrans && A.num_cols_ == B.num_cols_ && A.num_rows_ == num_rows_ && B.num_rows_ == num_cols_)
276  || (transA == kTrans && transB == kTrans && A.num_rows_ == B.num_cols_ && A.num_cols_ == num_rows_ && B.num_rows_ == num_cols_));
277  KALDI_ASSERT(&A != this && &B != this);
278 
279  // We iterate over the columns of B.
280 
281  MatrixIndexT Astride = A.stride_, Bstride = B.stride_, stride = this->stride_,
282  Arows = A.num_rows_, Acols = A.num_cols_;
283  Real *data = this->data_, *Adata = A.data_, *Bdata = B.data_;
284  MatrixIndexT num_cols = this->num_cols_;
285  if (transB == kNoTrans) {
286  // Iterate over the columns of *this and of B.
287  for (MatrixIndexT c = 0; c < num_cols; c++) {
288  // for each column of *this, do
289  // [this column] = [alpha * A * this column of B] + [beta * this column]
290  Xgemv_sparsevec(transA, Arows, Acols, alpha, Adata, Astride,
291  Bdata + c, Bstride, beta, data + c, stride);
292  }
293  } else {
294  // Iterate over the columns of *this and the rows of B.
295  for (MatrixIndexT c = 0; c < num_cols; c++) {
296  // for each column of *this, do
297  // [this column] = [alpha * A * this row of B] + [beta * this column]
298  Xgemv_sparsevec(transA, Arows, Acols, alpha, Adata, Astride,
299  Bdata + (c * Bstride), 1, beta, data + c, stride);
300  }
301  }
302 }
303 
304 template<typename Real>
305 void MatrixBase<Real>::AddSmatMat(const Real alpha,
306  const MatrixBase<Real> &A,
307  MatrixTransposeType transA,
308  const MatrixBase<Real> &B,
309  MatrixTransposeType transB,
310  const Real beta) {
311  KALDI_ASSERT((transA == kNoTrans && transB == kNoTrans && A.num_cols_ == B.num_rows_ && A.num_rows_ == num_rows_ && B.num_cols_ == num_cols_)
312  || (transA == kTrans && transB == kNoTrans && A.num_rows_ == B.num_rows_ && A.num_cols_ == num_rows_ && B.num_cols_ == num_cols_)
313  || (transA == kNoTrans && transB == kTrans && A.num_cols_ == B.num_cols_ && A.num_rows_ == num_rows_ && B.num_rows_ == num_cols_)
314  || (transA == kTrans && transB == kTrans && A.num_rows_ == B.num_cols_ && A.num_cols_ == num_rows_ && B.num_rows_ == num_cols_));
315  KALDI_ASSERT(&A != this && &B != this);
316 
317  MatrixIndexT Astride = A.stride_, Bstride = B.stride_, stride = this->stride_,
318  Brows = B.num_rows_, Bcols = B.num_cols_;
319  MatrixTransposeType invTransB = (transB == kTrans ? kNoTrans : kTrans);
320  Real *data = this->data_, *Adata = A.data_, *Bdata = B.data_;
321  MatrixIndexT num_rows = this->num_rows_;
322  if (transA == kNoTrans) {
323  // Iterate over the rows of *this and of A.
324  for (MatrixIndexT r = 0; r < num_rows; r++) {
325  // for each row of *this, do
326  // [this row] = [alpha * (this row of A) * B^T] + [beta * this row]
327  Xgemv_sparsevec(invTransB, Brows, Bcols, alpha, Bdata, Bstride,
328  Adata + (r * Astride), 1, beta, data + (r * stride), 1);
329  }
330  } else {
331  // Iterate over the rows of *this and the columns of A.
332  for (MatrixIndexT r = 0; r < num_rows; r++) {
333  // for each row of *this, do
334  // [this row] = [alpha * (this column of A) * B^T] + [beta * this row]
335  Xgemv_sparsevec(invTransB, Brows, Bcols, alpha, Bdata, Bstride,
336  Adata + r, Astride, beta, data + (r * stride), 1);
337  }
338  }
339 }
340 
341 template<typename Real>
342 void MatrixBase<Real>::AddSpSp(const Real alpha, const SpMatrix<Real> &A_in,
343  const SpMatrix<Real> &B_in, const Real beta) {
344  MatrixIndexT sz = num_rows_;
345  KALDI_ASSERT(sz == num_cols_ && sz == A_in.NumRows() && sz == B_in.NumRows());
346 
347  Matrix<Real> A(A_in), B(B_in);
348  // CblasLower or CblasUpper would work below as symmetric matrix is copied
349  // fully (to save work, we used the matrix constructor from SpMatrix).
350  // CblasLeft means A is on the left: C <-- alpha A B + beta C
351  if (sz == 0) return;
352  cblas_Xsymm(alpha, sz, A.data_, A.stride_, B.data_, B.stride_, beta, data_, stride_);
353 }
354 
355 template<typename Real>
356 void MatrixBase<Real>::AddMat(const Real alpha, const MatrixBase<Real>& A,
357  MatrixTransposeType transA) {
358  if (&A == this) {
359  if (transA == kNoTrans) {
360  Scale(alpha + 1.0);
361  } else {
362  KALDI_ASSERT(num_rows_ == num_cols_ && "AddMat: adding to self (transposed): not symmetric.");
363  Real *data = data_;
364  if (alpha == 1.0) { // common case-- handle separately.
365  for (MatrixIndexT row = 0; row < num_rows_; row++) {
366  for (MatrixIndexT col = 0; col < row; col++) {
367  Real *lower = data + (row * stride_) + col, *upper = data + (col
368  * stride_) + row;
369  Real sum = *lower + *upper;
370  *lower = *upper = sum;
371  }
372  *(data + (row * stride_) + row) *= 2.0; // diagonal.
373  }
374  } else {
375  for (MatrixIndexT row = 0; row < num_rows_; row++) {
376  for (MatrixIndexT col = 0; col < row; col++) {
377  Real *lower = data + (row * stride_) + col, *upper = data + (col
378  * stride_) + row;
379  Real lower_tmp = *lower;
380  *lower += alpha * *upper;
381  *upper += alpha * lower_tmp;
382  }
383  *(data + (row * stride_) + row) *= (1.0 + alpha); // diagonal.
384  }
385  }
386  }
387  } else {
388  int aStride = (int) A.stride_, stride = stride_;
389  Real *adata = A.data_, *data = data_;
390  if (transA == kNoTrans) {
391  KALDI_ASSERT(A.num_rows_ == num_rows_ && A.num_cols_ == num_cols_);
392  if (num_rows_ == 0) return;
393  for (MatrixIndexT row = 0; row < num_rows_; row++, adata += aStride,
394  data += stride) {
395  cblas_Xaxpy(num_cols_, alpha, adata, 1, data, 1);
396  }
397  } else {
398  KALDI_ASSERT(A.num_cols_ == num_rows_ && A.num_rows_ == num_cols_);
399  if (num_rows_ == 0) return;
400  for (MatrixIndexT row = 0; row < num_rows_; row++, adata++, data += stride)
401  cblas_Xaxpy(num_cols_, alpha, adata, aStride, data, 1);
402  }
403  }
404 }
405 
406 template<typename Real>
408  MatrixTransposeType trans) {
409  if (trans == kNoTrans) {
410  KALDI_ASSERT(NumRows() == A.NumRows());
411  KALDI_ASSERT(NumCols() == A.NumCols());
412  MatrixIndexT a_num_rows = A.NumRows();
413  for (MatrixIndexT i = 0; i < a_num_rows; ++i) {
414  const SparseVector<Real> &row = A.Row(i);
415  MatrixIndexT num_elems = row.NumElements();
416  for (MatrixIndexT id = 0; id < num_elems; ++id) {
417  (*this)(i, row.GetElement(id).first) += alpha
418  * row.GetElement(id).second;
419  }
420  }
421  } else {
422  KALDI_ASSERT(NumRows() == A.NumCols());
423  KALDI_ASSERT(NumCols() == A.NumRows());
424  MatrixIndexT a_num_rows = A.NumRows();
425  for (MatrixIndexT i = 0; i < a_num_rows; ++i) {
426  const SparseVector<Real> &row = A.Row(i);
427  MatrixIndexT num_elems = row.NumElements();
428  for (MatrixIndexT id = 0; id < num_elems; ++id) {
429  (*this)(row.GetElement(id).first, i) += alpha
430  * row.GetElement(id).second;
431  }
432  }
433  }
434 }
435 
436 template<typename Real>
438  MatrixTransposeType transA,
439  const MatrixBase<Real> &B, Real beta) {
440  if (transA == kNoTrans) {
441  KALDI_ASSERT(NumRows() == A.NumRows());
442  KALDI_ASSERT(NumCols() == B.NumCols());
443  KALDI_ASSERT(A.NumCols() == B.NumRows());
444 
445  this->Scale(beta);
446  MatrixIndexT a_num_rows = A.NumRows(),
447  this_num_cols = this->NumCols();
448  for (MatrixIndexT i = 0; i < a_num_rows; ++i) {
449  Real *this_row_i = this->RowData(i);
450  const SparseVector<Real> &A_row_i = A.Row(i);
451  MatrixIndexT num_elems = A_row_i.NumElements();
452  for (MatrixIndexT e = 0; e < num_elems; ++e) {
453  const std::pair<MatrixIndexT, Real> &p = A_row_i.GetElement(e);
454  MatrixIndexT k = p.first;
455  Real alpha_A_ik = alpha * p.second;
456  const Real *b_row_k = B.RowData(k);
457  cblas_Xaxpy(this_num_cols, alpha_A_ik, b_row_k, 1,
458  this_row_i, 1);
459  //for (MatrixIndexT j = 0; j < this_num_cols; ++j)
460  // this_row_i[j] += alpha_A_ik * b_row_k[j];
461  }
462  }
463  } else {
464  KALDI_ASSERT(NumRows() == A.NumCols());
465  KALDI_ASSERT(NumCols() == B.NumCols());
466  KALDI_ASSERT(A.NumRows() == B.NumRows());
467 
468  this->Scale(beta);
469  Matrix<Real> buf(NumRows(), NumCols(), kSetZero);
470  MatrixIndexT a_num_rows = A.NumRows(),
471  this_num_cols = this->NumCols();
472  for (int k = 0; k < a_num_rows; ++k) {
473  const Real *b_row_k = B.RowData(k);
474  const SparseVector<Real> &A_row_k = A.Row(k);
475  MatrixIndexT num_elems = A_row_k.NumElements();
476  for (MatrixIndexT e = 0; e < num_elems; ++e) {
477  const std::pair<MatrixIndexT, Real> &p = A_row_k.GetElement(e);
478  MatrixIndexT i = p.first;
479  Real alpha_A_ki = alpha * p.second;
480  Real *this_row_i = this->RowData(i);
481  cblas_Xaxpy(this_num_cols, alpha_A_ki, b_row_k, 1,
482  this_row_i, 1);
483  //for (MatrixIndexT j = 0; j < this_num_cols; ++j)
484  // this_row_i[j] += alpha_A_ki * b_row_k[j];
485  }
486  }
487  }
488 }
489 
490 template<typename Real>
492  const SparseMatrix<Real> &B,
493  MatrixTransposeType transB, Real beta) {
494  if (transB == kNoTrans) {
495  KALDI_ASSERT(NumRows() == A.NumRows());
496  KALDI_ASSERT(NumCols() == B.NumCols());
497  KALDI_ASSERT(A.NumCols() == B.NumRows());
498 
499  this->Scale(beta);
500  MatrixIndexT b_num_rows = B.NumRows(),
501  this_num_rows = this->NumRows();
502  // Iterate over the rows of sparse matrix B and columns of A.
503  for (MatrixIndexT k = 0; k < b_num_rows; ++k) {
504  const SparseVector<Real> &B_row_k = B.Row(k);
505  MatrixIndexT num_elems = B_row_k.NumElements();
506  const Real *a_col_k = A.Data() + k;
507  for (MatrixIndexT e = 0; e < num_elems; ++e) {
508  const std::pair<MatrixIndexT, Real> &p = B_row_k.GetElement(e);
509  MatrixIndexT j = p.first;
510  Real alpha_B_kj = alpha * p.second;
511  Real *this_col_j = this->Data() + j;
512  // Add to entire 'j'th column of *this at once using cblas_Xaxpy.
513  // pass stride to write a colmun as matrices are stored in row major order.
514  cblas_Xaxpy(this_num_rows, alpha_B_kj, a_col_k, A.stride_,
515  this_col_j, this->stride_);
516  //for (MatrixIndexT i = 0; i < this_num_rows; ++i)
517  // this_col_j[i*this->stride_] += alpha_B_kj * a_col_k[i*A.stride_];
518  }
519  }
520  } else {
521  KALDI_ASSERT(NumRows() == A.NumRows());
522  KALDI_ASSERT(NumCols() == B.NumRows());
523  KALDI_ASSERT(A.NumCols() == B.NumCols());
524 
525  this->Scale(beta);
526  MatrixIndexT b_num_rows = B.NumRows(),
527  this_num_rows = this->NumRows();
528  // Iterate over the rows of sparse matrix B and columns of *this.
529  for (MatrixIndexT j = 0; j < b_num_rows; ++j) {
530  const SparseVector<Real> &B_row_j = B.Row(j);
531  MatrixIndexT num_elems = B_row_j.NumElements();
532  Real *this_col_j = this->Data() + j;
533  for (MatrixIndexT e = 0; e < num_elems; ++e) {
534  const std::pair<MatrixIndexT, Real> &p = B_row_j.GetElement(e);
535  MatrixIndexT k = p.first;
536  Real alpha_B_jk = alpha * p.second;
537  const Real *a_col_k = A.Data() + k;
538  // Add to entire 'j'th column of *this at once using cblas_Xaxpy.
539  // pass stride to write a column as matrices are stored in row major order.
540  cblas_Xaxpy(this_num_rows, alpha_B_jk, a_col_k, A.stride_,
541  this_col_j, this->stride_);
542  //for (MatrixIndexT i = 0; i < this_num_rows; ++i)
543  // this_col_j[i*this->stride_] += alpha_B_jk * a_col_k[i*A.stride_];
544  }
545  }
546  }
547 }
548 
549 template<typename Real>
550 template<typename OtherReal>
551 void MatrixBase<Real>::AddSp(const Real alpha, const SpMatrix<OtherReal> &S) {
552  KALDI_ASSERT(S.NumRows() == NumRows() && S.NumRows() == NumCols());
553  Real *data = data_; const OtherReal *sdata = S.Data();
554  MatrixIndexT num_rows = NumRows(), stride = Stride();
555  for (MatrixIndexT i = 0; i < num_rows; i++) {
556  for (MatrixIndexT j = 0; j < i; j++, sdata++) {
557  data[i*stride + j] += alpha * *sdata;
558  data[j*stride + i] += alpha * *sdata;
559  }
560  data[i*stride + i] += alpha * *sdata++;
561  }
562 }
563 
564 // instantiate the template above.
565 template
566 void MatrixBase<float>::AddSp(const float alpha, const SpMatrix<float> &S);
567 template
568 void MatrixBase<double>::AddSp(const double alpha, const SpMatrix<double> &S);
569 template
570 void MatrixBase<float>::AddSp(const float alpha, const SpMatrix<double> &S);
571 template
572 void MatrixBase<double>::AddSp(const double alpha, const SpMatrix<float> &S);
573 
574 
575 template<typename Real>
577  const Real alpha, const VectorBase<Real> &v,
578  const MatrixBase<Real> &M,
579  MatrixTransposeType transM,
580  Real beta) {
581  if (beta != 1.0) this->Scale(beta);
582 
583  if (transM == kNoTrans) {
584  KALDI_ASSERT(SameDim(*this, M));
585  } else {
586  KALDI_ASSERT(M.NumRows() == NumCols() && M.NumCols() == NumRows());
587  }
588  KALDI_ASSERT(v.Dim() == this->NumRows());
589 
590  MatrixIndexT M_row_stride = M.Stride(), M_col_stride = 1, stride = stride_,
591  num_rows = num_rows_, num_cols = num_cols_;
592  if (transM == kTrans) std::swap(M_row_stride, M_col_stride);
593  Real *data = data_;
594  const Real *Mdata = M.Data(), *vdata = v.Data();
595  if (num_rows_ == 0) return;
596  for (MatrixIndexT i = 0; i < num_rows; i++, data += stride, Mdata += M_row_stride, vdata++)
597  cblas_Xaxpy(num_cols, alpha * *vdata, Mdata, M_col_stride, data, 1);
598 }
599 
600 template<typename Real>
602  const Real alpha,
603  const MatrixBase<Real> &M, MatrixTransposeType transM,
604  VectorBase<Real> &v,
605  Real beta) {
606 
607  if (beta != 1.0) this->Scale(beta);
608 
609  if (transM == kNoTrans) {
610  KALDI_ASSERT(SameDim(*this, M));
611  } else {
612  KALDI_ASSERT(M.NumRows() == NumCols() && M.NumCols() == NumRows());
613  }
614  KALDI_ASSERT(v.Dim() == this->NumCols());
615 
616  MatrixIndexT M_row_stride = M.Stride(),
617  M_col_stride = 1,
618  stride = stride_,
619  num_rows = num_rows_,
620  num_cols = num_cols_;
621 
622  if (transM == kTrans)
623  std::swap(M_row_stride, M_col_stride);
624 
625  Real *data = data_;
626  const Real *Mdata = M.Data(), *vdata = v.Data();
627  if (num_rows_ == 0) return;
628  for (MatrixIndexT i = 0; i < num_rows; i++){
629  for(MatrixIndexT j = 0; j < num_cols; j ++ ){
630  data[i*stride + j] += alpha * vdata[j] * Mdata[i*M_row_stride + j*M_col_stride];
631  }
632  }
633 }
634 
635 template<typename Real>
637  const MatrixBase<Real>& A,
638  const MatrixBase<Real>& B,
639  const Real beta) {
640  KALDI_ASSERT(A.NumRows() == B.NumRows() && A.NumCols() == B.NumCols());
641  KALDI_ASSERT(A.NumRows() == NumRows() && A.NumCols() == NumCols());
642  Real *data = data_;
643  const Real *dataA = A.Data();
644  const Real *dataB = B.Data();
645 
646  for (MatrixIndexT i = 0; i < num_rows_; i++) {
647  for (MatrixIndexT j = 0; j < num_cols_; j++) {
648  data[j] = beta*data[j] + alpha*dataA[j]*dataB[j];
649  }
650  data += Stride();
651  dataA += A.Stride();
652  dataB += B.Stride();
653  }
654 }
655 
656 #if !defined(HAVE_ATLAS) && !defined(USE_KALDI_SVD)
657 // ****************************************************************************
658 // ****************************************************************************
659 template<typename Real>
661  MatrixBase<Real> *V_in) {
662  KALDI_ASSERT(s != NULL && U_in != this && V_in != this);
663 
664  Matrix<Real> tmpU, tmpV;
665  if (U_in == NULL) tmpU.Resize(this->num_rows_, 1); // work-space if U_in empty.
666  if (V_in == NULL) tmpV.Resize(1, this->num_cols_); // work-space if V_in empty.
667 
671 
672  KaldiBlasInt M = num_cols_;
673  KaldiBlasInt N = num_rows_;
674  KaldiBlasInt LDA = Stride();
675 
676  KALDI_ASSERT(N>=M); // NumRows >= columns.
677 
678  if (U_in) {
679  KALDI_ASSERT((int)U_in->num_rows_ == N && (int)U_in->num_cols_ == M);
680  }
681  if (V_in) {
682  KALDI_ASSERT((int)V_in->num_rows_ == M && (int)V_in->num_cols_ == M);
683  }
684  KALDI_ASSERT((int)s->Dim() == std::min(M, N));
685 
686  MatrixBase<Real> *U = (U_in ? U_in : &tmpU);
687  MatrixBase<Real> *V = (V_in ? V_in : &tmpV);
688 
689  KaldiBlasInt V_stride = V->Stride();
690  KaldiBlasInt U_stride = U->Stride();
691 
692  // Original LAPACK recipe
693  // KaldiBlasInt l_work = std::max(std::max<long int>
694  // (1, 3*std::min(M, N)+std::max(M, N)), 5*std::min(M, N))*2;
695  KaldiBlasInt l_work = -1;
696  Real work_query;
697  KaldiBlasInt result;
698 
699  // query for work space
700  char *u_job = const_cast<char*>(U_in ? "s" : "N"); // "s" == skinny, "N" == "none."
701  char *v_job = const_cast<char*>(V_in ? "s" : "N"); // "s" == skinny, "N" == "none."
702  clapack_Xgesvd(v_job, u_job,
703  &M, &N, data_, &LDA,
704  s->Data(),
705  V->Data(), &V_stride,
706  U->Data(), &U_stride,
707  &work_query, &l_work,
708  &result);
709 
710  KALDI_ASSERT(result >= 0 && "Call to CLAPACK dgesvd_ called with wrong arguments");
711 
712  l_work = static_cast<KaldiBlasInt>(work_query);
713  Real *p_work;
714  void *temp;
715  if ((p_work = static_cast<Real*>(
716  KALDI_MEMALIGN(16, sizeof(Real)*l_work, &temp))) == NULL)
717  throw std::bad_alloc();
718 
719  // perform svd
720  clapack_Xgesvd(v_job, u_job,
721  &M, &N, data_, &LDA,
722  s->Data(),
723  V->Data(), &V_stride,
724  U->Data(), &U_stride,
725  p_work, &l_work,
726  &result);
727 
728  KALDI_ASSERT(result >= 0 && "Call to CLAPACK dgesvd_ called with wrong arguments");
729 
730  if (result != 0) {
731  KALDI_WARN << "CLAPACK sgesvd_ : some weird convergence not satisfied";
732  }
733  KALDI_MEMALIGN_FREE(p_work);
734 }
735 
736 #endif
737 
738 // Copy constructor. Copies data to newly allocated memory.
739 template<typename Real>
741  MatrixTransposeType trans/*=kNoTrans*/)
742  : MatrixBase<Real>() {
743  if (trans == kNoTrans) {
744  Resize(M.num_rows_, M.num_cols_);
745  this->CopyFromMat(M);
746  } else {
747  Resize(M.num_cols_, M.num_rows_);
748  this->CopyFromMat(M, kTrans);
749  }
750 }
751 
752 // Copy constructor. Copies data to newly allocated memory.
753 template<typename Real>
755  MatrixBase<Real>() {
756  Resize(M.num_rows_, M.num_cols_);
757  this->CopyFromMat(M);
758 }
759 
761 template<typename Real>
762 template<typename OtherReal>
764  MatrixTransposeType trans) : MatrixBase<Real>() {
765  if (trans == kNoTrans) {
766  Resize(M.NumRows(), M.NumCols());
767  this->CopyFromMat(M);
768  } else {
769  Resize(M.NumCols(), M.NumRows());
770  this->CopyFromMat(M, kTrans);
771  }
772 }
773 
774 // Instantiate this constructor for float->double and double->float.
775 template
777  MatrixTransposeType trans);
778 template
780  MatrixTransposeType trans);
781 
782 template<typename Real>
783 inline void Matrix<Real>::Init(const MatrixIndexT rows,
784  const MatrixIndexT cols,
785  const MatrixStrideType stride_type) {
786  if (rows * cols == 0) {
787  KALDI_ASSERT(rows == 0 && cols == 0);
788  this->num_rows_ = 0;
789  this->num_cols_ = 0;
790  this->stride_ = 0;
791  this->data_ = NULL;
792  return;
793  }
794  KALDI_ASSERT(rows > 0 && cols > 0);
795  MatrixIndexT skip, stride;
796  size_t size;
797  void *data; // aligned memory block
798  void *temp; // memory block to be really freed
799 
800  // compute the size of skip and real cols
801  skip = ((16 / sizeof(Real)) - cols % (16 / sizeof(Real)))
802  % (16 / sizeof(Real));
803  stride = cols + skip;
804  size = static_cast<size_t>(rows) * static_cast<size_t>(stride)
805  * sizeof(Real);
806 
807  // allocate the memory and set the right dimensions and parameters
808  if (NULL != (data = KALDI_MEMALIGN(16, size, &temp))) {
809  MatrixBase<Real>::data_ = static_cast<Real *> (data);
812  MatrixBase<Real>::stride_ = (stride_type == kDefaultStride ? stride : cols);
813  } else {
814  throw std::bad_alloc();
815  }
816 }
817 
818 template<typename Real>
820  const MatrixIndexT cols,
821  MatrixResizeType resize_type,
822  MatrixStrideType stride_type) {
823  // the next block uses recursion to handle what we have to do if
824  // resize_type == kCopyData.
825  if (resize_type == kCopyData) {
826  if (this->data_ == NULL || rows == 0) resize_type = kSetZero; // nothing to copy.
827  else if (rows == this->num_rows_ && cols == this->num_cols_ &&
828  (stride_type == kDefaultStride || this->stride_ == this->num_cols_)) { return; } // nothing to do.
829  else {
830  // set tmp to a matrix of the desired size; if new matrix
831  // is bigger in some dimension, zero it.
832  MatrixResizeType new_resize_type =
833  (rows > this->num_rows_ || cols > this->num_cols_) ? kSetZero : kUndefined;
834  Matrix<Real> tmp(rows, cols, new_resize_type, stride_type);
835  MatrixIndexT rows_min = std::min(rows, this->num_rows_),
836  cols_min = std::min(cols, this->num_cols_);
837  tmp.Range(0, rows_min, 0, cols_min).
838  CopyFromMat(this->Range(0, rows_min, 0, cols_min));
839  tmp.Swap(this);
840  // and now let tmp go out of scope, deleting what was in *this.
841  return;
842  }
843  }
844  // At this point, resize_type == kSetZero or kUndefined.
845 
846  if (MatrixBase<Real>::data_ != NULL) {
847  if (rows == MatrixBase<Real>::num_rows_
848  && cols == MatrixBase<Real>::num_cols_) {
849  if (resize_type == kSetZero)
850  this->SetZero();
851  return;
852  }
853  else
854  Destroy();
855  }
856  Init(rows, cols, stride_type);
857  if (resize_type == kSetZero) MatrixBase<Real>::SetZero();
858 }
859 
860 template<typename Real>
861 template<typename OtherReal>
863  MatrixTransposeType Trans) {
864  if (sizeof(Real) == sizeof(OtherReal) &&
865  static_cast<const void*>(M.Data()) ==
866  static_cast<const void*>(this->Data())) {
867  // CopyFromMat called on same data. Nothing to do (except sanity checks).
868  KALDI_ASSERT(Trans == kNoTrans && M.NumRows() == NumRows() &&
869  M.NumCols() == NumCols() && M.Stride() == Stride());
870  return;
871  }
872  if (Trans == kNoTrans) {
873  KALDI_ASSERT(num_rows_ == M.NumRows() && num_cols_ == M.NumCols());
874  for (MatrixIndexT i = 0; i < num_rows_; i++)
875  (*this).Row(i).CopyFromVec(M.Row(i));
876  } else {
877  KALDI_ASSERT(num_cols_ == M.NumRows() && num_rows_ == M.NumCols());
878  int32 this_stride = stride_, other_stride = M.Stride();
879  Real *this_data = data_;
880  const OtherReal *other_data = M.Data();
881  for (MatrixIndexT i = 0; i < num_rows_; i++)
882  for (MatrixIndexT j = 0; j < num_cols_; j++)
883  this_data[i * this_stride + j] = other_data[j * other_stride + i];
884  }
885 }
886 
887 // template instantiations.
888 template
890  MatrixTransposeType Trans);
891 template
893  MatrixTransposeType Trans);
894 template
895 void MatrixBase<float>::CopyFromMat(const MatrixBase<float> & M,
896  MatrixTransposeType Trans);
897 template
898 void MatrixBase<double>::CopyFromMat(const MatrixBase<double> & M,
899  MatrixTransposeType Trans);
900 
901 // Specialize the template for CopyFromSp for float, float.
902 template<>
903 template<>
905  KALDI_ASSERT(num_rows_ == M.NumRows() && num_cols_ == num_rows_);
906  MatrixIndexT num_rows = num_rows_, stride = stride_;
907  const float *Mdata = M.Data();
908  float *row_data = data_, *col_data = data_;
909  for (MatrixIndexT i = 0; i < num_rows; i++) {
910  cblas_scopy(i+1, Mdata, 1, row_data, 1); // copy to the row.
911  cblas_scopy(i, Mdata, 1, col_data, stride); // copy to the column.
912  Mdata += i+1;
913  row_data += stride;
914  col_data += 1;
915  }
916 }
917 
918 // Specialize the template for CopyFromSp for double, double.
919 template<>
920 template<>
922  KALDI_ASSERT(num_rows_ == M.NumRows() && num_cols_ == num_rows_);
923  MatrixIndexT num_rows = num_rows_, stride = stride_;
924  const double *Mdata = M.Data();
925  double *row_data = data_, *col_data = data_;
926  for (MatrixIndexT i = 0; i < num_rows; i++) {
927  cblas_dcopy(i+1, Mdata, 1, row_data, 1); // copy to the row.
928  cblas_dcopy(i, Mdata, 1, col_data, stride); // copy to the column.
929  Mdata += i+1;
930  row_data += stride;
931  col_data += 1;
932  }
933 }
934 
935 
936 template<typename Real>
937 template<typename OtherReal>
939  KALDI_ASSERT(num_rows_ == M.NumRows() && num_cols_ == num_rows_);
940  // MORE EFFICIENT IF LOWER TRIANGULAR! Reverse code otherwise.
941  for (MatrixIndexT i = 0; i < num_rows_; i++) {
942  for (MatrixIndexT j = 0; j < i; j++) {
943  (*this)(j, i) = (*this)(i, j) = M(i, j);
944  }
945  (*this)(i, i) = M(i, i);
946  }
947 }
948 
949 // Instantiate this function
950 template
952 template
954 
955 
956 template<typename Real>
957 template<typename OtherReal>
959  MatrixTransposeType Trans) {
960  if (Trans == kNoTrans) {
961  KALDI_ASSERT(num_rows_ == M.NumRows() && num_cols_ == num_rows_);
962  SetZero();
963  Real *out_i = data_;
964  const OtherReal *in_i = M.Data();
965  for (MatrixIndexT i = 0; i < num_rows_; i++, out_i += stride_, in_i += i) {
966  for (MatrixIndexT j = 0; j <= i; j++)
967  out_i[j] = in_i[j];
968  }
969  } else {
970  SetZero();
971  KALDI_ASSERT(num_rows_ == M.NumRows() && num_cols_ == num_rows_);
972  MatrixIndexT stride = stride_;
973  Real *out_i = data_;
974  const OtherReal *in_i = M.Data();
975  for (MatrixIndexT i = 0; i < num_rows_; i++, out_i ++, in_i += i) {
976  for (MatrixIndexT j = 0; j <= i; j++)
977  out_i[j*stride] = in_i[j];
978  }
979  }
980 }
981 
982 template
984  MatrixTransposeType trans);
985 template
986 void MatrixBase<float>::CopyFromTp(const TpMatrix<double> & M,
987  MatrixTransposeType trans);
988 template
990  MatrixTransposeType trans);
991 template
992 void MatrixBase<double>::CopyFromTp(const TpMatrix<double> & M,
993  MatrixTransposeType trans);
994 
995 
996 template<typename Real>
998  if (rv.Dim() == num_rows_*num_cols_) {
999  if (stride_ == num_cols_) {
1000  // one big copy operation.
1001  const Real *rv_data = rv.Data();
1002  std::memcpy(data_, rv_data, sizeof(Real)*num_rows_*num_cols_);
1003  } else {
1004  const Real *rv_data = rv.Data();
1005  for (MatrixIndexT r = 0; r < num_rows_; r++) {
1006  Real *row_data = RowData(r);
1007  for (MatrixIndexT c = 0; c < num_cols_; c++) {
1008  row_data[c] = rv_data[c];
1009  }
1010  rv_data += num_cols_;
1011  }
1012  }
1013  } else if (rv.Dim() == num_cols_) {
1014  const Real *rv_data = rv.Data();
1015  for (MatrixIndexT r = 0; r < num_rows_; r++)
1016  std::memcpy(RowData(r), rv_data, sizeof(Real)*num_cols_);
1017  } else {
1018  KALDI_ERR << "Wrong sized arguments";
1019  }
1020 }
1021 
1022 template<typename Real>
1023 template<typename OtherReal>
1025  if (rv.Dim() == num_rows_*num_cols_) {
1026  const OtherReal *rv_data = rv.Data();
1027  for (MatrixIndexT r = 0; r < num_rows_; r++) {
1028  Real *row_data = RowData(r);
1029  for (MatrixIndexT c = 0; c < num_cols_; c++) {
1030  row_data[c] = static_cast<Real>(rv_data[c]);
1031  }
1032  rv_data += num_cols_;
1033  }
1034  } else if (rv.Dim() == num_cols_) {
1035  const OtherReal *rv_data = rv.Data();
1036  Real *first_row_data = RowData(0);
1037  for (MatrixIndexT c = 0; c < num_cols_; c++)
1038  first_row_data[c] = rv_data[c];
1039  for (MatrixIndexT r = 1; r < num_rows_; r++)
1040  std::memcpy(RowData(r), first_row_data, sizeof(Real)*num_cols_);
1041  } else {
1042  KALDI_ERR << "Wrong sized arguments.";
1043  }
1044 }
1045 
1046 
1047 template
1049 template
1051 
1052 template<typename Real>
1054  if (rv.Dim() == num_rows_*num_cols_) {
1055  const Real *v_inc_data = rv.Data();
1056  Real *m_inc_data = data_;
1057 
1058  for (MatrixIndexT c = 0; c < num_cols_; c++) {
1059  for (MatrixIndexT r = 0; r < num_rows_; r++) {
1060  m_inc_data[r * stride_] = v_inc_data[r];
1061  }
1062  v_inc_data += num_rows_;
1063  m_inc_data ++;
1064  }
1065  } else if (rv.Dim() == num_rows_) {
1066  const Real *v_inc_data = rv.Data();
1067  Real *m_inc_data = data_;
1068  for (MatrixIndexT r = 0; r < num_rows_; r++) {
1069  Real value = *(v_inc_data++);
1070  for (MatrixIndexT c = 0; c < num_cols_; c++)
1071  m_inc_data[c] = value;
1072  m_inc_data += stride_;
1073  }
1074  } else {
1075  KALDI_ERR << "Wrong size of arguments.";
1076  }
1077 }
1078 
1079 
1080 template<typename Real>
1082  KALDI_ASSERT(rv.Dim() == num_cols_ &&
1083  static_cast<UnsignedMatrixIndexT>(row) <
1084  static_cast<UnsignedMatrixIndexT>(num_rows_));
1085 
1086  const Real *rv_data = rv.Data();
1087  Real *row_data = RowData(row);
1088 
1089  std::memcpy(row_data, rv_data, num_cols_ * sizeof(Real));
1090 }
1091 
1092 template<typename Real>
1094  KALDI_ASSERT(rv.Dim() == std::min(num_cols_, num_rows_));
1095  const Real *rv_data = rv.Data(), *rv_end = rv_data + rv.Dim();
1096  Real *my_data = this->Data();
1097  for (; rv_data != rv_end; rv_data++, my_data += (this->stride_+1))
1098  *my_data = *rv_data;
1099 }
1100 
1101 template<typename Real>
1103  const MatrixIndexT col) {
1104  KALDI_ASSERT(rv.Dim() == num_rows_ &&
1105  static_cast<UnsignedMatrixIndexT>(col) <
1106  static_cast<UnsignedMatrixIndexT>(num_cols_));
1107 
1108  const Real *rv_data = rv.Data();
1109  Real *col_data = data_ + col;
1110 
1111  for (MatrixIndexT r = 0; r < num_rows_; r++)
1112  col_data[r * stride_] = rv_data[r];
1113 }
1114 
1115 
1116 
1117 template<typename Real>
1119  KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(i) <
1120  static_cast<UnsignedMatrixIndexT>(MatrixBase<Real>::num_rows_)
1121  && "Access out of matrix");
1122  for (MatrixIndexT j = i + 1; j < MatrixBase<Real>::num_rows_; j++)
1123  MatrixBase<Real>::Row(j-1).CopyFromVec( MatrixBase<Real>::Row(j));
1125 }
1126 
1127 template<typename Real>
1129  // we need to free the data block if it was defined
1130  if (NULL != MatrixBase<Real>::data_)
1132  MatrixBase<Real>::data_ = NULL;
1135 }
1136 
1137 
1138 
1139 template<typename Real>
1141  KALDI_ASSERT(a.NumRows() == num_rows_ && a.NumCols() == num_cols_);
1142 
1143  if (num_cols_ == stride_ && num_cols_ == a.stride_) {
1144  mul_elements(num_rows_ * num_cols_, a.data_, data_);
1145  } else {
1146  MatrixIndexT a_stride = a.stride_, stride = stride_;
1147  Real *data = data_, *a_data = a.data_;
1148  for (MatrixIndexT i = 0; i < num_rows_; i++) {
1149  mul_elements(num_cols_, a_data, data);
1150  a_data += a_stride;
1151  data += stride;
1152  }
1153  }
1154 }
1155 
1156 template<typename Real>
1158  KALDI_ASSERT(a.NumRows() == num_rows_ && a.NumCols() == num_cols_);
1159  MatrixIndexT i;
1160  MatrixIndexT j;
1161 
1162  for (i = 0; i < num_rows_; i++) {
1163  for (j = 0; j < num_cols_; j++) {
1164  (*this)(i, j) /= a(i, j);
1165  }
1166  }
1167 }
1168 
1169 template<typename Real>
1171  double sum = 0.0;
1172 
1173  for (MatrixIndexT i = 0; i < num_rows_; i++) {
1174  for (MatrixIndexT j = 0; j < num_cols_; j++) {
1175  sum += (*this)(i, j);
1176  }
1177  }
1178 
1179  return (Real)sum;
1180 }
1181 
1182 template<typename Real> void MatrixBase<Real>::Max(const MatrixBase<Real> &A) {
1183  KALDI_ASSERT(A.NumRows() == NumRows() && A.NumCols() == NumCols());
1184  for (MatrixIndexT row = 0; row < num_rows_; row++) {
1185  Real *row_data = RowData(row);
1186  const Real *other_row_data = A.RowData(row);
1187  MatrixIndexT num_cols = num_cols_;
1188  for (MatrixIndexT col = 0; col < num_cols; col++) {
1189  row_data[col] = std::max(row_data[col],
1190  other_row_data[col]);
1191  }
1192  }
1193 }
1194 
1195 template<typename Real> void MatrixBase<Real>::Min(const MatrixBase<Real> &A) {
1196  KALDI_ASSERT(A.NumRows() == NumRows() && A.NumCols() == NumCols());
1197  for (MatrixIndexT row = 0; row < num_rows_; row++) {
1198  Real *row_data = RowData(row);
1199  const Real *other_row_data = A.RowData(row);
1200  MatrixIndexT num_cols = num_cols_;
1201  for (MatrixIndexT col = 0; col < num_cols; col++) {
1202  row_data[col] = std::min(row_data[col],
1203  other_row_data[col]);
1204  }
1205  }
1206 }
1207 
1208 
1209 template<typename Real> void MatrixBase<Real>::Scale(Real alpha) {
1210  if (alpha == 1.0) return;
1211  if (num_rows_ == 0) return;
1212  if (num_cols_ == stride_) {
1213  cblas_Xscal(static_cast<size_t>(num_rows_) * static_cast<size_t>(num_cols_),
1214  alpha, data_,1);
1215  } else {
1216  Real *data = data_;
1217  for (MatrixIndexT i = 0; i < num_rows_; ++i, data += stride_) {
1218  cblas_Xscal(num_cols_, alpha, data,1);
1219  }
1220  }
1221 }
1222 
1223 template<typename Real> // scales each row by scale[i].
1225  KALDI_ASSERT(scale.Dim() == num_rows_);
1226  MatrixIndexT M = num_rows_, N = num_cols_;
1227 
1228  for (MatrixIndexT i = 0; i < M; i++) {
1229  Real this_scale = scale(i);
1230  for (MatrixIndexT j = 0; j < N; j++) {
1231  (*this)(i, j) *= this_scale;
1232  }
1233  }
1234 }
1235 
1236 
1237 template<typename Real>
1239  KALDI_ASSERT(src.NumRows() == this->NumRows() &&
1240  this->NumCols() % src.NumCols() == 0);
1241  int32 group_size = this->NumCols() / src.NumCols(),
1242  num_groups = this->NumCols() / group_size,
1243  num_rows = this->NumRows();
1244 
1245  for (MatrixIndexT i = 0; i < num_rows; i++) {
1246  Real *data = this->RowData(i);
1247  for (MatrixIndexT j = 0; j < num_groups; j++, data += group_size) {
1248  Real scale = src(i, j);
1249  cblas_Xscal(group_size, scale, data, 1);
1250  }
1251  }
1252 }
1253 
1254 template<typename Real>
1256  const MatrixBase<Real> &output,
1257  Real power) {
1258  KALDI_ASSERT(input.NumCols() == this->NumCols() && input.NumRows() == this->NumRows());
1259  KALDI_ASSERT(this->NumCols() % output.NumCols() == 0 &&
1260  this->NumRows() == output.NumRows());
1261 
1262  int group_size = this->NumCols() / output.NumCols(),
1263  num_rows = this->NumRows(), num_cols = this->NumCols();
1264 
1265  if (power == 1.0) {
1266  for (MatrixIndexT i = 0; i < num_rows; i++) {
1267  for (MatrixIndexT j = 0; j < num_cols; j++) {
1268  Real input_val = input(i, j);
1269  (*this)(i, j) = (input_val == 0 ? 0 : (input_val > 0 ? 1 : -1));
1270  }
1271  }
1272  } else if (power == std::numeric_limits<Real>::infinity()) {
1273  for (MatrixIndexT i = 0; i < num_rows; i++) {
1274  for (MatrixIndexT j = 0; j < num_cols; j++) {
1275  Real output_val = output(i, j / group_size), input_val = input(i, j);
1276  if (output_val == 0)
1277  (*this)(i, j) = 0;
1278  else
1279  (*this)(i, j) = (std::abs(input_val) == output_val ? 1.0 : 0.0)
1280  * (input_val >= 0 ? 1 : -1);
1281  }
1282  }
1283  } else {
1284  for (MatrixIndexT i = 0; i < num_rows; i++) {
1285  for (MatrixIndexT j = 0; j < num_cols; j++) {
1286  Real output_val = output(i, j / group_size),
1287  input_val = input(i, j);
1288  if (output_val == 0)
1289  (*this)(i, j) = 0;
1290  else
1291  (*this)(i, j) = pow(std::abs(input_val), power - 1) *
1292  pow(output_val, 1 - power) * (input_val >= 0 ? 1 : -1) ;
1293  }
1294  }
1295  }
1296 }
1297 
1298 template<typename Real>
1300  const MatrixBase<Real> &output) {
1301  KALDI_ASSERT(input.NumCols() == this->NumCols() &&
1302  input.NumRows() == this->NumRows());
1303  KALDI_ASSERT(this->NumCols() % output.NumCols() == 0 &&
1304  this->NumRows() == output.NumRows());
1305 
1306  int group_size = this->NumCols() / output.NumCols(),
1307  num_rows = this->NumRows(), num_cols = this->NumCols();
1308 
1309  for (MatrixIndexT i = 0; i < num_rows; i++) {
1310  for (MatrixIndexT j = 0; j < num_cols; j++) {
1311  Real input_val = input(i, j);
1312  Real output_val = output(i, j / group_size);
1313  (*this)(i, j) = (input_val == output_val ? 1 : 0);
1314  }
1315  }
1316 }
1317 
1318 template<typename Real> // scales each column by scale[i].
1320  KALDI_ASSERT(scale.Dim() == num_cols_);
1321  for (MatrixIndexT i = 0; i < num_rows_; i++) {
1322  for (MatrixIndexT j = 0; j < num_cols_; j++) {
1323  Real this_scale = scale(j);
1324  (*this)(i, j) *= this_scale;
1325  }
1326  }
1327 }
1328 
1329 template<typename Real>
1331  if (num_cols_ == stride_)
1332  memset(data_, 0, sizeof(Real)*num_rows_*num_cols_);
1333  else
1334  for (MatrixIndexT row = 0; row < num_rows_; row++)
1335  memset(data_ + row*stride_, 0, sizeof(Real)*num_cols_);
1336 }
1337 
1338 template<typename Real>
1339 void MatrixBase<Real>::Set(Real value) {
1340  for (MatrixIndexT row = 0; row < num_rows_; row++) {
1341  for (MatrixIndexT col = 0; col < num_cols_; col++) {
1342  (*this)(row, col) = value;
1343  }
1344  }
1345 }
1346 
1347 template<typename Real>
1349  SetZero();
1350  for (MatrixIndexT row = 0; row < std::min(num_rows_, num_cols_); row++)
1351  (*this)(row, row) = 1.0;
1352 }
1353 
1354 template<typename Real>
1356  kaldi::RandomState rstate;
1357  for (MatrixIndexT row = 0; row < num_rows_; row++) {
1358  Real *row_data = this->RowData(row);
1359  MatrixIndexT nc = (num_cols_ % 2 == 1) ? num_cols_ - 1 : num_cols_;
1360  for (MatrixIndexT col = 0; col < nc; col += 2) {
1361  kaldi::RandGauss2(row_data + col, row_data + col + 1, &rstate);
1362  }
1363  if (nc != num_cols_) row_data[nc] = static_cast<Real>(kaldi::RandGauss(&rstate));
1364  }
1365 }
1366 
1367 template<typename Real>
1369  kaldi::RandomState rstate;
1370  for (MatrixIndexT row = 0; row < num_rows_; row++) {
1371  Real *row_data = this->RowData(row);
1372  for (MatrixIndexT col = 0; col < num_cols_; col++, row_data++) {
1373  *row_data = static_cast<Real>(kaldi::RandUniform(&rstate));
1374  }
1375  }
1376 }
1377 
1378 template<typename Real>
1379 void MatrixBase<Real>::Write(std::ostream &os, bool binary) const {
1380  if (!os.good()) {
1381  KALDI_ERR << "Failed to write matrix to stream: stream not good";
1382  }
1383  if (binary) { // Use separate binary and text formats,
1384  // since in binary mode we need to know if it's float or double.
1385  std::string my_token = (sizeof(Real) == 4 ? "FM" : "DM");
1386 
1387  WriteToken(os, binary, my_token);
1388  {
1389  int32 rows = this->num_rows_; // make the size 32-bit on disk.
1390  int32 cols = this->num_cols_;
1391  KALDI_ASSERT(this->num_rows_ == (MatrixIndexT) rows);
1392  KALDI_ASSERT(this->num_cols_ == (MatrixIndexT) cols);
1393  WriteBasicType(os, binary, rows);
1394  WriteBasicType(os, binary, cols);
1395  }
1396  if (Stride() == NumCols())
1397  os.write(reinterpret_cast<const char*> (Data()), sizeof(Real)
1398  * static_cast<size_t>(num_rows_) * static_cast<size_t>(num_cols_));
1399  else
1400  for (MatrixIndexT i = 0; i < num_rows_; i++)
1401  os.write(reinterpret_cast<const char*> (RowData(i)), sizeof(Real)
1402  * num_cols_);
1403  if (!os.good()) {
1404  KALDI_ERR << "Failed to write matrix to stream";
1405  }
1406  } else { // text mode.
1407  if (num_cols_ == 0) {
1408  os << " [ ]\n";
1409  } else {
1410  os << " [";
1411  for (MatrixIndexT i = 0; i < num_rows_; i++) {
1412  os << "\n ";
1413  for (MatrixIndexT j = 0; j < num_cols_; j++)
1414  os << (*this)(i, j) << " ";
1415  }
1416  os << "]\n";
1417  }
1418  }
1419 }
1420 
1421 
1422 template<typename Real>
1423 void MatrixBase<Real>::Read(std::istream & is, bool binary, bool add) {
1424  if (add) {
1425  Matrix<Real> tmp(num_rows_, num_cols_);
1426  tmp.Read(is, binary, false); // read without adding.
1427  if (tmp.num_rows_ != this->num_rows_ || tmp.num_cols_ != this->num_cols_)
1428  KALDI_ERR << "MatrixBase::Read, size mismatch "
1429  << this->num_rows_ << ", " << this->num_cols_
1430  << " vs. " << tmp.num_rows_ << ", " << tmp.num_cols_;
1431  this->AddMat(1.0, tmp);
1432  return;
1433  }
1434  // now assume add == false.
1435 
1436  // In order to avoid rewriting this, we just declare a Matrix and
1437  // use it to read the data, then copy.
1438  Matrix<Real> tmp;
1439  tmp.Read(is, binary, false);
1440  if (tmp.NumRows() != NumRows() || tmp.NumCols() != NumCols()) {
1441  KALDI_ERR << "MatrixBase<Real>::Read, size mismatch "
1442  << NumRows() << " x " << NumCols() << " versus "
1443  << tmp.NumRows() << " x " << tmp.NumCols();
1444  }
1445  CopyFromMat(tmp);
1446 }
1447 
1448 
1449 template<typename Real>
1450 void Matrix<Real>::Read(std::istream & is, bool binary, bool add) {
1451  if (add) {
1452  Matrix<Real> tmp;
1453  tmp.Read(is, binary, false); // read without adding.
1454  if (this->num_rows_ == 0) this->Resize(tmp.num_rows_, tmp.num_cols_);
1455  else {
1456  if (this->num_rows_ != tmp.num_rows_ || this->num_cols_ != tmp.num_cols_) {
1457  if (tmp.num_rows_ == 0) return; // do nothing in this case.
1458  else KALDI_ERR << "Matrix::Read, size mismatch "
1459  << this->num_rows_ << ", " << this->num_cols_
1460  << " vs. " << tmp.num_rows_ << ", " << tmp.num_cols_;
1461  }
1462  }
1463  this->AddMat(1.0, tmp);
1464  return;
1465  }
1466 
1467  // now assume add == false.
1468  MatrixIndexT pos_at_start = is.tellg();
1469  std::ostringstream specific_error;
1470 
1471  if (binary) { // Read in binary mode.
1472  int peekval = Peek(is, binary);
1473  if (peekval == 'C') {
1474  // This code enables us to read CompressedMatrix as a regular matrix.
1475  CompressedMatrix compressed_mat;
1476  compressed_mat.Read(is, binary); // at this point, add == false.
1477  this->Resize(compressed_mat.NumRows(), compressed_mat.NumCols());
1478  compressed_mat.CopyToMat(this);
1479  return;
1480  }
1481  const char *my_token = (sizeof(Real) == 4 ? "FM" : "DM");
1482  char other_token_start = (sizeof(Real) == 4 ? 'D' : 'F');
1483  if (peekval == other_token_start) { // need to instantiate the other type to read it.
1484  typedef typename OtherReal<Real>::Real OtherType; // if Real == float, OtherType == double, and vice versa.
1485  Matrix<OtherType> other(this->num_rows_, this->num_cols_);
1486  other.Read(is, binary, false); // add is false at this point anyway.
1487  this->Resize(other.NumRows(), other.NumCols());
1488  this->CopyFromMat(other);
1489  return;
1490  }
1491  std::string token;
1492  ReadToken(is, binary, &token);
1493  if (token != my_token) {
1494  if (token.length() > 20) token = token.substr(0, 17) + "...";
1495  specific_error << ": Expected token " << my_token << ", got " << token;
1496  goto bad;
1497  }
1498  int32 rows, cols;
1499  ReadBasicType(is, binary, &rows); // throws on error.
1500  ReadBasicType(is, binary, &cols); // throws on error.
1501  if ((MatrixIndexT)rows != this->num_rows_ || (MatrixIndexT)cols != this->num_cols_) {
1502  this->Resize(rows, cols);
1503  }
1504  if (this->Stride() == this->NumCols() && rows*cols!=0) {
1505  is.read(reinterpret_cast<char*>(this->Data()),
1506  sizeof(Real)*rows*cols);
1507  if (is.fail()) goto bad;
1508  } else {
1509  for (MatrixIndexT i = 0; i < (MatrixIndexT)rows; i++) {
1510  is.read(reinterpret_cast<char*>(this->RowData(i)), sizeof(Real)*cols);
1511  if (is.fail()) goto bad;
1512  }
1513  }
1514  if (is.eof()) return;
1515  if (is.fail()) goto bad;
1516  return;
1517  } else { // Text mode.
1518  std::string str;
1519  is >> str; // get a token
1520  if (is.fail()) { specific_error << ": Expected \"[\", got EOF"; goto bad; }
1521  // if ((str.compare("DM") == 0) || (str.compare("FM") == 0)) { // Back compatibility.
1522  // is >> str; // get #rows
1523  // is >> str; // get #cols
1524  // is >> str; // get "["
1525  // }
1526  if (str == "[]") { Resize(0, 0); return; } // Be tolerant of variants.
1527  else if (str != "[") {
1528  if (str.length() > 20) str = str.substr(0, 17) + "...";
1529  specific_error << ": Expected \"[\", got \"" << str << '"';
1530  goto bad;
1531  }
1532  // At this point, we have read "[".
1533  std::vector<std::vector<Real>* > data;
1534  std::vector<Real> *cur_row = new std::vector<Real>;
1535  while (1) {
1536  int i = is.peek();
1537  if (i == -1) { specific_error << "Got EOF while reading matrix data"; goto cleanup; }
1538  else if (static_cast<char>(i) == ']') { // Finished reading matrix.
1539  is.get(); // eat the "]".
1540  i = is.peek();
1541  if (static_cast<char>(i) == '\r') {
1542  is.get();
1543  is.get(); // get \r\n (must eat what we wrote)
1544  } else if (static_cast<char>(i) == '\n') { is.get(); } // get \n (must eat what we wrote)
1545  if (is.fail()) {
1546  KALDI_WARN << "After end of matrix data, read error.";
1547  // we got the data we needed, so just warn for this error.
1548  }
1549  // Now process the data.
1550  if (!cur_row->empty()) data.push_back(cur_row);
1551  else delete(cur_row);
1552  cur_row = NULL;
1553  if (data.empty()) { this->Resize(0, 0); return; }
1554  else {
1555  int32 num_rows = data.size(), num_cols = data[0]->size();
1556  this->Resize(num_rows, num_cols);
1557  for (int32 i = 0; i < num_rows; i++) {
1558  if (static_cast<int32>(data[i]->size()) != num_cols) {
1559  specific_error << "Matrix has inconsistent #cols: " << num_cols
1560  << " vs." << data[i]->size() << " (processing row"
1561  << i << ")";
1562  goto cleanup;
1563  }
1564  for (int32 j = 0; j < num_cols; j++)
1565  (*this)(i, j) = (*(data[i]))[j];
1566  delete data[i];
1567  data[i] = NULL;
1568  }
1569  }
1570  return;
1571  } else if (static_cast<char>(i) == '\n' || static_cast<char>(i) == ';') {
1572  // End of matrix row.
1573  is.get();
1574  if (cur_row->size() != 0) {
1575  data.push_back(cur_row);
1576  cur_row = new std::vector<Real>;
1577  cur_row->reserve(data.back()->size());
1578  }
1579  } else if ( (i >= '0' && i <= '9') || i == '-' ) { // A number...
1580  Real r;
1581  is >> r;
1582  if (is.fail()) {
1583  specific_error << "Stream failure/EOF while reading matrix data.";
1584  goto cleanup;
1585  }
1586  cur_row->push_back(r);
1587  } else if (isspace(i)) {
1588  is.get(); // eat the space and do nothing.
1589  } else { // NaN or inf or error.
1590  std::string str;
1591  is >> str;
1592  if (!KALDI_STRCASECMP(str.c_str(), "inf") ||
1593  !KALDI_STRCASECMP(str.c_str(), "infinity")) {
1594  cur_row->push_back(std::numeric_limits<Real>::infinity());
1595  KALDI_WARN << "Reading infinite value into matrix.";
1596  } else if (!KALDI_STRCASECMP(str.c_str(), "nan")) {
1597  cur_row->push_back(std::numeric_limits<Real>::quiet_NaN());
1598  KALDI_WARN << "Reading NaN value into matrix.";
1599  } else {
1600  if (str.length() > 20) str = str.substr(0, 17) + "...";
1601  specific_error << "Expecting numeric matrix data, got " << str;
1602  goto cleanup;
1603  }
1604  }
1605  }
1606  // Note, we never leave the while () loop before this
1607  // line (we return from it.)
1608  cleanup: // We only reach here in case of error in the while loop above.
1609  if(cur_row != NULL)
1610  delete cur_row;
1611  for (size_t i = 0; i < data.size(); i++)
1612  if(data[i] != NULL)
1613  delete data[i];
1614  // and then go on to "bad" below, where we print error.
1615  }
1616 bad:
1617  KALDI_ERR << "Failed to read matrix from stream. " << specific_error.str()
1618  << " File position at start is "
1619  << pos_at_start << ", currently " << is.tellg();
1620 }
1621 
1622 
1623 // Constructor... note that this is not const-safe as it would
1624 // be quite complicated to implement a "const SubMatrix" class that
1625 // would not allow its contents to be changed.
1626 template<typename Real>
1628  const MatrixIndexT ro,
1629  const MatrixIndexT r,
1630  const MatrixIndexT co,
1631  const MatrixIndexT c) {
1632  if (r == 0 || c == 0) {
1633  // we support the empty sub-matrix as a special case.
1634  KALDI_ASSERT(c == 0 && r == 0);
1635  this->data_ = NULL;
1636  this->num_cols_ = 0;
1637  this->num_rows_ = 0;
1638  this->stride_ = 0;
1639  return;
1640  }
1641  KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(ro) <
1642  static_cast<UnsignedMatrixIndexT>(M.num_rows_) &&
1643  static_cast<UnsignedMatrixIndexT>(co) <
1644  static_cast<UnsignedMatrixIndexT>(M.num_cols_) &&
1645  static_cast<UnsignedMatrixIndexT>(r) <=
1646  static_cast<UnsignedMatrixIndexT>(M.num_rows_ - ro) &&
1647  static_cast<UnsignedMatrixIndexT>(c) <=
1648  static_cast<UnsignedMatrixIndexT>(M.num_cols_ - co));
1649  // point to the begining of window
1654  static_cast<size_t>(co) +
1655  static_cast<size_t>(ro) * static_cast<size_t>(M.Stride());
1656 }
1657 
1658 
1659 template<typename Real>
1661  MatrixIndexT num_rows,
1662  MatrixIndexT num_cols,
1663  MatrixIndexT stride):
1664  MatrixBase<Real>(data, num_cols, num_rows, stride) { // caution: reversed order!
1665  if (data == NULL) {
1666  KALDI_ASSERT(num_rows * num_cols == 0);
1667  this->num_rows_ = 0;
1668  this->num_cols_ = 0;
1669  this->stride_ = 0;
1670  } else {
1671  KALDI_ASSERT(this->stride_ >= this->num_cols_);
1672  }
1673 }
1674 
1675 
1676 template<typename Real>
1677 void MatrixBase<Real>::Add(const Real alpha) {
1678  Real *data = data_;
1679  MatrixIndexT stride = stride_;
1680  for (MatrixIndexT r = 0; r < num_rows_; r++)
1681  for (MatrixIndexT c = 0; c < num_cols_; c++)
1682  data[c + stride*r] += alpha;
1683 }
1684 
1685 template<typename Real>
1686 void MatrixBase<Real>::AddToDiag(const Real alpha) {
1687  Real *data = data_;
1688  MatrixIndexT this_stride = stride_ + 1,
1689  num_to_add = std::min(num_rows_, num_cols_);
1690  for (MatrixIndexT r = 0; r < num_to_add; r++)
1691  data[r * this_stride] += alpha;
1692 }
1693 
1694 
1695 template<typename Real>
1697  KALDI_ASSERT(num_rows_ > 0&&num_cols_ > 0);
1698  Vector<Real> singular_values(std::min(num_rows_, num_cols_));
1699  Svd(&singular_values); // Get singular values...
1700  Real min = singular_values(0), max = singular_values(0); // both absolute values...
1701  for (MatrixIndexT i = 1;i < singular_values.Dim();i++) {
1702  min = std::min((Real)std::abs(singular_values(i)), min); max = std::max((Real)std::abs(singular_values(i)), max);
1703  }
1704  if (min > 0) return max/min;
1705  else return std::numeric_limits<Real>::infinity();
1706 }
1707 
1708 template<typename Real>
1709 Real MatrixBase<Real>::Trace(bool check_square) const {
1710  KALDI_ASSERT(!check_square || num_rows_ == num_cols_);
1711  Real ans = 0.0;
1712  for (MatrixIndexT r = 0;r < std::min(num_rows_, num_cols_);r++) ans += data_ [r + stride_*r];
1713  return ans;
1714 }
1715 
1716 template<typename Real>
1718  KALDI_ASSERT(num_rows_ > 0 && num_cols_ > 0);
1719  Real ans= *data_;
1720  for (MatrixIndexT r = 0; r < num_rows_; r++)
1721  for (MatrixIndexT c = 0; c < num_cols_; c++)
1722  if (data_[c + stride_*r] > ans)
1723  ans = data_[c + stride_*r];
1724  return ans;
1725 }
1726 
1727 template<typename Real>
1729  KALDI_ASSERT(num_rows_ > 0 && num_cols_ > 0);
1730  Real ans= *data_;
1731  for (MatrixIndexT r = 0; r < num_rows_; r++)
1732  for (MatrixIndexT c = 0; c < num_cols_; c++)
1733  if (data_[c + stride_*r] < ans)
1734  ans = data_[c + stride_*r];
1735  return ans;
1736 }
1737 
1738 
1739 
1740 template <typename Real>
1742  const MatrixBase<Real> &A, MatrixTransposeType transA,
1743  const MatrixBase<Real> &B, MatrixTransposeType transB,
1744  const MatrixBase<Real> &C, MatrixTransposeType transC,
1745  Real beta) {
1746  // Note on time taken with different orders of computation. Assume not transposed in this /
1747  // discussion. Firstly, normalize expressions using A.NumCols == B.NumRows and B.NumCols == C.NumRows, prefer
1748  // rows where there is a choice.
1749  // time taken for (AB) is: A.NumRows*B.NumRows*C.Rows
1750  // time taken for (AB)C is A.NumRows*C.NumRows*C.Cols
1751  // so this order is A.NumRows*B.NumRows*C.NumRows + A.NumRows*C.NumRows*C.NumCols.
1752 
1753  // time taken for (BC) is: B.NumRows*C.NumRows*C.Cols
1754  // time taken for A(BC) is: A.NumRows*B.NumRows*C.Cols
1755  // so this order is B.NumRows*C.NumRows*C.NumCols + A.NumRows*B.NumRows*C.Cols
1756 
1757  MatrixIndexT ARows = A.num_rows_, ACols = A.num_cols_, BRows = B.num_rows_, BCols = B.num_cols_,
1758  CRows = C.num_rows_, CCols = C.num_cols_;
1759  if (transA == kTrans) std::swap(ARows, ACols);
1760  if (transB == kTrans) std::swap(BRows, BCols);
1761  if (transC == kTrans) std::swap(CRows, CCols);
1762 
1763  MatrixIndexT AB_C_time = ARows*BRows*CRows + ARows*CRows*CCols;
1764  MatrixIndexT A_BC_time = BRows*CRows*CCols + ARows*BRows*CCols;
1765 
1766  if (AB_C_time < A_BC_time) {
1767  Matrix<Real> AB(ARows, BCols);
1768  AB.AddMatMat(1.0, A, transA, B, transB, 0.0); // AB = A * B.
1769  (*this).AddMatMat(alpha, AB, kNoTrans, C, transC, beta);
1770  } else {
1771  Matrix<Real> BC(BRows, CCols);
1772  BC.AddMatMat(1.0, B, transB, C, transC, 0.0); // BC = B * C.
1773  (*this).AddMatMat(alpha, A, transA, BC, kNoTrans, beta);
1774  }
1775 }
1776 
1777 
1778 
1779 
1780 template<typename Real>
1782  // Svd, *this = U*diag(s)*Vt.
1783  // With (*this).num_rows_ == m, (*this).num_cols_ == n,
1784  // Support only skinny Svd with m>=n (NumRows>=NumCols), and zero sizes for U and Vt mean
1785  // we do not want that output. We expect that s.Dim() == m,
1786  // U is either 0 by 0 or m by n, and rv is either 0 by 0 or n by n.
1787  // Throws exception on error.
1788 
1789  KALDI_ASSERT(num_rows_>=num_cols_ && "Svd requires that #rows by >= #cols."); // For compatibility with JAMA code.
1790  KALDI_ASSERT(s->Dim() == num_cols_); // s should be the smaller dim.
1791  KALDI_ASSERT(U == NULL || (U->num_rows_ == num_rows_&&U->num_cols_ == num_cols_));
1792  KALDI_ASSERT(Vt == NULL || (Vt->num_rows_ == num_cols_&&Vt->num_cols_ == num_cols_));
1793 
1794  Real prescale = 1.0;
1795  if ( std::abs((*this)(0, 0) ) < 1.0e-30) { // Very tiny value... can cause problems in Svd.
1796  Real max_elem = LargestAbsElem();
1797  if (max_elem != 0) {
1798  prescale = 1.0 / max_elem;
1799  if (std::abs(prescale) == std::numeric_limits<Real>::infinity()) { prescale = 1.0e+40; }
1800  (*this).Scale(prescale);
1801  }
1802  }
1803 
1804 #if !defined(HAVE_ATLAS) && !defined(USE_KALDI_SVD)
1805  // "S" == skinny Svd (only one we support because of compatibility with Jama one which is only skinny),
1806  // "N"== no eigenvectors wanted.
1807  LapackGesvd(s, U, Vt);
1808 #else
1809  /* if (num_rows_ > 1 && num_cols_ > 1 && (*this)(0, 0) == (*this)(1, 1)
1810  && Max() == Min() && (*this)(0, 0) != 0.0) { // special case that JamaSvd sometimes crashes on.
1811  KALDI_WARN << "Jama SVD crashes on this type of matrix, perturbing it to prevent crash.";
1812  for(int32 i = 0; i < NumRows(); i++)
1813  (*this)(i, i) *= 1.00001;
1814  }*/
1815  bool ans = JamaSvd(s, U, Vt);
1816  if (Vt != NULL) Vt->Transpose(); // possibly to do: change this and also the transpose inside the JamaSvd routine. note, Vt is square.
1817  if (!ans) {
1818  KALDI_ERR << "Error doing Svd"; // This one will be caught.
1819  }
1820 #endif
1821  if (prescale != 1.0) s->Scale(1.0/prescale);
1822 }
1823 
1824 template<typename Real>
1826  try {
1827  if (num_rows_ >= num_cols_) {
1828  Matrix<Real> tmp(*this);
1829  tmp.DestructiveSvd(s, U, Vt);
1830  } else {
1831  Matrix<Real> tmp(*this, kTrans); // transpose of *this.
1832  // rVt will have different dim so cannot transpose in-place --> use a temp matrix.
1833  Matrix<Real> Vt_Trans(Vt ? Vt->num_cols_ : 0, Vt ? Vt->num_rows_ : 0);
1834  // U will be transpose
1835  tmp.DestructiveSvd(s, Vt ? &Vt_Trans : NULL, U);
1836  if (U) U->Transpose();
1837  if (Vt) Vt->CopyFromMat(Vt_Trans, kTrans); // copy with transpose.
1838  }
1839  } catch (...) {
1840  KALDI_ERR << "Error doing Svd (did not converge), first part of matrix is\n"
1841  << SubMatrix<Real>(*this, 0, std::min((MatrixIndexT)10, num_rows_),
1842  0, std::min((MatrixIndexT)10, num_cols_))
1843  << ", min and max are: " << Min() << ", " << Max();
1844  }
1845 }
1846 
1847 template<typename Real>
1848 bool MatrixBase<Real>::IsSymmetric(Real cutoff) const {
1849  MatrixIndexT R = num_rows_, C = num_cols_;
1850  if (R != C) return false;
1851  Real bad_sum = 0.0, good_sum = 0.0;
1852  for (MatrixIndexT i = 0;i < R;i++) {
1853  for (MatrixIndexT j = 0;j < i;j++) {
1854  Real a = (*this)(i, j), b = (*this)(j, i), avg = 0.5*(a+b), diff = 0.5*(a-b);
1855  good_sum += std::abs(avg); bad_sum += std::abs(diff);
1856  }
1857  good_sum += std::abs((*this)(i, i));
1858  }
1859  if (bad_sum > cutoff*good_sum) return false;
1860  return true;
1861 }
1862 
1863 template<typename Real>
1864 bool MatrixBase<Real>::IsDiagonal(Real cutoff) const{
1865  MatrixIndexT R = num_rows_, C = num_cols_;
1866  Real bad_sum = 0.0, good_sum = 0.0;
1867  for (MatrixIndexT i = 0;i < R;i++) {
1868  for (MatrixIndexT j = 0;j < C;j++) {
1869  if (i == j) good_sum += std::abs((*this)(i, j));
1870  else bad_sum += std::abs((*this)(i, j));
1871  }
1872  }
1873  return (!(bad_sum > good_sum * cutoff));
1874 }
1875 
1876 // This does nothing, it's designed to trigger Valgrind errors
1877 // if any memory is uninitialized.
1878 template<typename Real>
1880  MatrixIndexT R = num_rows_, C = num_cols_, positive = 0;
1881  for (MatrixIndexT i = 0; i < R; i++)
1882  for (MatrixIndexT j = 0; j < C; j++)
1883  if ((*this)(i, j) > 0.0) positive++;
1884  if (positive > R * C)
1885  KALDI_ERR << "Error....";
1886 }
1887 
1888 
1889 template<typename Real>
1890 bool MatrixBase<Real>::IsUnit(Real cutoff) const {
1891  MatrixIndexT R = num_rows_, C = num_cols_;
1892  Real bad_max = 0.0;
1893  for (MatrixIndexT i = 0; i < R;i++)
1894  for (MatrixIndexT j = 0; j < C;j++)
1895  bad_max = std::max(bad_max, static_cast<Real>(std::abs( (*this)(i, j) - (i == j?1.0:0.0))));
1896  return (bad_max <= cutoff);
1897 }
1898 
1899 template<typename Real>
1900 bool MatrixBase<Real>::IsZero(Real cutoff)const {
1901  MatrixIndexT R = num_rows_, C = num_cols_;
1902  Real bad_max = 0.0;
1903  for (MatrixIndexT i = 0;i < R;i++)
1904  for (MatrixIndexT j = 0;j < C;j++)
1905  bad_max = std::max(bad_max, static_cast<Real>(std::abs( (*this)(i, j) )));
1906  return (bad_max <= cutoff);
1907 }
1908 
1909 template<typename Real>
1911  return std::sqrt(TraceMatMat(*this, *this, kTrans));
1912 }
1913 
1914 template<typename Real>
1915 bool MatrixBase<Real>::ApproxEqual(const MatrixBase<Real> &other, float tol) const {
1916  if (num_rows_ != other.num_rows_ || num_cols_ != other.num_cols_)
1917  KALDI_ERR << "ApproxEqual: size mismatch.";
1918  Matrix<Real> tmp(*this);
1919  tmp.AddMat(-1.0, other);
1920  return (tmp.FrobeniusNorm() <= static_cast<Real>(tol) *
1921  this->FrobeniusNorm());
1922 }
1923 
1924 template<typename Real>
1925 bool MatrixBase<Real>::Equal(const MatrixBase<Real> &other) const {
1926  if (num_rows_ != other.num_rows_ || num_cols_ != other.num_cols_)
1927  KALDI_ERR << "Equal: size mismatch.";
1928  for (MatrixIndexT i = 0; i < num_rows_; i++)
1929  for (MatrixIndexT j = 0; j < num_cols_; j++)
1930  if ( (*this)(i, j) != other(i, j))
1931  return false;
1932  return true;
1933 }
1934 
1935 
1936 template<typename Real>
1938  MatrixIndexT R = num_rows_, C = num_cols_;
1939  Real largest = 0.0;
1940  for (MatrixIndexT i = 0;i < R;i++)
1941  for (MatrixIndexT j = 0;j < C;j++)
1942  largest = std::max(largest, (Real)std::abs((*this)(i, j)));
1943  return largest;
1944 }
1945 
1946 
1947 template<typename Real>
1949  KALDI_ASSERT(NumRows() <= NumCols());
1950  MatrixIndexT num_rows = num_rows_;
1951  for (MatrixIndexT i = 0; i < num_rows; i++) {
1952  int32 counter = 0;
1953  while (1) {
1954  Real start_prod = VecVec(this->Row(i), this->Row(i));
1955  if (start_prod - start_prod != 0.0 || start_prod == 0.0) {
1956  KALDI_WARN << "Self-product of row " << i << " of matrix is "
1957  << start_prod << ", randomizing.";
1958  this->Row(i).SetRandn();
1959  counter++;
1960  continue; // loop again.
1961  }
1962  for (MatrixIndexT j = 0; j < i; j++) {
1963  Real prod = VecVec(this->Row(i), this->Row(j));
1964  this->Row(i).AddVec(-prod, this->Row(j));
1965  }
1966  Real end_prod = VecVec(this->Row(i), this->Row(i));
1967  if (end_prod <= 0.01 * start_prod) { // We removed
1968  // almost all of the vector during orthogonalization,
1969  // so we have reason to doubt (for roundoff reasons)
1970  // that it's still orthogonal to the other vectors.
1971  // We need to orthogonalize again.
1972  if (end_prod == 0.0) { // Row is exactly zero:
1973  // generate random direction.
1974  this->Row(i).SetRandn();
1975  }
1976  counter++;
1977  if (counter > 100)
1978  KALDI_ERR << "Loop detected while orthogalizing matrix.";
1979  } else {
1980  this->Row(i).Scale(1.0 / std::sqrt(end_prod));
1981  break;
1982  }
1983  }
1984  }
1985 }
1986 
1987 
1988 // Uses Svd to compute the eigenvalue decomposition of a symmetric positive semidefinite
1989 // matrix:
1990 // (*this) = rU * diag(rs) * rU^T, with rU an orthogonal matrix so rU^{-1} = rU^T.
1991 // Does this by computing svd (*this) = U diag(rs) V^T ... answer is just U diag(rs) U^T.
1992 // Throws exception if this failed to within supplied precision (typically because *this was not
1993 // symmetric positive definite).
1994 
1995 template<typename Real>
1996 void MatrixBase<Real>::SymPosSemiDefEig(VectorBase<Real> *rs, MatrixBase<Real> *rU, Real check_thresh) // e.g. check_thresh = 0.001
1997 {
1998  const MatrixIndexT D = num_rows_;
1999 
2001  KALDI_ASSERT(IsSymmetric() && "SymPosSemiDefEig: expecting input to be symmetrical.");
2002  KALDI_ASSERT(rU->num_rows_ == D && rU->num_cols_ == D && rs->Dim() == D);
2003 
2004  Matrix<Real> Vt(D, D);
2005  Svd(rs, rU, &Vt);
2006 
2007  // First just zero any singular values if the column of U and V do not have +ve dot product--
2008  // this may mean we have small negative eigenvalues, and if we zero them the result will be closer to correct.
2009  for (MatrixIndexT i = 0;i < D;i++) {
2010  Real sum = 0.0;
2011  for (MatrixIndexT j = 0;j < D;j++) sum += (*rU)(j, i) * Vt(i, j);
2012  if (sum < 0.0) (*rs)(i) = 0.0;
2013  }
2014 
2015  {
2016  Matrix<Real> tmpU(*rU); Vector<Real> tmps(*rs); tmps.ApplyPow(0.5);
2017  tmpU.MulColsVec(tmps);
2018  SpMatrix<Real> tmpThis(D);
2019  tmpThis.AddMat2(1.0, tmpU, kNoTrans, 0.0);
2020  Matrix<Real> tmpThisFull(tmpThis);
2021  float new_norm = tmpThisFull.FrobeniusNorm();
2022  float old_norm = (*this).FrobeniusNorm();
2023  tmpThisFull.AddMat(-1.0, (*this));
2024 
2025  if (!(old_norm == 0 && new_norm == 0)) {
2026  float diff_norm = tmpThisFull.FrobeniusNorm();
2027  if (std::abs(new_norm-old_norm) > old_norm*check_thresh || diff_norm > old_norm*check_thresh) {
2028  KALDI_WARN << "SymPosSemiDefEig seems to have failed " << diff_norm << " !<< "
2029  << check_thresh << "*" << old_norm << ", maybe matrix was not "
2030  << "positive semi definite. Continuing anyway.";
2031  }
2032  }
2033  }
2034 }
2035 
2036 
2037 template<typename Real>
2038 Real MatrixBase<Real>::LogDet(Real *det_sign) const {
2039  Real log_det;
2040  Matrix<Real> tmp(*this);
2041  tmp.Invert(&log_det, det_sign, false); // false== output not needed (saves some computation).
2042  return log_det;
2043 }
2044 
2045 template<typename Real>
2046 void MatrixBase<Real>::InvertDouble(Real *log_det, Real *det_sign,
2047  bool inverse_needed) {
2048  double log_det_tmp, det_sign_tmp;
2049  Matrix<double> dmat(*this);
2050  dmat.Invert(&log_det_tmp, &det_sign_tmp, inverse_needed);
2051  if (inverse_needed) (*this).CopyFromMat(dmat);
2052  if (log_det) *log_det = log_det_tmp;
2053  if (det_sign) *det_sign = det_sign_tmp;
2054 }
2055 
2056 template<class Real>
2058  mat.CopyToMat(this);
2059 }
2060 
2061 template<class Real>
2063  Resize(M.NumRows(), M.NumCols(), kUndefined);
2064  M.CopyToMat(this);
2065 }
2066 
2067 
2068 
2069 template<typename Real>
2071  for (MatrixIndexT r = 0; r < num_rows_; r++) {
2072  for (MatrixIndexT c = 0; c < num_cols_; c++) {
2073  (*this)(r, c) = static_cast<Real>(1.0 / (*this)(r, c));
2074  }
2075  }
2076 }
2077 
2078 template<typename Real>
2080  KALDI_ASSERT(num_rows_ == num_cols_);
2081  MatrixIndexT M = num_rows_;
2082  for (MatrixIndexT i = 0;i < M;i++)
2083  for (MatrixIndexT j = 0;j < i;j++) {
2084  Real &a = (*this)(i, j), &b = (*this)(j, i);
2085  std::swap(a, b);
2086  }
2087 }
2088 
2089 
2090 template<typename Real>
2092  if (this->num_rows_ != this->num_cols_) {
2093  Matrix<Real> tmp(*this, kTrans);
2094  Resize(this->num_cols_, this->num_rows_);
2095  this->CopyFromMat(tmp);
2096  } else {
2097  (static_cast<MatrixBase<Real>&>(*this)).Transpose();
2098  }
2099 }
2100 
2101 template<typename Real>
2103  KALDI_ASSERT(SameDim(*this, src));
2104  MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
2105  Real *row_data = data_;
2106  const Real *src_row_data = src.Data();
2107  for (MatrixIndexT row = 0; row < num_rows;
2108  row++,row_data += stride_, src_row_data += src.stride_) {
2109  for (MatrixIndexT col = 0; col < num_cols; col++)
2110  row_data[col] = (src_row_data[col] > 0 ? 1.0 : 0.0);
2111  }
2112 }
2113 
2114 template<typename Real>
2116  KALDI_ASSERT(SameDim(*this, src));
2117  MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
2118  Real *row_data = data_;
2119  const Real *src_row_data = src.Data();
2120  for (MatrixIndexT row = 0; row < num_rows;
2121  row++,row_data += stride_, src_row_data += src.stride_) {
2122  for (MatrixIndexT col = 0; col < num_cols; col++)
2123  row_data[col] = kaldi::Exp(src_row_data[col]);
2124  }
2125 }
2126 
2127 template<typename Real>
2128 void MatrixBase<Real>::Pow(const MatrixBase<Real> &src, Real power) {
2129  KALDI_ASSERT(SameDim(*this, src));
2130  MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
2131  Real *row_data = data_;
2132  const Real *src_row_data = src.Data();
2133  for (MatrixIndexT row = 0; row < num_rows;
2134  row++,row_data += stride_, src_row_data += src.stride_) {
2135  for (MatrixIndexT col = 0; col < num_cols; col++) {
2136  row_data[col] = pow(src_row_data[col], power);
2137  }
2138  }
2139 }
2140 
2141 template<typename Real>
2142 void MatrixBase<Real>::PowAbs(const MatrixBase<Real> &src, Real power, bool include_sign) {
2143  KALDI_ASSERT(SameDim(*this, src));
2144  MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
2145  Real *row_data = data_;
2146  const Real *src_row_data = src.Data();
2147  for (MatrixIndexT row = 0; row < num_rows;
2148  row++,row_data += stride_, src_row_data += src.stride_) {
2149  for (MatrixIndexT col = 0; col < num_cols; col ++) {
2150  if (include_sign == true && src_row_data[col] < 0) {
2151  row_data[col] = -pow(std::abs(src_row_data[col]), power);
2152  } else {
2153  row_data[col] = pow(std::abs(src_row_data[col]), power);
2154  }
2155  }
2156  }
2157 }
2158 
2159 template<typename Real>
2160 void MatrixBase<Real>::Floor(const MatrixBase<Real> &src, Real floor_val) {
2161  KALDI_ASSERT(SameDim(*this, src));
2162  MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
2163  Real *row_data = data_;
2164  const Real *src_row_data = src.Data();
2165  for (MatrixIndexT row = 0; row < num_rows;
2166  row++,row_data += stride_, src_row_data += src.stride_) {
2167  for (MatrixIndexT col = 0; col < num_cols; col++)
2168  row_data[col] = (src_row_data[col] < floor_val ? floor_val : src_row_data[col]);
2169  }
2170 }
2171 
2172 template<typename Real>
2173 void MatrixBase<Real>::Ceiling(const MatrixBase<Real> &src, Real ceiling_val) {
2174  KALDI_ASSERT(SameDim(*this, src));
2175  MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
2176  Real *row_data = data_;
2177  const Real *src_row_data = src.Data();
2178  for (MatrixIndexT row = 0; row < num_rows;
2179  row++,row_data += stride_, src_row_data += src.stride_) {
2180  for (MatrixIndexT col = 0; col < num_cols; col++)
2181  row_data[col] = (src_row_data[col] > ceiling_val ? ceiling_val : src_row_data[col]);
2182  }
2183 }
2184 
2185 template<typename Real>
2187  KALDI_ASSERT(SameDim(*this, src));
2188  MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
2189  Real *row_data = data_;
2190  const Real *src_row_data = src.Data();
2191  for (MatrixIndexT row = 0; row < num_rows;
2192  row++,row_data += stride_, src_row_data += src.stride_) {
2193  for (MatrixIndexT col = 0; col < num_cols; col++)
2194  row_data[col] = kaldi::Log(src_row_data[col]);
2195  }
2196 }
2197 
2198 template<typename Real>
2200  KALDI_ASSERT(SameDim(*this, src));
2201  MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
2202  Real *row_data = data_;
2203  const Real *src_row_data = src.Data();
2204  for (MatrixIndexT row = 0; row < num_rows;
2205  row++,row_data += stride_, src_row_data += src.stride_) {
2206  for (MatrixIndexT col = 0; col < num_cols; col++)
2207  row_data[col] = (src_row_data[col] < Real(0) ? kaldi::Exp(src_row_data[col]) : (src_row_data[col] + Real(1)));
2208  }
2209 }
2210 
2211 template<typename Real>
2212 void MatrixBase<Real>::ExpLimited(const MatrixBase<Real> &src, Real lower_limit, Real upper_limit) {
2213  KALDI_ASSERT(SameDim(*this, src));
2214  MatrixIndexT num_rows = num_rows_, num_cols = num_cols_;
2215  Real *row_data = data_;
2216  const Real *src_row_data = src.Data();
2217  for (MatrixIndexT row = 0; row < num_rows;
2218  row++,row_data += stride_, src_row_data += src.stride_) {
2219  for (MatrixIndexT col = 0; col < num_cols; col++) {
2220  const Real x = src_row_data[col];
2221  if (!(x >= lower_limit))
2222  row_data[col] = kaldi::Exp(lower_limit);
2223  else if (x > upper_limit)
2224  row_data[col] = kaldi::Exp(upper_limit);
2225  else
2226  row_data[col] = kaldi::Exp(x);
2227  }
2228  }
2229 }
2230 
2231 template<typename Real>
2232 bool MatrixBase<Real>::Power(Real power) {
2233  KALDI_ASSERT(num_rows_ > 0 && num_rows_ == num_cols_);
2234  MatrixIndexT n = num_rows_;
2235  Matrix<Real> P(n, n);
2236  Vector<Real> re(n), im(n);
2237  this->Eig(&P, &re, &im);
2238  // Now attempt to take the complex eigenvalues to this power.
2239  for (MatrixIndexT i = 0; i < n; i++)
2240  if (!AttemptComplexPower(&(re(i)), &(im(i)), power))
2241  return false; // e.g. real and negative, or zero, eigenvalues.
2242 
2243  Matrix<Real> D(n, n); // D to the power.
2244  CreateEigenvalueMatrix(re, im, &D);
2245 
2246  Matrix<Real> tmp(n, n); // P times D
2247  tmp.AddMatMat(1.0, P, kNoTrans, D, kNoTrans, 0.0); // tmp := P*D
2248  P.Invert();
2249  // next line is: *this = tmp * P^{-1} = P * D * P^{-1}
2250  (*this).AddMatMat(1.0, tmp, kNoTrans, P, kNoTrans, 0.0);
2251  return true;
2252 }
2253 
2254 template<typename Real>
2256  std::swap(this->data_, other->data_);
2257  std::swap(this->num_cols_, other->num_cols_);
2258  std::swap(this->num_rows_, other->num_rows_);
2259  std::swap(this->stride_, other->stride_);
2260 }
2261 
2262 // Repeating this comment that appeared in the header:
2263 // Eigenvalue Decomposition of a square NxN matrix into the form (*this) = P D
2264 // P^{-1}. Be careful: the relationship of D to the eigenvalues we output is
2265 // slightly complicated, due to the need for P to be real. In the symmetric
2266 // case D is diagonal and real, but in
2267 // the non-symmetric case there may be complex-conjugate pairs of eigenvalues.
2268 // In this case, for the equation (*this) = P D P^{-1} to hold, D must actually
2269 // be block diagonal, with 2x2 blocks corresponding to any such pairs. If a
2270 // pair is lambda +- i*mu, D will have a corresponding 2x2 block
2271 // [lambda, mu; -mu, lambda].
2272 // Note that if the input matrix (*this) is non-invertible, P may not be invertible
2273 // so in this case instead of the equation (*this) = P D P^{-1} holding, we have
2274 // instead (*this) P = P D.
2275 //
2276 // By making the pointer arguments non-NULL or NULL, the user can choose to take
2277 // not to take the eigenvalues directly, and/or the matrix D which is block-diagonal
2278 // with 2x2 blocks.
2279 template<typename Real>
2281  VectorBase<Real> *r,
2282  VectorBase<Real> *i) const {
2283  EigenvalueDecomposition<Real> eig(*this);
2284  if (P) eig.GetV(P);
2285  if (r) eig.GetRealEigenvalues(r);
2286  if (i) eig.GetImagEigenvalues(i);
2287 }
2288 
2289 
2290 // Begin non-member function definitions.
2291 
2292 // /**
2293 // * @brief Extension of the HTK header
2294 // */
2295 // struct HtkHeaderExt
2296 // {
2297 // INT_32 mHeaderSize;
2298 // INT_32 mVersion;
2299 // INT_32 mSampSize;
2300 // };
2301 
2302 template<typename Real>
2303 bool ReadHtk(std::istream &is, Matrix<Real> *M_ptr, HtkHeader *header_ptr)
2304 {
2305  // check instantiated with double or float.
2307  Matrix<Real> &M = *M_ptr;
2308  HtkHeader htk_hdr;
2309 
2310  // TODO(arnab): this fails if the HTK file has CRC cheksum or is compressed.
2311  is.read((char*)&htk_hdr, sizeof(htk_hdr)); // we're being really POSIX here!
2312  if (is.fail()) {
2313  KALDI_WARN << "Could not read header from HTK feature file ";
2314  return false;
2315  }
2316 
2317  KALDI_SWAP4(htk_hdr.mNSamples);
2318  KALDI_SWAP4(htk_hdr.mSamplePeriod);
2319  KALDI_SWAP2(htk_hdr.mSampleSize);
2320  KALDI_SWAP2(htk_hdr.mSampleKind);
2321 
2322  bool has_checksum = false;
2323  {
2324  // See HParm.h in HTK code for sources of these things.
2325  enum BaseParmKind{
2326  Waveform, Lpc, Lprefc, Lpcepstra, Lpdelcep,
2327  Irefc, Mfcc, Fbank, Melspec, User, Discrete, Plp, Anon };
2328 
2329  const int32 IsCompressed = 02000, HasChecksum = 010000, HasVq = 040000,
2330  Problem = IsCompressed | HasVq;
2331  int32 base_parm = htk_hdr.mSampleKind & (077);
2332  has_checksum = (base_parm & HasChecksum) != 0;
2333  htk_hdr.mSampleKind &= ~HasChecksum; // We don't support writing with
2334  // checksum so turn it off.
2335  if (htk_hdr.mSampleKind & Problem)
2336  KALDI_ERR << "Code to read HTK features does not support compressed "
2337  "features, or features with VQ.";
2338  if (base_parm == Waveform || base_parm == Irefc || base_parm == Discrete)
2339  KALDI_ERR << "Attempting to read HTK features from unsupported type "
2340  "(e.g. waveform or discrete features.";
2341  }
2342 
2343  KALDI_VLOG(3) << "HTK header: Num Samples: " << htk_hdr.mNSamples
2344  << "; Sample period: " << htk_hdr.mSamplePeriod
2345  << "; Sample size: " << htk_hdr.mSampleSize
2346  << "; Sample kind: " << htk_hdr.mSampleKind;
2347 
2348  M.Resize(htk_hdr.mNSamples, htk_hdr.mSampleSize / sizeof(float));
2349 
2350  MatrixIndexT i;
2351  MatrixIndexT j;
2352  if (sizeof(Real) == sizeof(float)) {
2353  for (i = 0; i< M.NumRows(); i++) {
2354  is.read((char*)M.RowData(i), sizeof(float)*M.NumCols());
2355  if (is.fail()) {
2356  KALDI_WARN << "Could not read data from HTK feature file ";
2357  return false;
2358  }
2359  if (MachineIsLittleEndian()) {
2360  MatrixIndexT C = M.NumCols();
2361  for (j = 0; j < C; j++) {
2362  KALDI_SWAP4((M(i, j))); // The HTK standard is big-endian!
2363  }
2364  }
2365  }
2366  } else {
2367  float *pmem = new float[M.NumCols()];
2368  for (i = 0; i < M.NumRows(); i++) {
2369  is.read((char*)pmem, sizeof(float)*M.NumCols());
2370  if (is.fail()) {
2371  KALDI_WARN << "Could not read data from HTK feature file ";
2372  delete [] pmem;
2373  return false;
2374  }
2375  MatrixIndexT C = M.NumCols();
2376  for (j = 0; j < C; j++) {
2377  if (MachineIsLittleEndian()) // HTK standard is big-endian!
2378  KALDI_SWAP4(pmem[j]);
2379  M(i, j) = static_cast<Real>(pmem[j]);
2380  }
2381  }
2382  delete [] pmem;
2383  }
2384  if (header_ptr) *header_ptr = htk_hdr;
2385  if (has_checksum) {
2386  int16 checksum;
2387  is.read((char*)&checksum, sizeof(checksum));
2388  if (is.fail())
2389  KALDI_WARN << "Could not read checksum from HTK feature file ";
2390  // We ignore the checksum.
2391  }
2392  return true;
2393 }
2394 
2395 
2396 template
2397 bool ReadHtk(std::istream &is, Matrix<float> *M, HtkHeader *header_ptr);
2398 
2399 template
2400 bool ReadHtk(std::istream &is, Matrix<double> *M, HtkHeader *header_ptr);
2401 
2402 template<typename Real>
2403 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.
2404 {
2405  KALDI_ASSERT(M.NumRows() == static_cast<MatrixIndexT>(htk_hdr.mNSamples));
2406  KALDI_ASSERT(M.NumCols() == static_cast<MatrixIndexT>(htk_hdr.mSampleSize) /
2407  static_cast<MatrixIndexT>(sizeof(float)));
2408 
2409  KALDI_SWAP4(htk_hdr.mNSamples);
2410  KALDI_SWAP4(htk_hdr.mSamplePeriod);
2411  KALDI_SWAP2(htk_hdr.mSampleSize);
2412  KALDI_SWAP2(htk_hdr.mSampleKind);
2413 
2414  os.write((char*)&htk_hdr, sizeof(htk_hdr));
2415  if (os.fail()) goto bad;
2416 
2417  MatrixIndexT i;
2418  MatrixIndexT j;
2419  if (sizeof(Real) == sizeof(float) && !MachineIsLittleEndian()) {
2420  for (i = 0; i< M.NumRows(); i++) { // Unlikely to reach here ever!
2421  os.write((char*)M.RowData(i), sizeof(float)*M.NumCols());
2422  if (os.fail()) goto bad;
2423  }
2424  } else {
2425  float *pmem = new float[M.NumCols()];
2426 
2427  for (i = 0; i < M.NumRows(); i++) {
2428  const Real *rowData = M.RowData(i);
2429  for (j = 0;j < M.NumCols();j++)
2430  pmem[j] = static_cast<float> ( rowData[j] );
2431  if (MachineIsLittleEndian())
2432  for (j = 0;j < M.NumCols();j++)
2433  KALDI_SWAP4(pmem[j]);
2434  os.write((char*)pmem, sizeof(float)*M.NumCols());
2435  if (os.fail()) {
2436  delete [] pmem;
2437  goto bad;
2438  }
2439  }
2440  delete [] pmem;
2441  }
2442  return true;
2443 bad:
2444  KALDI_WARN << "Could not write to HTK feature file ";
2445  return false;
2446 }
2447 
2448 template
2449 bool WriteHtk(std::ostream &os, const MatrixBase<float> &M, HtkHeader htk_hdr);
2450 
2451 template
2452 bool WriteHtk(std::ostream &os, const MatrixBase<double> &M, HtkHeader htk_hdr);
2453 
2454 template<class Real>
2455 bool WriteSphinx(std::ostream &os, const MatrixBase<Real> &M)
2456 {
2457  // CMUSphinx mfc file header contains count of the floats, followed
2458  // by the data in float little endian format.
2459 
2460  int size = M.NumRows() * M.NumCols();
2461  os.write((char*)&size, sizeof(int));
2462  if (os.fail()) goto bad;
2463 
2464  MatrixIndexT i;
2465  MatrixIndexT j;
2466  if (sizeof(Real) == sizeof(float) && MachineIsLittleEndian()) {
2467  for (i = 0; i< M.NumRows(); i++) { // Unlikely to reach here ever!
2468  os.write((char*)M.RowData(i), sizeof(float)*M.NumCols());
2469  if (os.fail()) goto bad;
2470  }
2471  } else {
2472  float *pmem = new float[M.NumCols()];
2473 
2474  for (i = 0; i < M.NumRows(); i++) {
2475  const Real *rowData = M.RowData(i);
2476  for (j = 0;j < M.NumCols();j++)
2477  pmem[j] = static_cast<float> ( rowData[j] );
2478  if (!MachineIsLittleEndian())
2479  for (j = 0;j < M.NumCols();j++)
2480  KALDI_SWAP4(pmem[j]);
2481  os.write((char*)pmem, sizeof(float)*M.NumCols());
2482  if (os.fail()) {
2483  delete [] pmem;
2484  goto bad;
2485  }
2486  }
2487  delete [] pmem;
2488  }
2489  return true;
2490 bad:
2491  KALDI_WARN << "Could not write to Sphinx feature file";
2492  return false;
2493 }
2494 
2495 template
2496 bool WriteSphinx(std::ostream &os, const MatrixBase<float> &M);
2497 
2498 template
2499 bool WriteSphinx(std::ostream &os, const MatrixBase<double> &M);
2500 
2501 template <typename Real>
2503  const MatrixBase<Real> &B, MatrixTransposeType transB,
2504  const MatrixBase<Real> &C, MatrixTransposeType transC) {
2505  MatrixIndexT ARows = A.NumRows(), ACols = A.NumCols(), BRows = B.NumRows(), BCols = B.NumCols(),
2506  CRows = C.NumRows(), CCols = C.NumCols();
2507  if (transA == kTrans) std::swap(ARows, ACols);
2508  if (transB == kTrans) std::swap(BRows, BCols);
2509  if (transC == kTrans) std::swap(CRows, CCols);
2510  KALDI_ASSERT( CCols == ARows && ACols == BRows && BCols == CRows && "TraceMatMatMat: args have mismatched dimensions.");
2511  if (ARows*BCols < std::min(BRows*CCols, CRows*ACols)) {
2512  Matrix<Real> AB(ARows, BCols);
2513  AB.AddMatMat(1.0, A, transA, B, transB, 0.0); // AB = A * B.
2514  return TraceMatMat(AB, C, transC);
2515  } else if ( BRows*CCols < CRows*ACols) {
2516  Matrix<Real> BC(BRows, CCols);
2517  BC.AddMatMat(1.0, B, transB, C, transC, 0.0); // BC = B * C.
2518  return TraceMatMat(BC, A, transA);
2519  } else {
2520  Matrix<Real> CA(CRows, ACols);
2521  CA.AddMatMat(1.0, C, transC, A, transA, 0.0); // CA = C * A
2522  return TraceMatMat(CA, B, transB);
2523  }
2524 }
2525 
2526 template
2528  const MatrixBase<float> &B, MatrixTransposeType transB,
2529  const MatrixBase<float> &C, MatrixTransposeType transC);
2530 
2531 template
2532 double TraceMatMatMat(const MatrixBase<double> &A, MatrixTransposeType transA,
2533  const MatrixBase<double> &B, MatrixTransposeType transB,
2534  const MatrixBase<double> &C, MatrixTransposeType transC);
2535 
2536 
2537 template <typename Real>
2539  const MatrixBase<Real> &B, MatrixTransposeType transB,
2540  const MatrixBase<Real> &C, MatrixTransposeType transC,
2541  const MatrixBase<Real> &D, MatrixTransposeType transD) {
2542  MatrixIndexT ARows = A.NumRows(), ACols = A.NumCols(), BRows = B.NumRows(), BCols = B.NumCols(),
2543  CRows = C.NumRows(), CCols = C.NumCols(), DRows = D.NumRows(), DCols = D.NumCols();
2544  if (transA == kTrans) std::swap(ARows, ACols);
2545  if (transB == kTrans) std::swap(BRows, BCols);
2546  if (transC == kTrans) std::swap(CRows, CCols);
2547  if (transD == kTrans) std::swap(DRows, DCols);
2548  KALDI_ASSERT( DCols == ARows && ACols == BRows && BCols == CRows && CCols == DRows && "TraceMatMatMat: args have mismatched dimensions.");
2549  if (ARows*BCols < std::min(BRows*CCols, std::min(CRows*DCols, DRows*ACols))) {
2550  Matrix<Real> AB(ARows, BCols);
2551  AB.AddMatMat(1.0, A, transA, B, transB, 0.0); // AB = A * B.
2552  return TraceMatMatMat(AB, kNoTrans, C, transC, D, transD);
2553  } else if ((BRows*CCols) < std::min(CRows*DCols, DRows*ACols)) {
2554  Matrix<Real> BC(BRows, CCols);
2555  BC.AddMatMat(1.0, B, transB, C, transC, 0.0); // BC = B * C.
2556  return TraceMatMatMat(BC, kNoTrans, D, transD, A, transA);
2557  } else if (CRows*DCols < DRows*ACols) {
2558  Matrix<Real> CD(CRows, DCols);
2559  CD.AddMatMat(1.0, C, transC, D, transD, 0.0); // CD = C * D
2560  return TraceMatMatMat(CD, kNoTrans, A, transA, B, transB);
2561  } else {
2562  Matrix<Real> DA(DRows, ACols);
2563  DA.AddMatMat(1.0, D, transD, A, transA, 0.0); // DA = D * A
2564  return TraceMatMatMat(DA, kNoTrans, B, transB, C, transC);
2565  }
2566 }
2567 
2568 template
2570  const MatrixBase<float> &B, MatrixTransposeType transB,
2571  const MatrixBase<float> &C, MatrixTransposeType transC,
2572  const MatrixBase<float> &D, MatrixTransposeType transD);
2573 
2574 template
2576  const MatrixBase<double> &B, MatrixTransposeType transB,
2577  const MatrixBase<double> &C, MatrixTransposeType transC,
2578  const MatrixBase<double> &D, MatrixTransposeType transD);
2579 
2580 template<typename Real> void SortSvd(VectorBase<Real> *s, MatrixBase<Real> *U,
2581  MatrixBase<Real> *Vt, bool sort_on_absolute_value) {
2583  MatrixIndexT num_singval = s->Dim();
2584  KALDI_ASSERT(U == NULL || U->NumCols() == num_singval);
2585  KALDI_ASSERT(Vt == NULL || Vt->NumRows() == num_singval);
2586 
2587  std::vector<std::pair<Real, MatrixIndexT> > vec(num_singval);
2588  // negative because we want revese order.
2589  for (MatrixIndexT d = 0; d < num_singval; d++) {
2590  Real val = (*s)(d),
2591  sort_val = -(sort_on_absolute_value ? std::abs(val) : val);
2592  vec[d] = std::pair<Real, MatrixIndexT>(sort_val, d);
2593  }
2594  std::sort(vec.begin(), vec.end());
2595  Vector<Real> s_copy(*s);
2596  for (MatrixIndexT d = 0; d < num_singval; d++)
2597  (*s)(d) = s_copy(vec[d].second);
2598  if (U != NULL) {
2599  Matrix<Real> Utmp(*U);
2600  MatrixIndexT dim = Utmp.NumRows();
2601  for (MatrixIndexT d = 0; d < num_singval; d++) {
2602  MatrixIndexT oldidx = vec[d].second;
2603  for (MatrixIndexT e = 0; e < dim; e++)
2604  (*U)(e, d) = Utmp(e, oldidx);
2605  }
2606  }
2607  if (Vt != NULL) {
2608  Matrix<Real> Vttmp(*Vt);
2609  for (MatrixIndexT d = 0; d < num_singval; d++)
2610  (*Vt).Row(d).CopyFromVec(Vttmp.Row(vec[d].second));
2611  }
2612 }
2613 
2614 template
2616  MatrixBase<float> *Vt, bool);
2617 
2618 template
2620  MatrixBase<double> *Vt, bool);
2621 
2622 template<typename Real>
2624  MatrixBase<Real> *D) {
2625  MatrixIndexT n = re.Dim();
2626  KALDI_ASSERT(im.Dim() == n && D->NumRows() == n && D->NumCols() == n);
2627 
2628  MatrixIndexT j = 0;
2629  D->SetZero();
2630  while (j < n) {
2631  if (im(j) == 0) { // Real eigenvalue
2632  (*D)(j, j) = re(j);
2633  j++;
2634  } else { // First of a complex pair
2635  KALDI_ASSERT(j+1 < n && ApproxEqual(im(j+1), -im(j))
2636  && ApproxEqual(re(j+1), re(j)));
2638  Real lambda = re(j), mu = im(j);
2639  // create 2x2 block [lambda, mu; -mu, lambda]
2640  (*D)(j, j) = lambda;
2641  (*D)(j, j+1) = mu;
2642  (*D)(j+1, j) = -mu;
2643  (*D)(j+1, j+1) = lambda;
2644  j += 2;
2645  }
2646  }
2647 }
2648 
2649 template
2651  MatrixBase<float> *D);
2652 template
2654  MatrixBase<double> *D);
2655 
2656 
2657 
2658 template<typename Real>
2659 bool AttemptComplexPower(Real *x_re, Real *x_im, Real power) {
2660  // Used in Matrix<Real>::Power().
2661  // Attempts to take the complex value x to the power "power",
2662  // assuming that power is fractional (i.e. we don't treat integers as a
2663  // special case). Returns false if this is not possible, either
2664  // because x is negative and real (hence there is no obvious answer
2665  // that is "closest to 1", and anyway this case does not make sense
2666  // in the Matrix<Real>::Power() routine);
2667  // or because power is negative, and x is zero.
2668 
2669  // First solve for r and theta in
2670  // x_re = r*cos(theta), x_im = r*sin(theta)
2671  if (*x_re < 0.0 && *x_im == 0.0) return false; // can't do
2672  // it for negative real values.
2673  Real r = std::sqrt((*x_re * *x_re) + (*x_im * *x_im)); // r == radius.
2674  if (power < 0.0 && r == 0.0) return false;
2675  Real theta = std::atan2(*x_im, *x_re);
2676  // Take the power.
2677  r = std::pow(r, power);
2678  theta *= power;
2679  *x_re = r * std::cos(theta);
2680  *x_im = r * std::sin(theta);
2681  return true;
2682 }
2683 
2684 template
2685 bool AttemptComplexPower(float *x_re, float *x_im, float power);
2686 template
2687 bool AttemptComplexPower(double *x_re, double *x_im, double power);
2688 
2689 
2690 
2691 template <typename Real>
2693  const MatrixBase<Real> &B,
2694  MatrixTransposeType trans) { // tr(A B), equivalent to sum of each element of A times same element in B'
2695  MatrixIndexT aStride = A.stride_, bStride = B.stride_;
2696  if (trans == kNoTrans) {
2697  KALDI_ASSERT(A.NumRows() == B.NumCols() && A.NumCols() == B.NumRows());
2698  Real ans = 0.0;
2699  Real *adata = A.data_, *bdata = B.data_;
2700  MatrixIndexT arows = A.NumRows(), acols = A.NumCols();
2701  for (MatrixIndexT row = 0;row < arows;row++, adata+=aStride, bdata++)
2702  ans += cblas_Xdot(acols, adata, 1, bdata, bStride);
2703  return ans;
2704  } else {
2705  KALDI_ASSERT(A.NumRows() == B.NumRows() && A.NumCols() == B.NumCols());
2706  Real ans = 0.0;
2707  Real *adata = A.data_, *bdata = B.data_;
2708  MatrixIndexT arows = A.NumRows(), acols = A.NumCols();
2709  for (MatrixIndexT row = 0;row < arows;row++, adata+=aStride, bdata+=bStride)
2710  ans += cblas_Xdot(acols, adata, 1, bdata, 1);
2711  return ans;
2712  }
2713 }
2714 
2715 
2716 // Instantiate the template above for float and double.
2717 template
2718 float TraceMatMat(const MatrixBase<float> &A,
2719  const MatrixBase<float> &B,
2720  MatrixTransposeType trans);
2721 template
2722 double TraceMatMat(const MatrixBase<double> &A,
2723  const MatrixBase<double> &B,
2724  MatrixTransposeType trans);
2725 
2726 
2727 template<typename Real>
2728 Real MatrixBase<Real>::LogSumExp(Real prune) const {
2729  Real sum;
2730  if (sizeof(sum) == 8) sum = kLogZeroDouble;
2731  else sum = kLogZeroFloat;
2732  Real max_elem = Max(), cutoff;
2733  if (sizeof(Real) == 4) cutoff = max_elem + kMinLogDiffFloat;
2734  else cutoff = max_elem + kMinLogDiffDouble;
2735  if (prune > 0.0 && max_elem - prune > cutoff) // explicit pruning...
2736  cutoff = max_elem - prune;
2737 
2738  double sum_relto_max_elem = 0.0;
2739 
2740  for (MatrixIndexT i = 0; i < num_rows_; i++) {
2741  for (MatrixIndexT j = 0; j < num_cols_; j++) {
2742  BaseFloat f = (*this)(i, j);
2743  if (f >= cutoff)
2744  sum_relto_max_elem += kaldi::Exp(f - max_elem);
2745  }
2746  }
2747  return max_elem + kaldi::Log(sum_relto_max_elem);
2748 }
2749 
2750 template<typename Real>
2752  Real max = this->Max(), sum = 0.0;
2753  // the 'max' helps to get in good numeric range.
2754  for (MatrixIndexT i = 0; i < num_rows_; i++)
2755  for (MatrixIndexT j = 0; j < num_cols_; j++)
2756  sum += ((*this)(i, j) = kaldi::Exp((*this)(i, j) - max));
2757  this->Scale(1.0 / sum);
2758  return max + kaldi::Log(sum);
2759 }
2760 
2761 template<typename Real>
2763  KALDI_ASSERT(SameDim(*this, src));
2764 
2765  if (num_cols_ == stride_ && src.num_cols_ == src.stride_) {
2766  SubVector<Real> src_vec(src.data_, num_rows_ * num_cols_),
2767  dst_vec(this->data_, num_rows_ * num_cols_);
2768  dst_vec.Tanh(src_vec);
2769  } else {
2770  for (MatrixIndexT r = 0; r < num_rows_; r++) {
2771  SubVector<Real> src_vec(src, r), dest_vec(*this, r);
2772  dest_vec.Tanh(src_vec);
2773  }
2774  }
2775 }
2776 
2777 template<typename Real>
2779  KALDI_ASSERT(SameDim(*this, src));
2780  int32 num_rows = num_rows_, num_cols = num_cols_;
2781  for (MatrixIndexT r = 0; r < num_rows; r++) {
2782  Real *row_data = this->RowData(r);
2783  const Real *src_row_data = src.RowData(r);
2784  for (MatrixIndexT c = 0; c < num_cols; c++) {
2785  Real x = src_row_data[c], y;
2786  if (x > 10.0) y = x; // avoid exponentiating large numbers; function
2787  // approaches y=x.
2788  else y = Log1p(kaldi::Exp(x)); // these defined in kaldi-math.h
2789  row_data[c] = y;
2790  }
2791  }
2792 }
2793 
2794 template<typename Real>
2795 void MatrixBase<Real>::GroupPnorm(const MatrixBase<Real> &src, Real power) {
2796  KALDI_ASSERT(src.NumCols() % this->NumCols() == 0 &&
2797  src.NumRows() == this->NumRows());
2798  int group_size = src.NumCols() / this->NumCols(),
2799  num_rows = this->NumRows(), num_cols = this->NumCols();
2800  for (MatrixIndexT i = 0; i < num_rows; i++)
2801  for (MatrixIndexT j = 0; j < num_cols; j++)
2802  (*this)(i, j) = src.Row(i).Range(j * group_size, group_size).Norm(power);
2803 }
2804 
2805 template<typename Real>
2807  KALDI_ASSERT(src.NumCols() % this->NumCols() == 0 &&
2808  src.NumRows() == this->NumRows());
2809  int group_size = src.NumCols() / this->NumCols(),
2810  num_rows = this->NumRows(), num_cols = this->NumCols();
2811  for (MatrixIndexT i = 0; i < num_rows; i++) {
2812  const Real *src_row_data = src.RowData(i);
2813  for (MatrixIndexT j = 0; j < num_cols; j++) {
2814  Real max_val = -1e20;
2815  for (MatrixIndexT k = 0; k < group_size; k++) {
2816  Real src_data = src_row_data[j * group_size + k];
2817  if (src_data > max_val)
2818  max_val = src_data;
2819  }
2820  (*this)(i, j) = max_val;
2821  }
2822  }
2823 }
2824 
2825 template<typename Real>
2827  const MatrixIndexT *indices) {
2828  KALDI_ASSERT(NumRows() == src.NumRows());
2829  MatrixIndexT num_rows = num_rows_, num_cols = num_cols_,
2830  this_stride = stride_, src_stride = src.stride_;
2831  Real *this_data = this->data_;
2832  const Real *src_data = src.data_;
2833 #ifdef KALDI_PARANOID
2834  MatrixIndexT src_cols = src.NumCols();
2835  for (MatrixIndexT i = 0; i < num_cols; i++)
2836  KALDI_ASSERT(indices[i] >= -1 && indices[i] < src_cols);
2837 #endif
2838 
2839  // For the sake of memory locality we do this row by row, rather
2840  // than doing it column-wise using cublas_Xcopy
2841  for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride, src_data += src_stride) {
2842  const MatrixIndexT *index_ptr = &(indices[0]);
2843  for (MatrixIndexT c = 0; c < num_cols; c++, index_ptr++) {
2844  if (*index_ptr < 0) this_data[c] = 0;
2845  else this_data[c] = src_data[*index_ptr];
2846  }
2847  }
2848 }
2849 
2850 
2851 template<typename Real>
2853  const MatrixIndexT *indices) {
2854  KALDI_ASSERT(NumRows() == src.NumRows());
2855  MatrixIndexT num_rows = num_rows_, num_cols = num_cols_,
2856  this_stride = stride_, src_stride = src.stride_;
2857  Real *this_data = this->data_;
2858  const Real *src_data = src.data_;
2859 #ifdef KALDI_PARANOID
2860  MatrixIndexT src_cols = src.NumCols();
2861  for (MatrixIndexT i = 0; i < num_cols; i++)
2862  KALDI_ASSERT(indices[i] >= -1 && indices[i] < src_cols);
2863 #endif
2864 
2865  // For the sake of memory locality we do this row by row, rather
2866  // than doing it column-wise using cublas_Xcopy
2867  for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride, src_data += src_stride) {
2868  const MatrixIndexT *index_ptr = &(indices[0]);
2869  for (MatrixIndexT c = 0; c < num_cols; c++, index_ptr++) {
2870  if (*index_ptr >= 0)
2871  this_data[c] += src_data[*index_ptr];
2872  }
2873  }
2874 }
2875 
2876 template<typename Real>
2878  const MatrixIndexT *indices) {
2879  KALDI_ASSERT(NumCols() == src.NumCols());
2880  MatrixIndexT num_rows = num_rows_, num_cols = num_cols_,
2881  this_stride = stride_;
2882  Real *this_data = this->data_;
2883 
2884  for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride) {
2885  MatrixIndexT index = indices[r];
2886  if (index < 0) memset(this_data, 0, sizeof(Real) * num_cols_);
2887  else cblas_Xcopy(num_cols, src.RowData(index), 1, this_data, 1);
2888  }
2889 }
2890 
2891 template<typename Real>
2892 void MatrixBase<Real>::CopyRows(const Real *const *src) {
2893  MatrixIndexT num_rows = num_rows_,
2894  num_cols = num_cols_, this_stride = stride_;
2895  Real *this_data = this->data_;
2896 
2897  for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride) {
2898  const Real *const src_data = src[r];
2899  if (src_data == NULL) memset(this_data, 0, sizeof(Real) * num_cols);
2900  else cblas_Xcopy(num_cols, src_data, 1, this_data, 1);
2901  }
2902 }
2903 
2904 template<typename Real>
2905 void MatrixBase<Real>::CopyToRows(Real *const *dst) const {
2906  MatrixIndexT num_rows = num_rows_,
2907  num_cols = num_cols_, this_stride = stride_;
2908  const Real *this_data = this->data_;
2909 
2910  for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride) {
2911  Real *const dst_data = dst[r];
2912  if (dst_data != NULL)
2913  cblas_Xcopy(num_cols, this_data, 1, dst_data, 1);
2914  }
2915 }
2916 
2917 template<typename Real>
2919  const MatrixBase<Real> &src,
2920  const MatrixIndexT *indexes) {
2921  KALDI_ASSERT(NumCols() == src.NumCols());
2922  MatrixIndexT num_rows = num_rows_,
2923  num_cols = num_cols_, this_stride = stride_;
2924  Real *this_data = this->data_;
2925 
2926  for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride) {
2927  MatrixIndexT index = indexes[r];
2928  KALDI_ASSERT(index >= -1 && index < src.NumRows());
2929  if (index != -1)
2930  cblas_Xaxpy(num_cols, alpha, src.RowData(index), 1, this_data, 1);
2931  }
2932 }
2933 
2934 template<typename Real>
2935 void MatrixBase<Real>::AddRows(Real alpha, const Real *const *src) {
2936  MatrixIndexT num_rows = num_rows_,
2937  num_cols = num_cols_, this_stride = stride_;
2938  Real *this_data = this->data_;
2939 
2940  for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride) {
2941  const Real *const src_data = src[r];
2942  if (src_data != NULL)
2943  cblas_Xaxpy(num_cols, alpha, src_data, 1, this_data, 1);
2944  }
2945 }
2946 
2947 template<typename Real>
2949  const MatrixIndexT *indexes,
2950  MatrixBase<Real> *dst) const {
2951  KALDI_ASSERT(NumCols() == dst->NumCols());
2952  MatrixIndexT num_rows = num_rows_,
2953  num_cols = num_cols_, this_stride = stride_;
2954  Real *this_data = this->data_;
2955 
2956  for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride) {
2957  MatrixIndexT index = indexes[r];
2958  KALDI_ASSERT(index >= -1 && index < dst->NumRows());
2959  if (index != -1)
2960  cblas_Xaxpy(num_cols, alpha, this_data, 1, dst->RowData(index), 1);
2961  }
2962 }
2963 
2964 template<typename Real>
2965 void MatrixBase<Real>::AddToRows(Real alpha, Real *const *dst) const {
2966  MatrixIndexT num_rows = num_rows_,
2967  num_cols = num_cols_, this_stride = stride_;
2968  const Real *this_data = this->data_;
2969 
2970  for (MatrixIndexT r = 0; r < num_rows; r++, this_data += this_stride) {
2971  Real *const dst_data = dst[r];
2972  if (dst_data != NULL)
2973  cblas_Xaxpy(num_cols, alpha, this_data, 1, dst_data, 1);
2974  }
2975 }
2976 
2977 template<typename Real>
2979  KALDI_ASSERT(SameDim(*this, src));
2980 
2981  if (num_cols_ == stride_ && src.num_cols_ == src.stride_) {
2982  SubVector<Real> src_vec(src.data_, num_rows_ * num_cols_),
2983  dst_vec(this->data_, num_rows_ * num_cols_);
2984  dst_vec.Sigmoid(src_vec);
2985  } else {
2986  for (MatrixIndexT r = 0; r < num_rows_; r++) {
2987  SubVector<Real> src_vec(src, r), dest_vec(*this, r);
2988  dest_vec.Sigmoid(src_vec);
2989  }
2990  }
2991 }
2992 
2993 template<typename Real>
2995  const MatrixBase<Real> &diff) {
2996  KALDI_ASSERT(SameDim(*this, value) && SameDim(*this, diff));
2997  MatrixIndexT num_rows = num_rows_, num_cols = num_cols_,
2998  stride = stride_, value_stride = value.stride_, diff_stride = diff.stride_;
2999  Real *data = data_;
3000  const Real *value_data = value.data_, *diff_data = diff.data_;
3001  for (MatrixIndexT r = 0; r < num_rows; r++) {
3002  for (MatrixIndexT c = 0; c < num_cols; c++)
3003  data[c] = diff_data[c] * value_data[c] * (1.0 - value_data[c]);
3004  data += stride;
3005  value_data += value_stride;
3006  diff_data += diff_stride;
3007  }
3008 }
3009 
3010 template<typename Real>
3012  const MatrixBase<Real> &diff) {
3013  KALDI_ASSERT(SameDim(*this, value) && SameDim(*this, diff));
3014  MatrixIndexT num_rows = num_rows_, num_cols = num_cols_,
3015  stride = stride_, value_stride = value.stride_, diff_stride = diff.stride_;
3016  Real *data = data_;
3017  const Real *value_data = value.data_, *diff_data = diff.data_;
3018  for (MatrixIndexT r = 0; r < num_rows; r++) {
3019  for (MatrixIndexT c = 0; c < num_cols; c++)
3020  data[c] = diff_data[c] * (1.0 - (value_data[c] * value_data[c]));
3021  data += stride;
3022  value_data += value_stride;
3023  diff_data += diff_stride;
3024  }
3025 }
3026 
3027 
3028 template<typename Real>
3029 template<typename OtherReal>
3030 void MatrixBase<Real>::AddVecToRows(const Real alpha, const VectorBase<OtherReal> &v) {
3031  const MatrixIndexT num_rows = num_rows_, num_cols = num_cols_,
3032  stride = stride_;
3033  KALDI_ASSERT(v.Dim() == num_cols);
3034  if(num_cols <= 64) {
3035  Real *data = data_;
3036  const OtherReal *vdata = v.Data();
3037  for (MatrixIndexT i = 0; i < num_rows; i++, data += stride) {
3038  for (MatrixIndexT j = 0; j < num_cols; j++)
3039  data[j] += alpha * vdata[j];
3040  }
3041 
3042  } else {
3043  Vector<OtherReal> ones(num_rows);
3044  ones.Set(1.0);
3045  this->AddVecVec(alpha, ones, v);
3046  }
3047 }
3048 
3049 template void MatrixBase<float>::AddVecToRows(const float alpha,
3050  const VectorBase<float> &v);
3051 template void MatrixBase<float>::AddVecToRows(const float alpha,
3052  const VectorBase<double> &v);
3053 template void MatrixBase<double>::AddVecToRows(const double alpha,
3054  const VectorBase<float> &v);
3055 template void MatrixBase<double>::AddVecToRows(const double alpha,
3056  const VectorBase<double> &v);
3057 
3058 
3059 template<typename Real>
3060 template<typename OtherReal>
3061 void MatrixBase<Real>::AddVecToCols(const Real alpha, const VectorBase<OtherReal> &v) {
3062  const MatrixIndexT num_rows = num_rows_, num_cols = num_cols_,
3063  stride = stride_;
3064  KALDI_ASSERT(v.Dim() == num_rows);
3065 
3066  if (num_rows <= 64) {
3067  Real *data = data_;
3068  const OtherReal *vdata = v.Data();
3069  for (MatrixIndexT i = 0; i < num_rows; i++, data += stride) {
3070  Real to_add = alpha * vdata[i];
3071  for (MatrixIndexT j = 0; j < num_cols; j++)
3072  data[j] += to_add;
3073  }
3074 
3075  } else {
3076  Vector<OtherReal> ones(num_cols);
3077  ones.Set(1.0);
3078  this->AddVecVec(alpha, v, ones);
3079  }
3080 }
3081 
3082 template void MatrixBase<float>::AddVecToCols(const float alpha,
3083  const VectorBase<float> &v);
3084 template void MatrixBase<float>::AddVecToCols(const float alpha,
3085  const VectorBase<double> &v);
3086 template void MatrixBase<double>::AddVecToCols(const double alpha,
3087  const VectorBase<float> &v);
3088 template void MatrixBase<double>::AddVecToCols(const double alpha,
3089  const VectorBase<double> &v);
3090 
3091 //Explicit instantiation of the classes
3092 //Apparently, it seems to be necessary that the instantiation
3093 //happens at the end of the file. Otherwise, not all the member
3094 //functions will get instantiated.
3095 
3096 template class Matrix<float>;
3097 template class Matrix<double>;
3098 template class MatrixBase<float>;
3099 template class MatrixBase<double>;
3100 template class SubMatrix<float>;
3101 template class SubMatrix<double>;
3102 
3103 } // namespace kaldi
void AddMat2(const Real alpha, const MatrixBase< Real > &M, MatrixTransposeType transM, const Real beta)
rank-N update: if (transM == kNoTrans) (*this) = beta*(*this) + alpha * M * M^T, or (if transM == kTr...
Definition: sp-matrix.cc:1110
This code computes Goodness of Pronunciation (GOP) and extracts phone-level pronunciation feature for...
Definition: chain.dox:20
double Exp(double x)
Definition: kaldi-math.h:83
void Add(const Real alpha)
Add a scalar to each element.
Matrix()
Empty constructor.
void InvertDouble(Real *LogDet=NULL, Real *det_sign=NULL, bool inverse_needed=true)
matrix inverse [double].
Real Min() const
Returns minimum element of matrix.
This class provides a way for switching between double and float types.
Definition: matrix-common.h:84
#define KALDI_SWAP2(a)
Definition: kaldi-utils.h:114
template double TraceMatMat(const MatrixBase< double > &A, const MatrixBase< double > &B, MatrixTransposeType trans)
Packed symetric matrix class.
Definition: matrix-common.h:62
Real LargestAbsElem() const
largest absolute value.
int16 mSampleSize
Sample size.
Definition: kaldi-matrix.h:961
template void CreateEigenvalueMatrix(const VectorBase< double > &re, const VectorBase< double > &im, MatrixBase< double > *D)
void cblas_Xsyrk(const MatrixTransposeType trans, const MatrixIndexT dim_c, const MatrixIndexT other_dim_a, const float alpha, const float *A, const MatrixIndexT a_stride, const float beta, float *C, const MatrixIndexT c_stride)
bool IsDiagonal(Real cutoff=1.0e-05) const
Returns true if matrix is Diagonal.
MatrixResizeType
Definition: matrix-common.h:37
OfflineFeatureTpl< MfccComputer > Mfcc
Definition: feature-mfcc.h:147
void clapack_Xgetri2(KaldiBlasInt *num_rows, float *Mdata, KaldiBlasInt *stride, KaldiBlasInt *pivot, float *p_work, KaldiBlasInt *l_work, KaldiBlasInt *result)
void Sigmoid(const VectorBase< Real > &src)
Sets each element of *this to the sigmoid function of the corresponding element of "src"...
MatrixIndexT stride_
< Number of rows
Definition: kaldi-matrix.h:816
float RandUniform(struct RandomState *state=NULL)
Returns a random number strictly between 0 and 1.
Definition: kaldi-math.h:151
MatrixIndexT NumCols() const
Returns number of columns (or zero for empty matrix).
Definition: kaldi-matrix.h:67
Real Cond() const
Returns condition number by computing Svd.
Base class which provides matrix operations not involving resizing or allocation. ...
Definition: kaldi-matrix.h:49
const Real * Data() const
Gives pointer to raw data (const).
Definition: kaldi-matrix.h:79
void ReadBasicType(std::istream &is, bool binary, T *t)
ReadBasicType is the name of the read function for bool, integer types, and floating-point types...
Definition: io-funcs-inl.h:55
const std::pair< MatrixIndexT, Real > & GetElement(MatrixIndexT i) const
get an indexed element (0 <= i < NumElements()).
Definition: sparse-matrix.h:77
void Tanh(const VectorBase< Real > &src)
Sets each element of *this to the tanh of the corresponding element of "src".
Real * data_
data memory area
Definition: kaldi-matrix.h:808
OfflineFeatureTpl< PlpComputer > Plp
Definition: feature-plp.h:169
Real Max() const
Returns maximum element of matrix.
bool ApproxEqual(const MatrixBase< Real > &other, float tol=0.01) const
Returns true if ((*this)-other).FrobeniusNorm() <= tol * (*this).FrobeniusNorm(). ...
void Xgemv_sparsevec(MatrixTransposeType trans, MatrixIndexT num_rows, MatrixIndexT num_cols, Real alpha, const Real *Mdata, MatrixIndexT stride, const Real *xdata, MatrixIndexT incX, Real beta, Real *ydata, MatrixIndexT incY)
void SymPosSemiDefEig(VectorBase< Real > *s, MatrixBase< Real > *P, Real check_thresh=0.001)
Uses Svd to compute the eigenvalue decomposition of a symmetric positive semi-definite matrix: (*this...
Real Trace(bool check_square=true) const
Returns trace of matrix.
Real * RowData(MatrixIndexT i)
Returns pointer to data for one row (non-const)
Definition: kaldi-matrix.h:87
void clapack_Xgetrf2(KaldiBlasInt *num_rows, KaldiBlasInt *num_cols, float *Mdata, KaldiBlasInt *stride, KaldiBlasInt *pivot, KaldiBlasInt *result)
void AddMat(const Real alpha, const MatrixBase< Real > &M, MatrixTransposeType transA=kNoTrans)
*this += alpha * M [or M^T]
void swap(basic_filebuf< CharT, Traits > &x, basic_filebuf< CharT, Traits > &y)
void LapackGesvd(VectorBase< Real > *s, MatrixBase< Real > *U, MatrixBase< Real > *Vt)
float RandGauss(struct RandomState *state=NULL)
Definition: kaldi-math.h:155
kaldi::int32 int32
void AddToDiag(const Real alpha)
Add a scalar to each diagonal element.
void ReadToken(std::istream &is, bool binary, std::string *str)
ReadToken gets the next token and puts it in str (exception on failure).
Definition: io-funcs.cc:154
Real * Data_workaround() const
A workaround that allows SubMatrix to get a pointer to non-const data for const Matrix.
Definition: kaldi-matrix.h:803
A class for storing matrices.
Definition: kaldi-matrix.h:823
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)
void cblas_Xcopy(const int N, const float *X, const int incX, float *Y, const int incY)
void Swap(Matrix< Real > *other)
Swaps the contents of *this and *other. Shallow swap.
void CopyFromMat(const MatrixBase< OtherReal > &M, MatrixTransposeType trans=kNoTrans)
Copy given matrix. (no resize is done).
MatrixIndexT NumRows() const
void clapack_Xgesvd(char *v, char *u, KaldiBlasInt *num_cols, KaldiBlasInt *num_rows, float *Mdata, KaldiBlasInt *stride, float *sv, float *Vdata, KaldiBlasInt *vstride, float *Udata, KaldiBlasInt *ustride, float *p_work, KaldiBlasInt *l_work, KaldiBlasInt *result)
void Read(std::istream &is, bool binary)
uint32 UnsignedMatrixIndexT
void mul_elements(const MatrixIndexT dim, const double *a, double *b)
This is not really a wrapper for CBLAS as CBLAS does not have this; in future we could extend this so...
int Peek(std::istream &is, bool binary)
Peek consumes whitespace (if binary == false) and then returns the peek() value of the stream...
Definition: io-funcs.cc:145
uint64 data_
template bool ReadHtk(std::istream &is, Matrix< double > *M, HtkHeader *header_ptr)
int32 mSamplePeriod
Sample period.
Definition: kaldi-matrix.h:959
template double TraceMatMatMat(const MatrixBase< double > &A, MatrixTransposeType transA, const MatrixBase< double > &B, MatrixTransposeType transB, const MatrixBase< double > &C, MatrixTransposeType transC)
void DestructiveSvd(VectorBase< Real > *s, MatrixBase< Real > *U, MatrixBase< Real > *Vt)
Singular value decomposition Major limitations: For nonsquare matrices, we assume m>=n (NumRows >= Nu...
bool SameDim(const MatrixBase< Real > &M, const MatrixBase< Real > &N)
void TestUninitialized() const
float cblas_Xdot(const int N, const float *const X, const int incX, const float *const Y, const int incY)
bool IsZero(Real cutoff=1.0e-05) const
Returns true if matrix is all zeros.
void GetImagEigenvalues(VectorBase< Real > *i_out)
Definition: jama-eig.h:61
MatrixIndexT NumCols() const
void Read(std::istream &in, bool binary, bool add=false)
read from stream.
bool IsSymmetric(Real cutoff=1.0e-05) const
Returns true if matrix is Symmetric.
MatrixIndexT Stride() const
Stride (distance in memory between each row). Will be >= NumCols.
Definition: kaldi-matrix.h:70
int32 MatrixIndexT
Definition: matrix-common.h:98
void cblas_Xscal(const int N, const float alpha, float *data, const int inc)
const SubVector< Real > Row(MatrixIndexT i) const
Return specific row of matrix [const].
Definition: kaldi-matrix.h:188
double Log(double x)
Definition: kaldi-math.h:100
MatrixIndexT num_rows_
< Number of columns
Definition: kaldi-matrix.h:813
#define KALDI_MEMALIGN(align, size, pp_orig)
Definition: kaldi-utils.h:58
static const double kMinLogDiffDouble
Definition: kaldi-math.h:124
void RandGauss2(float *a, float *b, RandomState *state)
Definition: kaldi-math.cc:139
struct rnnlm::@11::@12 n
template bool WriteHtk(std::ostream &os, const MatrixBase< double > &M, HtkHeader htk_hdr)
MatrixStrideType
Definition: matrix-common.h:44
OfflineFeatureTpl< FbankComputer > Fbank
void GetV(MatrixBase< Real > *V_out)
Definition: jama-eig.h:47
void Transpose()
Transpose the matrix.
void cblas_Xsymm(const float alpha, MatrixIndexT sz, const float *Adata, MatrixIndexT a_stride, const float *Bdata, MatrixIndexT b_stride, const float beta, float *Mdata, MatrixIndexT stride)
void AddMatMat(const Real alpha, const MatrixBase< Real > &A, MatrixTransposeType transA, const MatrixBase< Real > &B, MatrixTransposeType transB, const Real beta)
int MachineIsLittleEndian()
Definition: kaldi-utils.h:83
#define KALDI_ERR
Definition: kaldi-error.h:147
#define KALDI_MEMALIGN_FREE(x)
Definition: kaldi-utils.h:60
uint16 mSampleKind
Sample kind.
Definition: kaldi-matrix.h:963
void AddMatMatMat(const Real alpha, const MatrixBase< Real > &A, MatrixTransposeType transA, const MatrixBase< Real > &B, MatrixTransposeType transB, const MatrixBase< Real > &C, MatrixTransposeType transC, const Real beta)
this <– beta*this + alpha*A*B*C.
Packed symetric matrix class.
Definition: matrix-common.h:63
MatrixIndexT num_cols_
these attributes store the real matrix size as it is stored in memory including memalignment ...
Definition: kaldi-matrix.h:812
MatrixIndexT NumElements() const
Returns the number of nonzero elements.
Definition: sparse-matrix.h:74
#define KALDI_WARN
Definition: kaldi-error.h:150
Real TraceMatMat(const MatrixBase< Real > &A, const MatrixBase< Real > &B, MatrixTransposeType trans)
We need to declare this here as it will be a friend function.
Real * Data()
Returns a pointer to the start of the vector&#39;s data.
Definition: kaldi-vector.h:70
void WriteToken(std::ostream &os, bool binary, const char *token)
The WriteToken functions are for writing nonempty sequences of non-space characters.
Definition: io-funcs.cc:134
void cblas_Xgemm(const float alpha, MatrixTransposeType transA, const float *Adata, MatrixIndexT a_num_rows, MatrixIndexT a_num_cols, MatrixIndexT a_stride, MatrixTransposeType transB, const float *Bdata, MatrixIndexT b_stride, const float beta, float *Mdata, MatrixIndexT num_rows, MatrixIndexT num_cols, MatrixIndexT stride)
MatrixIndexT Dim() const
Returns the dimension of the vector.
Definition: kaldi-vector.h:64
void SetZero()
Sets matrix to zero.
void Scale(Real alpha)
Multiplies all elements by this constant.
template void SortSvd(VectorBase< double > *s, MatrixBase< double > *U, MatrixBase< double > *Vt, bool)
static const float kMinLogDiffFloat
Definition: kaldi-math.h:125
void MulColsVec(const VectorBase< Real > &scale)
Equivalent to (*this) = (*this) * diag(scale).
#define KALDI_ASSERT_IS_FLOATING_TYPE(F)
Definition: kaldi-utils.h:137
void GetRealEigenvalues(VectorBase< Real > *r_out)
Definition: jama-eig.h:55
A class representing a vector.
Definition: kaldi-vector.h:406
void cblas_Xaxpy(const int N, const float alpha, const float *X, const int incX, float *Y, const int incY)
MatrixIndexT NumRows() const
Returns number of rows (or zero for emtpy matrix).
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:185
MatrixIndexT NumRows() const
Returns number of rows (or zero for empty matrix).
Definition: kaldi-matrix.h:64
void Set(Real f)
Set all members of a vector to a specified value.
Real FrobeniusNorm() const
Frobenius norm, which is the sqrt of sum of square elements.
Real LogDet(Real *det_sign=NULL) const
Returns logdet of matrix.
void ApplyPow(Real power)
Take all elements of vector to a power.
Definition: kaldi-vector.h:179
MatrixIndexT NumRows() const
#define KALDI_STRCASECMP
Definition: kaldi-utils.h:147
void AddVecVec(const Real alpha, const VectorBase< OtherReal > &a, const VectorBase< OtherReal > &b)
*this += alpha * a * b^T
MatrixTransposeType
Definition: matrix-common.h:32
#define KALDI_VLOG(v)
Definition: kaldi-error.h:156
template bool WriteSphinx(std::ostream &os, const MatrixBase< double > &M)
SubMatrix< Real > Range(const MatrixIndexT row_offset, const MatrixIndexT num_rows, const MatrixIndexT col_offset, const MatrixIndexT num_cols) const
Return a sub-part of matrix.
Definition: kaldi-matrix.h:202
void WriteBasicType(std::ostream &os, bool binary, T t)
WriteBasicType is the name of the write function for bool, integer types, and floating-point types...
Definition: io-funcs-inl.h:34
bool Equal(const MatrixBase< Real > &other) const
Tests for exact equality. It&#39;s usually preferable to use ApproxEqual.
void Resize(const MatrixIndexT r, const MatrixIndexT c, MatrixResizeType resize_type=kSetZero, MatrixStrideType stride_type=kDefaultStride)
Sets matrix to a specified size (zero is OK as long as both r and c are zero).
const float kLogZeroFloat
Definition: kaldi-math.h:128
void Svd(VectorBase< Real > *s, MatrixBase< Real > *U, MatrixBase< Real > *Vt) const
Compute SVD (*this) = U diag(s) Vt.
double Log1p(double x)
Definition: kaldi-math.h:104
void OrthogonalizeRows()
This function orthogonalizes the rows of a matrix using the Gram-Schmidt process. ...
int32 mNSamples
Number of samples.
Definition: kaldi-matrix.h:957
void Invert(Real *log_det=NULL, Real *det_sign=NULL, bool inverse_needed=true)
matrix inverse.
Definition: kaldi-matrix.cc:38
MatrixIndexT NumCols() const
Returns number of columns (or zero for emtpy matrix).
bool IsUnit(Real cutoff=1.0e-05) const
Returns true if the matrix is all zeros, except for ones on diagonal.
void cblas_Xger(MatrixIndexT num_rows, MatrixIndexT num_cols, float alpha, const float *xdata, MatrixIndexT incX, const float *ydata, MatrixIndexT incY, float *Mdata, MatrixIndexT stride)
const double kLogZeroDouble
Definition: kaldi-math.h:129
Real VecVec(const VectorBase< Real > &a, const VectorBase< Real > &b)
Returns dot product between v1 and v2.
Definition: kaldi-vector.cc:37
void CopyToMat(MatrixBase< Real > *mat, MatrixTransposeType trans=kNoTrans) const
Copies contents to matrix.
A structure containing the HTK header.
Definition: kaldi-matrix.h:955
Sub-matrix representation.
Definition: kaldi-matrix.h:988
Represents a non-allocating general vector which can be defined as a sub-vector of higher-level vecto...
Definition: kaldi-vector.h:501
static bool ApproxEqual(float a, float b, float relative_tolerance=0.001)
return abs(a - b) <= relative_tolerance * (abs(a)+abs(b)).
Definition: kaldi-math.h:265
#define KALDI_SWAP4(a)
Definition: kaldi-utils.h:107
const SparseVector< Real > & Row(MatrixIndexT r) const
template bool AttemptComplexPower(double *x_re, double *x_im, double power)