kaldi-vector.cc
Go to the documentation of this file.
1 // matrix/kaldi-vector.cc
2 
3 // Copyright 2009-2011 Microsoft Corporation; Lukas Burget;
4 // Saarland University; Go Vivace Inc.; Ariya Rastrow;
5 // Petr Schwarz; Yanmin Qian; Jan Silovsky;
6 // Haihua Xu; Wei Shi
7 // 2015 Guoguo Chen
8 // 2017 Daniel Galvez
9 // 2019 Yiwen Shao
10 
11 // See ../../COPYING for clarification regarding multiple authors
12 //
13 // Licensed under the Apache License, Version 2.0 (the "License");
14 // you may not use this file except in compliance with the License.
15 // You may obtain a copy of the License at
16 //
17 // http://www.apache.org/licenses/LICENSE-2.0
18 //
19 // THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
20 // KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
21 // WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
22 // MERCHANTABLITY OR NON-INFRINGEMENT.
23 // See the Apache 2 License for the specific language governing permissions and
24 // limitations under the License.
25 
26 #include <algorithm>
27 #include <string>
28 #include "matrix/cblas-wrappers.h"
29 #include "matrix/kaldi-vector.h"
30 #include "matrix/kaldi-matrix.h"
31 #include "matrix/sp-matrix.h"
32 #include "matrix/sparse-matrix.h"
33 
34 namespace kaldi {
35 
36 template<typename Real>
37 Real VecVec(const VectorBase<Real> &a,
38  const VectorBase<Real> &b) {
39  MatrixIndexT adim = a.Dim();
40  KALDI_ASSERT(adim == b.Dim());
41  return cblas_Xdot(adim, a.Data(), 1, b.Data(), 1);
42 }
43 
44 template
45 float VecVec<>(const VectorBase<float> &a,
46  const VectorBase<float> &b);
47 template
48 double VecVec<>(const VectorBase<double> &a,
49  const VectorBase<double> &b);
50 
51 template<typename Real, typename OtherReal>
52 Real VecVec(const VectorBase<Real> &ra,
53  const VectorBase<OtherReal> &rb) {
54  MatrixIndexT adim = ra.Dim();
55  KALDI_ASSERT(adim == rb.Dim());
56  const Real *a_data = ra.Data();
57  const OtherReal *b_data = rb.Data();
58  Real sum = 0.0;
59  for (MatrixIndexT i = 0; i < adim; i++)
60  sum += a_data[i]*b_data[i];
61  return sum;
62 }
63 
64 // instantiate the template above.
65 template
66 float VecVec<>(const VectorBase<float> &ra,
67  const VectorBase<double> &rb);
68 template
69 double VecVec<>(const VectorBase<double> &ra,
70  const VectorBase<float> &rb);
71 
72 
73 template<>
74 template<>
75 void VectorBase<float>::AddVec(const float alpha,
76  const VectorBase<float> &v) {
77  KALDI_ASSERT(dim_ == v.dim_);
78  KALDI_ASSERT(&v != this);
79  cblas_Xaxpy(dim_, alpha, v.Data(), 1, data_, 1);
80 }
81 
82 template<>
83 template<>
84 void VectorBase<double>::AddVec(const double alpha,
85  const VectorBase<double> &v) {
86  KALDI_ASSERT(dim_ == v.dim_);
87  KALDI_ASSERT(&v != this);
88  cblas_Xaxpy(dim_, alpha, v.Data(), 1, data_, 1);
89 }
90 
91 template<typename Real>
92 void VectorBase<Real>::AddMatVec(const Real alpha,
93  const MatrixBase<Real> &M,
94  MatrixTransposeType trans,
95  const VectorBase<Real> &v,
96  const Real beta) {
97  KALDI_ASSERT((trans == kNoTrans && M.NumCols() == v.dim_ && M.NumRows() == dim_)
98  || (trans == kTrans && M.NumRows() == v.dim_ && M.NumCols() == dim_));
99  KALDI_ASSERT(&v != this);
100  cblas_Xgemv(trans, M.NumRows(), M.NumCols(), alpha, M.Data(), M.Stride(),
101  v.Data(), 1, beta, data_, 1);
102 }
103 
104 template<typename Real>
105 void VectorBase<Real>::AddMatSvec(const Real alpha,
106  const MatrixBase<Real> &M,
107  MatrixTransposeType trans,
108  const VectorBase<Real> &v,
109  const Real beta) {
110  KALDI_ASSERT((trans == kNoTrans && M.NumCols() == v.dim_ && M.NumRows() == dim_)
111  || (trans == kTrans && M.NumRows() == v.dim_ && M.NumCols() == dim_));
112  KALDI_ASSERT(&v != this);
113  Xgemv_sparsevec(trans, M.NumRows(), M.NumCols(), alpha, M.Data(), M.Stride(),
114  v.Data(), 1, beta, data_, 1);
115  return;
116  /*
117  MatrixIndexT this_dim = this->dim_, v_dim = v.dim_,
118  M_stride = M.Stride();
119  Real *this_data = this->data_;
120  const Real *M_data = M.Data(), *v_data = v.data_;
121  if (beta != 1.0) this->Scale(beta);
122  if (trans == kNoTrans) {
123  for (MatrixIndexT i = 0; i < v_dim; i++) {
124  Real v_i = v_data[i];
125  if (v_i == 0.0) continue;
126  // Add to *this, the i'th column of the Matrix, times v_i.
127  cblas_Xaxpy(this_dim, v_i * alpha, M_data + i, M_stride, this_data, 1);
128  }
129  } else { // The transposed case is slightly more efficient, I guess.
130  for (MatrixIndexT i = 0; i < v_dim; i++) {
131  Real v_i = v.data_[i];
132  if (v_i == 0.0) continue;
133  // Add to *this, the i'th row of the Matrix, times v_i.
134  cblas_Xaxpy(this_dim, v_i * alpha,
135  M_data + (i * M_stride), 1, this_data, 1);
136  }
137  }*/
138 }
139 
140 template<typename Real>
141 void VectorBase<Real>::AddSpVec(const Real alpha,
142  const SpMatrix<Real> &M,
143  const VectorBase<Real> &v,
144  const Real beta) {
145  KALDI_ASSERT(M.NumRows() == v.dim_ && dim_ == v.dim_);
146  KALDI_ASSERT(&v != this);
147  cblas_Xspmv(alpha, M.NumRows(), M.Data(), v.Data(), 1, beta, data_, 1);
148 }
149 
150 
151 template<typename Real>
153  const MatrixTransposeType trans) {
154  KALDI_ASSERT(M.NumRows() == dim_);
155  cblas_Xtpmv(trans,M.Data(),M.NumRows(),data_,1);
156 }
157 
158 template<typename Real>
160  const MatrixTransposeType trans) {
161  KALDI_ASSERT(M.NumRows() == dim_);
162  cblas_Xtpsv(trans, M.Data(), M.NumRows(), data_, 1);
163 }
164 
165 
166 template<typename Real>
167 inline void Vector<Real>::Init(const MatrixIndexT dim) {
168  KALDI_ASSERT(dim >= 0);
169  if (dim == 0) {
170  this->dim_ = 0;
171  this->data_ = NULL;
172  return;
173  }
174  MatrixIndexT size;
175  void *data;
176  void *free_data;
177 
178  size = dim * sizeof(Real);
179 
180  if ((data = KALDI_MEMALIGN(16, size, &free_data)) != NULL) {
181  this->data_ = static_cast<Real*> (data);
182  this->dim_ = dim;
183  } else {
184  throw std::bad_alloc();
185  }
186 }
187 
188 
189 template<typename Real>
190 void Vector<Real>::Resize(const MatrixIndexT dim, MatrixResizeType resize_type) {
191 
192  // the next block uses recursion to handle what we have to do if
193  // resize_type == kCopyData.
194  if (resize_type == kCopyData) {
195  if (this->data_ == NULL || dim == 0) resize_type = kSetZero; // nothing to copy.
196  else if (this->dim_ == dim) { return; } // nothing to do.
197  else {
198  // set tmp to a vector of the desired size.
199  Vector<Real> tmp(dim, kUndefined);
200  if (dim > this->dim_) {
201  memcpy(tmp.data_, this->data_, sizeof(Real)*this->dim_);
202  memset(tmp.data_+this->dim_, 0, sizeof(Real)*(dim-this->dim_));
203  } else {
204  memcpy(tmp.data_, this->data_, sizeof(Real)*dim);
205  }
206  tmp.Swap(this);
207  // and now let tmp go out of scope, deleting what was in *this.
208  return;
209  }
210  }
211  // At this point, resize_type == kSetZero or kUndefined.
212 
213  if (this->data_ != NULL) {
214  if (this->dim_ == dim) {
215  if (resize_type == kSetZero) this->SetZero();
216  return;
217  } else {
218  Destroy();
219  }
220  }
221  Init(dim);
222  if (resize_type == kSetZero) this->SetZero();
223 }
224 
225 
227 template<typename Real>
229  KALDI_ASSERT(Dim() == v.Dim());
230  if (data_ != v.data_) {
231  std::memcpy(this->data_, v.data_, dim_ * sizeof(Real));
232  }
233 }
234 
235 template<typename Real>
236 template<typename OtherReal>
239  this->CopyFromVec(v);
240 }
241 // instantiate the template.
242 template void VectorBase<float>::CopyFromPacked(const PackedMatrix<double> &other);
243 template void VectorBase<float>::CopyFromPacked(const PackedMatrix<float> &other);
244 template void VectorBase<double>::CopyFromPacked(const PackedMatrix<double> &other);
245 template void VectorBase<double>::CopyFromPacked(const PackedMatrix<float> &other);
246 
248 template<typename Real>
249 void VectorBase<Real>::CopyFromPtr(const Real *data, MatrixIndexT sz) {
250  KALDI_ASSERT(dim_ == sz);
251  std::memcpy(this->data_, data, Dim() * sizeof(Real));
252 }
253 
254 template<typename Real>
255 template<typename OtherReal>
257  KALDI_ASSERT(dim_ == other.Dim());
258  Real * __restrict__ ptr = data_;
259  const OtherReal * __restrict__ other_ptr = other.Data();
260  for (MatrixIndexT i = 0; i < dim_; i++)
261  ptr[i] = other_ptr[i];
262 }
263 
264 template void VectorBase<float>::CopyFromVec(const VectorBase<double> &other);
265 template void VectorBase<double>::CopyFromVec(const VectorBase<float> &other);
266 
267 // Remove element from the vector. The vector is not reallocated
268 template<typename Real>
270  KALDI_ASSERT(i < this->dim_ && "Access out of vector");
271  for (MatrixIndexT j = i + 1; j < this->dim_; j++)
272  this->data_[j-1] = this->data_[j];
273  this->dim_--;
274 }
275 
276 
278 template<typename Real>
281  if (this->data_ != NULL)
282  KALDI_MEMALIGN_FREE(this->data_);
283  this->data_ = NULL;
284  this->dim_ = 0;
285 }
286 
287 template<typename Real>
289  std::memset(data_, 0, dim_ * sizeof(Real));
290 }
291 
292 template<typename Real>
293 bool VectorBase<Real>::IsZero(Real cutoff) const {
294  Real abs_max = 0.0;
295  for (MatrixIndexT i = 0; i < Dim(); i++)
296  abs_max = std::max(std::abs(data_[i]), abs_max);
297  return (abs_max <= cutoff);
298 }
299 
300 template<typename Real>
302  kaldi::RandomState rstate;
303  MatrixIndexT last = (Dim() % 2 == 1) ? Dim() - 1 : Dim();
304  for (MatrixIndexT i = 0; i < last; i += 2) {
305  kaldi::RandGauss2(data_ + i, data_ + i + 1, &rstate);
306  }
307  if (Dim() != last) data_[last] = static_cast<Real>(kaldi::RandGauss(&rstate));
308 }
309 
310 template<typename Real>
312  kaldi::RandomState rstate;
313  for (MatrixIndexT i = 0; i < Dim(); i++) {
314  *(data_+i) = RandUniform(&rstate);
315  }
316 }
317 
318 template<typename Real>
320  kaldi::RandomState rstate;
321  Real sum = this->Sum();
322  KALDI_ASSERT(this->Min() >= 0.0 && sum > 0.0);
323  Real r = RandUniform(&rstate) * sum;
324  Real *data = this->data_;
325  MatrixIndexT dim = this->dim_;
326  Real running_sum = 0.0;
327  for (MatrixIndexT i = 0; i < dim; i++) {
328  running_sum += data[i];
329  if (r < running_sum) return i;
330  }
331  return dim_ - 1; // Should only happen if RandUniform()
332  // returns exactly 1, or due to roundoff.
333 }
334 
335 template<typename Real>
336 void VectorBase<Real>::Set(Real f) {
337  // Why not use memset here?
338  // The basic unit of memset is a byte.
339  // If f != 0 and sizeof(Real) > 1, then we cannot use memset.
340  if (f == 0) {
341  this->SetZero(); // calls std::memset
342  } else {
343  for (MatrixIndexT i = 0; i < dim_; i++) { data_[i] = f; }
344  }
345 }
346 
347 template<typename Real>
349  KALDI_ASSERT(dim_ == mat.NumCols() * mat.NumRows());
350 
351  Real *inc_data = data_;
352  const MatrixIndexT cols = mat.NumCols(), rows = mat.NumRows();
353 
354  if (mat.Stride() == mat.NumCols()) {
355  memcpy(inc_data, mat.Data(), cols*rows*sizeof(Real));
356  } else {
357  for (MatrixIndexT i = 0; i < rows; i++) {
358  // copy the data to the propper position
359  memcpy(inc_data, mat.RowData(i), cols * sizeof(Real));
360  // set new copy position
361  inc_data += cols;
362  }
363  }
364 }
365 
366 template<typename Real>
367 template<typename OtherReal>
369  KALDI_ASSERT(dim_ == mat.NumCols() * mat.NumRows());
370  Real *vec_data = data_;
371  const MatrixIndexT cols = mat.NumCols(),
372  rows = mat.NumRows();
373 
374  for (MatrixIndexT i = 0; i < rows; i++) {
375  const OtherReal *mat_row = mat.RowData(i);
376  for (MatrixIndexT j = 0; j < cols; j++) {
377  vec_data[j] = static_cast<Real>(mat_row[j]);
378  }
379  vec_data += cols;
380  }
381 }
382 
383 template
385 template
387 
388 
389 template<typename Real>
391  KALDI_ASSERT(dim_ == mat.NumCols() * mat.NumRows());
392 
393  Real* inc_data = data_;
394  const MatrixIndexT cols = mat.NumCols(), rows = mat.NumRows(), stride = mat.Stride();
395  const Real *mat_inc_data = mat.Data();
396 
397  for (MatrixIndexT i = 0; i < cols; i++) {
398  for (MatrixIndexT j = 0; j < rows; j++) {
399  inc_data[j] = mat_inc_data[j*stride];
400  }
401  mat_inc_data++;
402  inc_data += rows;
403  }
404 }
405 
406 template<typename Real>
408  KALDI_ASSERT(row < mat.NumRows());
409  KALDI_ASSERT(dim_ == mat.NumCols());
410  const Real *mat_row = mat.RowData(row);
411  memcpy(data_, mat_row, sizeof(Real)*dim_);
412 }
413 
414 template<typename Real>
415 template<typename OtherReal>
417  KALDI_ASSERT(row < mat.NumRows());
418  KALDI_ASSERT(dim_ == mat.NumCols());
419  const OtherReal *mat_row = mat.RowData(row);
420  for (MatrixIndexT i = 0; i < dim_; i++)
421  data_[i] = static_cast<Real>(mat_row[i]);
422 }
423 
424 template
426 template
428 
429 template<typename Real>
430 template<typename OtherReal>
432  KALDI_ASSERT(row < sp.NumRows());
433  KALDI_ASSERT(dim_ == sp.NumCols());
434 
435  const OtherReal *sp_data = sp.Data();
436 
437  sp_data += (row*(row+1)) / 2; // takes us to beginning of this row.
438  MatrixIndexT i;
439  for (i = 0; i < row; i++) // copy consecutive elements.
440  data_[i] = static_cast<Real>(*(sp_data++));
441  for(; i < dim_; ++i, sp_data += i)
442  data_[i] = static_cast<Real>(*sp_data);
443 }
444 
445 template
447 template
449 template
450 void VectorBase<float>::CopyRowFromSp(const SpMatrix<float> &mat, MatrixIndexT row);
451 template
452 void VectorBase<double>::CopyRowFromSp(const SpMatrix<double> &mat, MatrixIndexT row);
453 
454 
455 #ifdef HAVE_MKL
456 template<>
457 void VectorBase<float>::Pow(const VectorBase<float> &v, float power) {
458  vsPowx(dim_, data_, power, v.data_);
459 }
460 template<>
461 void VectorBase<double>::Pow(const VectorBase<double> &v, double power) {
462  vdPowx(dim_, data_, power, v.data_);
463 }
464 #else
465 
466 // takes elements to a power. Does not check output.
467 template<typename Real>
468 void VectorBase<Real>::Pow(const VectorBase<Real> &v, Real power) {
469  KALDI_ASSERT(dim_ == v.dim_);
470  for (MatrixIndexT i = 0; i < dim_; i++) {
471  data_[i] = pow(v.data_[i], power);
472  }
473 }
474 #endif
475 
476 // takes absolute value of the elements to a power.
477 // Throws exception if could not (but only for power != 1 and power != 2).
478 template<typename Real>
479 void VectorBase<Real>::ApplyPowAbs(Real power, bool include_sign) {
480  if (power == 1.0)
481  for (MatrixIndexT i = 0; i < dim_; i++)
482  data_[i] = (include_sign && data_[i] < 0 ? -1 : 1) * std::abs(data_[i]);
483  if (power == 2.0) {
484  for (MatrixIndexT i = 0; i < dim_; i++)
485  data_[i] = (include_sign && data_[i] < 0 ? -1 : 1) * data_[i] * data_[i];
486  } else if (power == 0.5) {
487  for (MatrixIndexT i = 0; i < dim_; i++) {
488  data_[i] = (include_sign && data_[i] < 0 ? -1 : 1) * std::sqrt(std::abs(data_[i]));
489  }
490  } else if (power < 0.0) {
491  for (MatrixIndexT i = 0; i < dim_; i++) {
492  data_[i] = (data_[i] == 0.0 ? 0.0 : pow(std::abs(data_[i]), power));
493  data_[i] *= (include_sign && data_[i] < 0 ? -1 : 1);
494  if (data_[i] == HUGE_VAL) { // HUGE_VAL is what errno returns on error.
495  KALDI_ERR << "Could not raise element " << i << "to power "
496  << power << ": returned value = " << data_[i];
497  }
498  }
499  } else {
500  for (MatrixIndexT i = 0; i < dim_; i++) {
501  data_[i] = (include_sign && data_[i] < 0 ? -1 : 1) * pow(std::abs(data_[i]), power);
502  if (data_[i] == HUGE_VAL) { // HUGE_VAL is what errno returns on error.
503  KALDI_ERR << "Could not raise element " << i << "to power "
504  << power << ": returned value = " << data_[i];
505  }
506  }
507  }
508 }
509 
510 // Computes the p-th norm. Throws exception if could not.
511 template<typename Real>
512 Real VectorBase<Real>::Norm(Real p) const {
513  KALDI_ASSERT(p >= 0.0);
514  Real sum = 0.0;
515  if (p == 0.0) {
516  for (MatrixIndexT i = 0; i < dim_; i++)
517  if (data_[i] != 0.0) sum += 1.0;
518  return sum;
519  } else if (p == 1.0) {
520  for (MatrixIndexT i = 0; i < dim_; i++)
521  sum += std::abs(data_[i]);
522  return sum;
523  } else if (p == 2.0) {
524  for (MatrixIndexT i = 0; i < dim_; i++)
525  sum += data_[i] * data_[i];
526  return std::sqrt(sum);
527  } else if (p == std::numeric_limits<Real>::infinity()){
528  for (MatrixIndexT i = 0; i < dim_; i++)
529  sum = std::max(sum, std::abs(data_[i]));
530  return sum;
531  } else {
532  Real tmp;
533  bool ok = true;
534  for (MatrixIndexT i = 0; i < dim_; i++) {
535  tmp = pow(std::abs(data_[i]), p);
536  if (tmp == HUGE_VAL) // HUGE_VAL is what pow returns on error.
537  ok = false;
538  sum += tmp;
539  }
540  tmp = pow(sum, static_cast<Real>(1.0/p));
541  KALDI_ASSERT(tmp != HUGE_VAL); // should not happen here.
542  if (ok) {
543  return tmp;
544  } else {
545  Real maximum = this->Max(), minimum = this->Min(),
546  max_abs = std::max(maximum, -minimum);
547  KALDI_ASSERT(max_abs > 0); // Or should not have reached here.
548  Vector<Real> tmp(*this);
549  tmp.Scale(1.0 / max_abs);
550  return tmp.Norm(p) * max_abs;
551  }
552  }
553 }
554 
555 template<typename Real>
556 bool VectorBase<Real>::ApproxEqual(const VectorBase<Real> &other, float tol) const {
557  if (dim_ != other.dim_) KALDI_ERR << "ApproxEqual: size mismatch "
558  << dim_ << " vs. " << other.dim_;
559  KALDI_ASSERT(tol >= 0.0);
560  if (tol != 0.0) {
561  Vector<Real> tmp(*this);
562  tmp.AddVec(-1.0, other);
563  return (tmp.Norm(2.0) <= static_cast<Real>(tol) * this->Norm(2.0));
564  } else { // Test for exact equality.
565  const Real *data = data_;
566  const Real *other_data = other.data_;
567  for (MatrixIndexT dim = dim_, i = 0; i < dim; i++)
568  if (data[i] != other_data[i]) return false;
569  return true;
570  }
571 }
572 
573 template<typename Real>
574 Real VectorBase<Real>::Max() const {
575  Real ans = - std::numeric_limits<Real>::infinity();
576  const Real *data = data_;
577  MatrixIndexT i, dim = dim_;
578  for (i = 0; i + 4 <= dim; i += 4) {
579  Real a1 = data[i], a2 = data[i+1], a3 = data[i+2], a4 = data[i+3];
580  if (a1 > ans || a2 > ans || a3 > ans || a4 > ans) {
581  Real b1 = (a1 > a2 ? a1 : a2), b2 = (a3 > a4 ? a3 : a4);
582  if (b1 > ans) ans = b1;
583  if (b2 > ans) ans = b2;
584  }
585  }
586  for (; i < dim; i++)
587  if (data[i] > ans) ans = data[i];
588  return ans;
589 }
590 
591 template<typename Real>
592 Real VectorBase<Real>::Max(MatrixIndexT *index_out) const {
593  if (dim_ == 0) KALDI_ERR << "Empty vector";
594  Real ans = - std::numeric_limits<Real>::infinity();
595  MatrixIndexT index = 0;
596  const Real *data = data_;
597  MatrixIndexT i, dim = dim_;
598  for (i = 0; i + 4 <= dim; i += 4) {
599  Real a1 = data[i], a2 = data[i+1], a3 = data[i+2], a4 = data[i+3];
600  if (a1 > ans || a2 > ans || a3 > ans || a4 > ans) {
601  if (a1 > ans) { ans = a1; index = i; }
602  if (a2 > ans) { ans = a2; index = i + 1; }
603  if (a3 > ans) { ans = a3; index = i + 2; }
604  if (a4 > ans) { ans = a4; index = i + 3; }
605  }
606  }
607  for (; i < dim; i++)
608  if (data[i] > ans) { ans = data[i]; index = i; }
609  *index_out = index;
610  return ans;
611 }
612 
613 template<typename Real>
614 Real VectorBase<Real>::Min() const {
615  Real ans = std::numeric_limits<Real>::infinity();
616  const Real *data = data_;
617  MatrixIndexT i, dim = dim_;
618  for (i = 0; i + 4 <= dim; i += 4) {
619  Real a1 = data[i], a2 = data[i+1], a3 = data[i+2], a4 = data[i+3];
620  if (a1 < ans || a2 < ans || a3 < ans || a4 < ans) {
621  Real b1 = (a1 < a2 ? a1 : a2), b2 = (a3 < a4 ? a3 : a4);
622  if (b1 < ans) ans = b1;
623  if (b2 < ans) ans = b2;
624  }
625  }
626  for (; i < dim; i++)
627  if (data[i] < ans) ans = data[i];
628  return ans;
629 }
630 
631 template<typename Real>
632 Real VectorBase<Real>::Min(MatrixIndexT *index_out) const {
633  if (dim_ == 0) KALDI_ERR << "Empty vector";
634  Real ans = std::numeric_limits<Real>::infinity();
635  MatrixIndexT index = 0;
636  const Real *data = data_;
637  MatrixIndexT i, dim = dim_;
638  for (i = 0; i + 4 <= dim; i += 4) {
639  Real a1 = data[i], a2 = data[i+1], a3 = data[i+2], a4 = data[i+3];
640  if (a1 < ans || a2 < ans || a3 < ans || a4 < ans) {
641  if (a1 < ans) { ans = a1; index = i; }
642  if (a2 < ans) { ans = a2; index = i + 1; }
643  if (a3 < ans) { ans = a3; index = i + 2; }
644  if (a4 < ans) { ans = a4; index = i + 3; }
645  }
646  }
647  for (; i < dim; i++)
648  if (data[i] < ans) { ans = data[i]; index = i; }
649  *index_out = index;
650  return ans;
651 }
652 
653 
654 template<typename Real>
655 template<typename OtherReal>
657  KALDI_ASSERT(col < mat.NumCols());
658  KALDI_ASSERT(dim_ == mat.NumRows());
659  for (MatrixIndexT i = 0; i < dim_; i++)
660  data_[i] = mat(i, col);
661  // can't do this very efficiently so don't really bother. could improve this though.
662 }
663 // instantiate the template above.
664 template
666 template
667 void VectorBase<float>::CopyColFromMat(const MatrixBase<double> &mat, MatrixIndexT col);
668 template
670 template
671 void VectorBase<double>::CopyColFromMat(const MatrixBase<double> &mat, MatrixIndexT col);
672 
673 template<typename Real>
675  KALDI_ASSERT(dim_ == std::min(M.NumRows(), M.NumCols()));
676  cblas_Xcopy(dim_, M.Data(), M.Stride() + 1, data_, 1);
677 }
678 
679 template<typename Real>
681  KALDI_ASSERT(dim_ == M.NumCols());
682  for (MatrixIndexT i = 0; i < dim_; i++)
683  data_[i] = M(i, i);
684  // could make this more efficient.
685 }
686 
687 template<typename Real>
688 Real VectorBase<Real>::Sum() const {
689  // Do a dot-product with a size-1 array with a stride of 0 to
690  // implement sum. This allows us to access SIMD operations in a
691  // cross-platform way via your BLAS library.
692  Real one(1);
693  return cblas_Xdot(dim_, data_, 1, &one, 0);
694 }
695 
696 template<typename Real>
698  double sum_log = 0.0;
699  double prod = 1.0;
700  for (MatrixIndexT i = 0; i < dim_; i++) {
701  prod *= data_[i];
702  // Possible future work (arnab): change these magic values to pre-defined
703  // constants
704  if (prod < 1.0e-10 || prod > 1.0e+10) {
705  sum_log += Log(prod);
706  prod = 1.0;
707  }
708  }
709  if (prod != 1.0) sum_log += Log(prod);
710  return sum_log;
711 }
712 
713 template<typename Real>
714 void VectorBase<Real>::AddRowSumMat(Real alpha, const MatrixBase<Real> &M, Real beta) {
715  KALDI_ASSERT(dim_ == M.NumCols());
716  MatrixIndexT num_rows = M.NumRows(), stride = M.Stride(), dim = dim_;
717  Real *data = data_;
718 
719  // implement the function according to a dimension cutoff for computation efficiency
720  if (num_rows <= 64) {
721  cblas_Xscal(dim, beta, data, 1);
722  const Real *m_data = M.Data();
723  for (MatrixIndexT i = 0; i < num_rows; i++, m_data += stride)
724  cblas_Xaxpy(dim, alpha, m_data, 1, data, 1);
725 
726  } else {
727  Vector<Real> ones(M.NumRows());
728  ones.Set(1.0);
729  this->AddMatVec(alpha, M, kTrans, ones, beta);
730  }
731 }
732 
733 template<typename Real>
734 void VectorBase<Real>::AddColSumMat(Real alpha, const MatrixBase<Real> &M, Real beta) {
735  KALDI_ASSERT(dim_ == M.NumRows());
736  MatrixIndexT num_cols = M.NumCols();
737 
738  // implement the function according to a dimension cutoff for computation efficiency
739  if (num_cols <= 64) {
740  for (MatrixIndexT i = 0; i < dim_; i++) {
741  double sum = 0.0;
742  const Real *src = M.RowData(i);
743  for (MatrixIndexT j = 0; j < num_cols; j++)
744  sum += src[j];
745  data_[i] = alpha * sum + beta * data_[i];
746  }
747  } else {
748  Vector<Real> ones(M.NumCols());
749  ones.Set(1.0);
750  this->AddMatVec(alpha, M, kNoTrans, ones, beta);
751  }
752 }
753 
754 template<typename Real>
755 Real VectorBase<Real>::LogSumExp(Real prune) const {
756  Real sum;
757  if (sizeof(sum) == 8) sum = kLogZeroDouble;
758  else sum = kLogZeroFloat;
759  Real max_elem = Max(), cutoff;
760  if (sizeof(Real) == 4) cutoff = max_elem + kMinLogDiffFloat;
761  else cutoff = max_elem + kMinLogDiffDouble;
762  if (prune > 0.0 && max_elem - prune > cutoff) // explicit pruning...
763  cutoff = max_elem - prune;
764 
765  double sum_relto_max_elem = 0.0;
766 
767  for (MatrixIndexT i = 0; i < dim_; i++) {
768  BaseFloat f = data_[i];
769  if (f >= cutoff)
770  sum_relto_max_elem += Exp(f - max_elem);
771  }
772  return max_elem + Log(sum_relto_max_elem);
773 }
774 
775 template<typename Real>
777  for (MatrixIndexT i = 0; i < dim_; i++) {
778  data_[i] = static_cast<Real>(1 / data_[i]);
779  }
780 }
781 
782 template<typename Real>
784  for (MatrixIndexT i = 0; i < dim_; i++) {
785  if (data_[i] < 0.0)
786  KALDI_ERR << "Trying to take log of a negative number.";
787  data_[i] = Log(data_[i]);
788  }
789 }
790 
791 template<typename Real>
793  KALDI_ASSERT(dim_ == v.Dim());
794  for (MatrixIndexT i = 0; i < dim_; i++) {
795  data_[i] = Log(v(i));
796  }
797 }
798 
799 template<typename Real>
801  for (MatrixIndexT i = 0; i < dim_; i++) {
802  data_[i] = Exp(data_[i]);
803  }
804 }
805 
806 template<typename Real>
808  for (MatrixIndexT i = 0; i < dim_; i++) { data_[i] = std::abs(data_[i]); }
809 }
810 
811 template<typename Real>
812 void VectorBase<Real>::Floor(const VectorBase<Real> &v, Real floor_val, MatrixIndexT *floored_count) {
813  KALDI_ASSERT(dim_ == v.dim_);
814  if (floored_count == nullptr) {
815  for (MatrixIndexT i = 0; i < dim_; i++) {
816  data_[i] = std::max(v.data_[i], floor_val);
817  }
818  } else {
819  MatrixIndexT num_floored = 0;
820  for (MatrixIndexT i = 0; i < dim_; i++) {
821  if (v.data_[i] < floor_val) {
822  data_[i] = floor_val;
823  num_floored++;
824  } else {
825  data_[i] = v.data_[i];
826  }
827  }
828  *floored_count = num_floored;
829  }
830 }
831 
832 template<typename Real>
833 void VectorBase<Real>::Ceiling(const VectorBase<Real> &v, Real ceil_val, MatrixIndexT *ceiled_count) {
834  KALDI_ASSERT(dim_ == v.dim_);
835  if (ceiled_count == nullptr) {
836  for (MatrixIndexT i = 0; i < dim_; i++) {
837  data_[i] = std::min(v.data_[i], ceil_val);
838  }
839  } else {
840  MatrixIndexT num_changed = 0;
841  for (MatrixIndexT i = 0; i < dim_; i++) {
842  if (v.data_[i] > ceil_val) {
843  data_[i] = ceil_val;
844  num_changed++;
845  } else {
846  data_[i] = v.data_[i];
847  }
848  }
849  *ceiled_count = num_changed;
850  }
851 }
852 
853 template<typename Real>
855  KALDI_ASSERT(floor_vec.Dim() == dim_);
856  MatrixIndexT num_floored = 0;
857  for (MatrixIndexT i = 0; i < dim_; i++) {
858  if (data_[i] < floor_vec(i)) {
859  data_[i] = floor_vec(i);
860  num_floored++;
861  }
862  }
863  return num_floored;
864 }
865 
866 template<typename Real>
868  Real max = this->Max(), sum = 0.0;
869  for (MatrixIndexT i = 0; i < dim_; i++) {
870  sum += (data_[i] = Exp(data_[i] - max));
871  }
872  this->Scale(1.0 / sum);
873  return max + Log(sum);
874 }
875 
876 template<typename Real>
878  Real max = this->Max(), sum = 0.0;
879  for (MatrixIndexT i = 0; i < dim_; i++) {
880  sum += Exp((data_[i] -= max));
881  }
882  sum = Log(sum);
883  this->Add(-1.0 * sum);
884  return max + sum;
885 }
886 
887 #ifdef HAVE_MKL
888 template<>
890  KALDI_ASSERT(dim_ == src.dim_);
891  vsTanh(dim_, src.data_, data_);
892 }
893 template<>
895  KALDI_ASSERT(dim_ == src.dim_);
896  vdTanh(dim_, src.data_, data_);
897 }
898 #else
899 template<typename Real>
901  KALDI_ASSERT(dim_ == src.dim_);
902  for (MatrixIndexT i = 0; i < dim_; i++) {
903  Real x = src.data_[i];
904  if (x > 0.0) {
905  Real inv_expx = Exp(-x);
906  x = -1.0 + 2.0 / (1.0 + inv_expx * inv_expx);
907  } else {
908  Real expx = Exp(x);
909  x = 1.0 - 2.0 / (1.0 + expx * expx);
910  }
911  data_[i] = x;
912  }
913 }
914 #endif
915 
916 #ifdef HAVE_MKL
917 // Implementing sigmoid based on tanh.
918 template<>
920  KALDI_ASSERT(dim_ == src.dim_);
921  this->CopyFromVec(src);
922  this->Scale(0.5);
923  vsTanh(dim_, data_, data_);
924  this->Add(1.0);
925  this->Scale(0.5);
926 }
927 template<>
929  KALDI_ASSERT(dim_ == src.dim_);
930  this->CopyFromVec(src);
931  this->Scale(0.5);
932  vdTanh(dim_, data_, data_);
933  this->Add(1.0);
934  this->Scale(0.5);
935 }
936 #else
937 template<typename Real>
939  KALDI_ASSERT(dim_ == src.dim_);
940  for (MatrixIndexT i = 0; i < dim_; i++) {
941  Real x = src.data_[i];
942  // We aim to avoid floating-point overflow here.
943  if (x > 0.0) {
944  x = 1.0 / (1.0 + Exp(-x));
945  } else {
946  Real ex = Exp(x);
947  x = ex / (ex + 1.0);
948  }
949  data_[i] = x;
950  }
951 }
952 #endif
953 
954 
955 template<typename Real>
956 void VectorBase<Real>::Add(Real c) {
957  for (MatrixIndexT i = 0; i < dim_; i++) {
958  data_[i] += c;
959  }
960 }
961 
962 template<typename Real>
963 void VectorBase<Real>::Scale(Real alpha) {
964  cblas_Xscal(dim_, alpha, data_, 1);
965 }
966 
967 template<typename Real>
969  KALDI_ASSERT(dim_ == v.dim_);
970  for (MatrixIndexT i = 0; i < dim_; i++) {
971  data_[i] *= v.data_[i];
972  }
973 }
974 
975 template<typename Real> // Set each element to y = (x == orig ? changed : x).
976 void VectorBase<Real>::ReplaceValue(Real orig, Real changed) {
977  Real *data = data_;
978  for (MatrixIndexT i = 0; i < dim_; i++)
979  if (data[i] == orig) data[i] = changed;
980 }
981 
982 
983 template<typename Real>
984 template<typename OtherReal>
986  KALDI_ASSERT(dim_ == v.Dim());
987  const OtherReal *other_ptr = v.Data();
988  for (MatrixIndexT i = 0; i < dim_; i++) {
989  data_[i] *= other_ptr[i];
990  }
991 }
992 // instantiate template.
993 template
995 template
997 
998 
999 template<typename Real>
1001  const VectorBase<Real> &r, Real beta) {
1002  KALDI_ASSERT(v.data_ != this->data_ && r.data_ != this->data_);
1003  // We pretend that v is a band-diagonal matrix.
1004  KALDI_ASSERT(dim_ == v.dim_ && dim_ == r.dim_);
1005  cblas_Xgbmv(kNoTrans, dim_, dim_, 0, 0, alpha, v.data_, 1,
1006  r.data_, 1, beta, this->data_, 1);
1007 }
1008 
1009 
1010 template<typename Real>
1012  KALDI_ASSERT(dim_ == v.dim_);
1013  for (MatrixIndexT i = 0; i < dim_; i++) {
1014  data_[i] /= v.data_[i];
1015  }
1016 }
1017 
1018 template<typename Real>
1019 template<typename OtherReal>
1021  KALDI_ASSERT(dim_ == v.Dim());
1022  const OtherReal *other_ptr = v.Data();
1023  for (MatrixIndexT i = 0; i < dim_; i++) {
1024  data_[i] /= other_ptr[i];
1025  }
1026 }
1027 // instantiate template.
1028 template
1030 template
1032 
1033 template<typename Real>
1035  const VectorBase<Real> &rr, Real beta) {
1036  KALDI_ASSERT((dim_ == v.dim_ && dim_ == rr.dim_));
1037  for (MatrixIndexT i = 0; i < dim_; i++) {
1038  data_[i] = alpha * v.data_[i]/rr.data_[i] + beta * data_[i] ;
1039  }
1040 }
1041 
1042 template<typename Real>
1043 template<typename OtherReal>
1044 void VectorBase<Real>::AddVec(const Real alpha, const VectorBase<OtherReal> &v) {
1045  KALDI_ASSERT(dim_ == v.dim_);
1046  // remove __restrict__ if it causes compilation problems.
1047  Real *__restrict__ data = data_;
1048  OtherReal *__restrict__ other_data = v.data_;
1049  MatrixIndexT dim = dim_;
1050  if (alpha != 1.0)
1051  for (MatrixIndexT i = 0; i < dim; i++)
1052  data[i] += alpha * other_data[i];
1053  else
1054  for (MatrixIndexT i = 0; i < dim; i++)
1055  data[i] += other_data[i];
1056 }
1057 
1058 template
1059 void VectorBase<float>::AddVec(const float alpha, const VectorBase<double> &v);
1060 template
1061 void VectorBase<double>::AddVec(const double alpha, const VectorBase<float> &v);
1062 
1063 template<typename Real>
1064 template<typename OtherReal>
1065 void VectorBase<Real>::AddVec2(const Real alpha, const VectorBase<OtherReal> &v) {
1066  KALDI_ASSERT(dim_ == v.dim_);
1067  // remove __restrict__ if it causes compilation problems.
1068  Real *__restrict__ data = data_;
1069  OtherReal *__restrict__ other_data = v.data_;
1070  MatrixIndexT dim = dim_;
1071  if (alpha != 1.0)
1072  for (MatrixIndexT i = 0; i < dim; i++)
1073  data[i] += alpha * other_data[i] * other_data[i];
1074  else
1075  for (MatrixIndexT i = 0; i < dim; i++)
1076  data[i] += other_data[i] * other_data[i];
1077 }
1078 
1079 template
1080 void VectorBase<float>::AddVec2(const float alpha, const VectorBase<double> &v);
1081 template
1082 void VectorBase<double>::AddVec2(const double alpha, const VectorBase<float> &v);
1083 
1084 
1085 template<typename Real>
1086 void VectorBase<Real>::Read(std::istream &is, bool binary, bool add) {
1087  if (add) {
1088  Vector<Real> tmp(Dim());
1089  tmp.Read(is, binary, false); // read without adding.
1090  if (this->Dim() != tmp.Dim()) {
1091  KALDI_ERR << "VectorBase::Read, size mismatch " << this->Dim()<<" vs. "<<tmp.Dim();
1092  }
1093  this->AddVec(1.0, tmp);
1094  return;
1095  } // now assume add == false.
1096 
1097  // In order to avoid rewriting this, we just declare a Vector and
1098  // use it to read the data, then copy.
1099  Vector<Real> tmp;
1100  tmp.Read(is, binary, false);
1101  if (tmp.Dim() != Dim())
1102  KALDI_ERR << "VectorBase<Real>::Read, size mismatch "
1103  << Dim() << " vs. " << tmp.Dim();
1104  CopyFromVec(tmp);
1105 }
1106 
1107 
1108 template<typename Real>
1109 void Vector<Real>::Read(std::istream &is, bool binary, bool add) {
1110  if (add) {
1111  Vector<Real> tmp(this->Dim());
1112  tmp.Read(is, binary, false); // read without adding.
1113  if (this->Dim() == 0) this->Resize(tmp.Dim());
1114  if (this->Dim() != tmp.Dim()) {
1115  KALDI_ERR << "Vector<Real>::Read, adding but dimensions mismatch "
1116  << this->Dim() << " vs. " << tmp.Dim();
1117  }
1118  this->AddVec(1.0, tmp);
1119  return;
1120  } // now assume add == false.
1121 
1122  std::ostringstream specific_error;
1123  MatrixIndexT pos_at_start = is.tellg();
1124 
1125  if (binary) {
1126  int peekval = Peek(is, binary);
1127  const char *my_token = (sizeof(Real) == 4 ? "FV" : "DV");
1128  char other_token_start = (sizeof(Real) == 4 ? 'D' : 'F');
1129  if (peekval == other_token_start) { // need to instantiate the other type to read it.
1130  typedef typename OtherReal<Real>::Real OtherType; // if Real == float, OtherType == double, and vice versa.
1131  Vector<OtherType> other(this->Dim());
1132  other.Read(is, binary, false); // add is false at this point.
1133  if (this->Dim() != other.Dim()) this->Resize(other.Dim());
1134  this->CopyFromVec(other);
1135  return;
1136  }
1137  std::string token;
1138  ReadToken(is, binary, &token);
1139  if (token != my_token) {
1140  if (token.length() > 20) token = token.substr(0, 17) + "...";
1141  specific_error << ": Expected token " << my_token << ", got " << token;
1142  goto bad;
1143  }
1144  int32 size;
1145  ReadBasicType(is, binary, &size); // throws on error.
1146  if ((MatrixIndexT)size != this->Dim()) this->Resize(size);
1147  if (size > 0)
1148  is.read(reinterpret_cast<char*>(this->data_), sizeof(Real)*size);
1149  if (is.fail()) {
1150  specific_error << "Error reading vector data (binary mode); truncated "
1151  "stream? (size = " << size << ")";
1152  goto bad;
1153  }
1154  return;
1155  } else { // Text mode reading; format is " [ 1.1 2.0 3.4 ]\n"
1156  std::string s;
1157  is >> s;
1158  // if ((s.compare("DV") == 0) || (s.compare("FV") == 0)) { // Back compatibility.
1159  // is >> s; // get dimension
1160  // is >> s; // get "["
1161  // }
1162  if (is.fail()) { specific_error << "EOF while trying to read vector."; goto bad; }
1163  if (s.compare("[]") == 0) { Resize(0); return; } // tolerate this variant.
1164  if (s.compare("[")) {
1165  if (s.length() > 20) s = s.substr(0, 17) + "...";
1166  specific_error << "Expected \"[\" but got " << s;
1167  goto bad;
1168  }
1169  std::vector<Real> data;
1170  while (1) {
1171  int i = is.peek();
1172  if (i == '-' || (i >= '0' && i <= '9')) { // common cases first.
1173  Real r;
1174  is >> r;
1175  if (is.fail()) { specific_error << "Failed to read number."; goto bad; }
1176  if (! std::isspace(is.peek()) && is.peek() != ']') {
1177  specific_error << "Expected whitespace after number."; goto bad;
1178  }
1179  data.push_back(r);
1180  // But don't eat whitespace... we want to check that it's not newlines
1181  // which would be valid only for a matrix.
1182  } else if (i == ' ' || i == '\t') {
1183  is.get();
1184  } else if (i == ']') {
1185  is.get(); // eat the ']'
1186  this->Resize(data.size());
1187  for (size_t j = 0; j < data.size(); j++)
1188  this->data_[j] = data[j];
1189  i = is.peek();
1190  if (static_cast<char>(i) == '\r') {
1191  is.get();
1192  is.get(); // get \r\n (must eat what we wrote)
1193  } else if (static_cast<char>(i) == '\n') { is.get(); } // get \n (must eat what we wrote)
1194  if (is.fail()) {
1195  KALDI_WARN << "After end of vector data, read error.";
1196  // we got the data we needed, so just warn for this error.
1197  }
1198  return; // success.
1199  } else if (i == -1) {
1200  specific_error << "EOF while reading vector data.";
1201  goto bad;
1202  } else if (i == '\n' || i == '\r') {
1203  specific_error << "Newline found while reading vector (maybe it's a matrix?)";
1204  goto bad;
1205  } else {
1206  is >> s; // read string.
1207  if (!KALDI_STRCASECMP(s.c_str(), "inf") ||
1208  !KALDI_STRCASECMP(s.c_str(), "infinity")) {
1209  data.push_back(std::numeric_limits<Real>::infinity());
1210  KALDI_WARN << "Reading infinite value into vector.";
1211  } else if (!KALDI_STRCASECMP(s.c_str(), "nan")) {
1212  data.push_back(std::numeric_limits<Real>::quiet_NaN());
1213  KALDI_WARN << "Reading NaN value into vector.";
1214  } else {
1215  if (s.length() > 20) s = s.substr(0, 17) + "...";
1216  specific_error << "Expecting numeric vector data, got " << s;
1217  goto bad;
1218  }
1219  }
1220  }
1221  }
1222  // we never reach this line (the while loop returns directly).
1223 bad:
1224  KALDI_ERR << "Failed to read vector from stream. " << specific_error.str()
1225  << " File position at start is "
1226  << pos_at_start<<", currently "<<is.tellg();
1227 }
1228 
1229 
1230 template<typename Real>
1231 void VectorBase<Real>::Write(std::ostream & os, bool binary) const {
1232  if (!os.good()) {
1233  KALDI_ERR << "Failed to write vector to stream: stream not good";
1234  }
1235  if (binary) {
1236  std::string my_token = (sizeof(Real) == 4 ? "FV" : "DV");
1237  WriteToken(os, binary, my_token);
1238 
1239  int32 size = Dim(); // make the size 32-bit on disk.
1240  KALDI_ASSERT(Dim() == (MatrixIndexT) size);
1241  WriteBasicType(os, binary, size);
1242  os.write(reinterpret_cast<const char*>(Data()), sizeof(Real) * size);
1243  } else {
1244  os << " [ ";
1245  for (MatrixIndexT i = 0; i < Dim(); i++)
1246  os << (*this)(i) << " ";
1247  os << "]\n";
1248  }
1249  if (!os.good())
1250  KALDI_ERR << "Failed to write vector to stream";
1251 }
1252 
1253 
1254 template<typename Real>
1255 void VectorBase<Real>::AddVec2(const Real alpha, const VectorBase<Real> &v) {
1256  KALDI_ASSERT(dim_ == v.dim_);
1257  for (MatrixIndexT i = 0; i < dim_; i++)
1258  data_[i] += alpha * v.data_[i] * v.data_[i];
1259 }
1260 
1261 // this <-- beta*this + alpha*M*v.
1262 template<typename Real>
1263 void VectorBase<Real>::AddTpVec(const Real alpha, const TpMatrix<Real> &M,
1264  const MatrixTransposeType trans,
1265  const VectorBase<Real> &v,
1266  const Real beta) {
1267  KALDI_ASSERT(dim_ == v.dim_ && dim_ == M.NumRows());
1268  if (beta == 0.0) {
1269  if (&v != this) CopyFromVec(v);
1270  MulTp(M, trans);
1271  if (alpha != 1.0) Scale(alpha);
1272  } else {
1273  Vector<Real> tmp(v);
1274  tmp.MulTp(M, trans);
1275  if (beta != 1.0) Scale(beta); // *this <-- beta * *this
1276  AddVec(alpha, tmp); // *this += alpha * M * v
1277  }
1278 }
1279 
1280 template<typename Real>
1282  const VectorBase<Real> &v2) {
1283  KALDI_ASSERT(v1.Dim() == M.NumRows() && v2.Dim() == M.NumCols());
1284  Vector<Real> vtmp(M.NumRows());
1285  vtmp.AddMatVec(1.0, M, kNoTrans, v2, 0.0);
1286  return VecVec(v1, vtmp);
1287 }
1288 
1289 template
1290 float VecMatVec(const VectorBase<float> &v1, const MatrixBase<float> &M,
1291  const VectorBase<float> &v2);
1292 template
1293 double VecMatVec(const VectorBase<double> &v1, const MatrixBase<double> &M,
1294  const VectorBase<double> &v2);
1295 
1296 template<typename Real>
1298  std::swap(this->data_, other->data_);
1299  std::swap(this->dim_, other->dim_);
1300 }
1301 
1302 
1303 template<typename Real>
1305  Real alpha, const MatrixBase<Real> &M,
1306  MatrixTransposeType trans, Real beta) {
1307  if (trans == kNoTrans) {
1308  KALDI_ASSERT(this->dim_ == M.NumRows());
1309  MatrixIndexT rows = this->dim_, cols = M.NumCols(),
1310  mat_stride = M.Stride();
1311  Real *data = this->data_;
1312  const Real *mat_data = M.Data();
1313  for (MatrixIndexT i = 0; i < rows; i++, mat_data += mat_stride, data++)
1314  *data = beta * *data + alpha * cblas_Xdot(cols,mat_data,1,mat_data,1);
1315  } else {
1316  KALDI_ASSERT(this->dim_ == M.NumCols());
1317  MatrixIndexT rows = M.NumRows(), cols = this->dim_,
1318  mat_stride = M.Stride();
1319  Real *data = this->data_;
1320  const Real *mat_data = M.Data();
1321  for (MatrixIndexT i = 0; i < cols; i++, mat_data++, data++)
1322  *data = beta * *data + alpha * cblas_Xdot(rows, mat_data, mat_stride,
1323  mat_data, mat_stride);
1324  }
1325 }
1326 
1327 template<typename Real>
1329  Real alpha,
1330  const MatrixBase<Real> &M, MatrixTransposeType transM,
1331  const MatrixBase<Real> &N, MatrixTransposeType transN,
1332  Real beta) {
1333  MatrixIndexT dim = this->dim_,
1334  M_col_dim = (transM == kTrans ? M.NumRows() : M.NumCols()),
1335  N_row_dim = (transN == kTrans ? N.NumCols() : N.NumRows());
1336  KALDI_ASSERT(M_col_dim == N_row_dim); // this is the dimension we sum over
1337  MatrixIndexT M_row_stride = M.Stride(), M_col_stride = 1;
1338  if (transM == kTrans) std::swap(M_row_stride, M_col_stride);
1339  MatrixIndexT N_row_stride = N.Stride(), N_col_stride = 1;
1340  if (transN == kTrans) std::swap(N_row_stride, N_col_stride);
1341 
1342  Real *data = this->data_;
1343  const Real *Mdata = M.Data(), *Ndata = N.Data();
1344  for (MatrixIndexT i = 0; i < dim; i++, Mdata += M_row_stride, Ndata += N_col_stride, data++) {
1345  *data = beta * *data + alpha * cblas_Xdot(M_col_dim, Mdata, M_col_stride, Ndata, N_row_stride);
1346  }
1347 }
1348 
1349 
1350 template class Vector<float>;
1351 template class Vector<double>;
1352 template class VectorBase<float>;
1353 template class VectorBase<double>;
1354 
1355 } // namespace kaldi
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
This class provides a way for switching between double and float types.
Definition: matrix-common.h:84
Packed symetric matrix class.
Definition: matrix-common.h:62
template double VecMatVec(const VectorBase< double > &v1, const MatrixBase< double > &M, const VectorBase< double > &v2)
MatrixResizeType
Definition: matrix-common.h:37
void cblas_Xtpmv(MatrixTransposeType trans, const float *Mdata, const int num_rows, float *y, const int y_inc)
float RandUniform(struct RandomState *state=NULL)
Returns a random number strictly between 0 and 1.
Definition: kaldi-math.h:151
void cblas_Xtpsv(MatrixTransposeType trans, const float *Mdata, const int num_rows, float *y, const int y_inc)
MatrixIndexT NumCols() const
Returns number of columns (or zero for empty matrix).
Definition: kaldi-matrix.h:67
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
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)
Real * RowData(MatrixIndexT i)
Returns pointer to data for one row (non-const)
Definition: kaldi-matrix.h:87
void swap(basic_filebuf< CharT, Traits > &x, basic_filebuf< CharT, Traits > &y)
float RandGauss(struct RandomState *state=NULL)
Definition: kaldi-math.h:155
kaldi::int32 int32
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
void MulTp(const TpMatrix< Real > &M, const MatrixTransposeType trans)
Multiplies this vector by lower-triangular matrix: *this <– *this *M.
void cblas_Xcopy(const int N, const float *X, const int incX, float *Y, const int incY)
MatrixIndexT NumRows() const
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_
Real Norm(Real p) const
Compute the p-th norm of the vector.
float cblas_Xdot(const int N, const float *const X, const int incX, const float *const Y, const int incY)
MatrixIndexT NumCols() const
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
MatrixIndexT dim_
dimension of vector
Definition: kaldi-vector.h:397
Packed matrix: base class for triangular and symmetric matrices.
Definition: matrix-common.h:64
void cblas_Xscal(const int N, const float alpha, float *data, const int inc)
double Log(double x)
Definition: kaldi-math.h:100
void Swap(Vector< Real > *other)
Swaps the contents of *this and *other. Shallow swap.
#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
#define KALDI_ERR
Definition: kaldi-error.h:147
#define KALDI_MEMALIGN_FREE(x)
Definition: kaldi-utils.h:60
Packed symetric matrix class.
Definition: matrix-common.h:63
#define KALDI_WARN
Definition: kaldi-error.h:150
template double VecVec(const VectorBase< double > &ra, const VectorBase< float > &rb)
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
MatrixIndexT Dim() const
Returns the dimension of the vector.
Definition: kaldi-vector.h:64
void Scale(Real alpha)
Multiplies all elements by this constant.
void AddMatVec(const Real alpha, const MatrixBase< Real > &M, const MatrixTransposeType trans, const VectorBase< Real > &v, const Real beta)
Add matrix times vector : this <– beta*this + alpha*M*v.
Definition: kaldi-vector.cc:92
static const float kMinLogDiffFloat
Definition: kaldi-math.h:125
Real * data_
data memory area
Definition: kaldi-vector.h:395
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)
#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.
#define KALDI_STRCASECMP
Definition: kaldi-utils.h:147
void cblas_Xgemv(MatrixTransposeType trans, MatrixIndexT num_rows, MatrixIndexT num_cols, float alpha, const float *Mdata, MatrixIndexT stride, const float *xdata, MatrixIndexT incX, float beta, float *ydata, MatrixIndexT incY)
MatrixTransposeType
Definition: matrix-common.h:32
void cblas_Xspmv(const float alpha, const int num_rows, const float *Mdata, const float *v, const int v_inc, const float beta, float *y, const int y_inc)
void cblas_Xgbmv(MatrixTransposeType trans, MatrixIndexT num_rows, MatrixIndexT num_cols, MatrixIndexT num_below, MatrixIndexT num_above, float alpha, const float *Mdata, MatrixIndexT stride, const float *xdata, MatrixIndexT incX, float beta, float *ydata, MatrixIndexT incY)
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
const float kLogZeroFloat
Definition: kaldi-math.h:128
Provides a vector abstraction class.
Definition: kaldi-vector.h:41
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 AddVec(const Real alpha, const VectorBase< OtherReal > &v)
Add vector : *this = *this + alpha * rv (with casting between floats and doubles) ...
void Read(std::istream &in, bool binary, bool add=false)
Read function using C++ streams.
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