tp-matrix.cc
4.66 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
// matrix/tp-matrix.cc
// Copyright 2009-2011 Ondrej Glembek; Lukas Burget; Microsoft Corporation
// Saarland University; Yanmin Qian; Haihua Xu
// See ../../COPYING for clarification regarding multiple authors
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
// http://www.apache.org/licenses/LICENSE-2.0
// THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
// KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
// WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
// MERCHANTABLITY OR NON-INFRINGEMENT.
// See the Apache 2 License for the specific language governing permissions and
// limitations under the License.
#include "matrix/tp-matrix.h"
#include "matrix/sp-matrix.h"
#include "matrix/kaldi-matrix.h"
#include "matrix/cblas-wrappers.h"
namespace kaldi {
#ifndef HAVE_ATLAS
template<typename Real>
void TpMatrix<Real>::Invert() {
// these are CLAPACK types
KaldiBlasInt result;
KaldiBlasInt rows = static_cast<int>(this->num_rows_);
// clapack call
// NOTE: Even though "U" is for upper, lapack assumes column-wise storage
// of the data. We have a row-wise storage, therefore, we need to "invert"
clapack_Xtptri(&rows, this->data_, &result);
if (result < 0) {
KALDI_ERR << "Call to CLAPACK stptri_ function failed";
} else if (result > 0) {
KALDI_ERR << "Matrix is singular";
}
}
#else
template<typename Real>
void TpMatrix<Real>::Invert() {
// ATLAS doesn't implement triangular matrix inversion in packed
// format, so we temporarily put in non-packed format.
Matrix<Real> tmp(*this);
int rows = static_cast<int>(this->num_rows_);
// ATLAS call. It's really row-major ordering and a lower triangular matrix,
// but there is some weirdness with Fortran-style indexing that we need to
// take account of, so everything gets swapped.
int result = clapack_Xtrtri( rows, tmp.Data(), tmp.Stride());
// Let's hope ATLAS has the same return value conventions as clapack.
// I couldn't find any documentation online.
if (result < 0) {
KALDI_ERR << "Call to ATLAS strtri function failed";
} else if (result > 0) {
KALDI_ERR << "Matrix is singular";
}
(*this).CopyFromMat(tmp);
}
#endif
template<typename Real>
Real TpMatrix<Real>::Determinant() {
double det = 1.0;
for (MatrixIndexT i = 0; i<this->NumRows(); i++) {
det *= (*this)(i, i);
}
return static_cast<Real>(det);
}
template<typename Real>
void TpMatrix<Real>::Swap(TpMatrix<Real> *other) {
std::swap(this->data_, other->data_);
std::swap(this->num_rows_, other->num_rows_);
}
template<typename Real>
void TpMatrix<Real>::Cholesky(const SpMatrix<Real> &orig) {
KALDI_ASSERT(orig.NumRows() == this->NumRows());
MatrixIndexT n = this->NumRows();
this->SetZero();
Real *data = this->data_, *jdata = data; // start of j'th row of matrix.
const Real *orig_jdata = orig.Data(); // start of j'th row of matrix.
for (MatrixIndexT j = 0; j < n; j++, jdata += j, orig_jdata += j) {
Real *kdata = data; // start of k'th row of matrix.
Real d(0.0);
for (MatrixIndexT k = 0; k < j; k++, kdata += k) {
Real s = cblas_Xdot(k, kdata, 1, jdata, 1);
// (*this)(j, k) = s = (orig(j, k) - s)/(*this)(k, k);
jdata[k] = s = (orig_jdata[k] - s)/kdata[k];
d = d + s*s;
}
// d = orig(j, j) - d;
d = orig_jdata[j] - d;
if (d >= 0.0) {
// (*this)(j, j) = std::sqrt(d);
jdata[j] = std::sqrt(d);
} else {
KALDI_ERR << "Cholesky decomposition failed. Maybe matrix "
"is not positive definite.";
}
}
}
template<typename Real>
void TpMatrix<Real>::CopyFromMat(const MatrixBase<Real> &M,
MatrixTransposeType Trans) {
if (Trans == kNoTrans) {
KALDI_ASSERT(this->NumRows() == M.NumRows() && M.NumRows() == M.NumCols());
MatrixIndexT D = this->NumRows();
const Real *in_i = M.Data();
MatrixIndexT stride = M.Stride();
Real *out_i = this->data_;
for (MatrixIndexT i = 0; i < D; i++, in_i += stride, out_i += i)
for (MatrixIndexT j = 0; j <= i; j++)
out_i[j] = in_i[j];
} else {
KALDI_ASSERT(this->NumRows() == M.NumRows() && M.NumRows() == M.NumCols());
MatrixIndexT D = this->NumRows();
const Real *in_i = M.Data();
MatrixIndexT stride = M.Stride();
Real *out_i = this->data_;
for (MatrixIndexT i = 0; i < D; i++, in_i++, out_i += i) {
for (MatrixIndexT j = 0; j <= i; j++)
out_i[j] = in_i[stride*j];
}
}
}
template class TpMatrix<float>;
template class TpMatrix<double>;
} // namespace kaldi