6 using std::ostringstream;
7 using std::invalid_argument;
15 _dimen =
Dimen(data[0].size(), data.size());
19 _data =
Table(dimen.second,
Vector(dimen.first, val));
29 _dimen =
Dimen(data[0].size(), data.size());
33 _data =
Table(dimen.second,
Vector(dimen.first, val));
39 _data =
Table(dimen.second);
42 for (
int i = 0; i < dimen.first; ++i)
43 for (
int j = 0; j < dimen.second; ++j)
44 _data[j].push_back(data[k++]);
48 for (
int i = 0; i < dimen.second; ++i)
49 for (
int j = 0; j < dimen.first; ++j)
50 _data[i].push_back(data[k++]);
63 if (!_data.empty())
return false;
78 int m = mat._dimen.first, n = mat._dimen.second;
79 for (
int i = 0; i < m; ++i)
80 for (
int j = 0; j < n; ++j)
81 mat(i, j) =
prec(mat(i, j), ndec);
86 return _dimen.first == _dimen.second;
90 if (!
square())
return false;
96 for (
int i = 0; i < _dimen.second; ++i) row[i] = _data[i][m];
105 for (
int i = 0; i < _dimen.second; ++i) _data[i][m] = row[i];
128 for (
int i = 0; i < _dimen.first; ++i) tbl.push_back(
getRow(i));
134 int m = mat._dimen.first;
135 int n = mat._dimen.second;
136 int limit = m < n ? m : n;
138 for (
int i = 0; i < n; ++i) {
139 if (eps(mat(r, i))) {
140 bool exchanged =
false;
141 for (
int j = r + 1; j < m; ++j) {
142 if (!eps(mat(j, i))) {
148 if (!exchanged)
continue;
150 for (
int k = r + 1; k < m; ++k) {
151 if (!eps(mat(k, i))) {
152 double s = mat(k, i) / mat(r, i);
157 if (r == limit)
break;
163 int m = lrow - frow + 1, n = lcol - fcol + 1;
165 for (
int i = 0; i < m; ++i)
166 for (
int j = 0; j < n; ++j)
167 tbl[j][i] = _data[fcol + j][frow + i];
173 for (
int n = 0; n < _dimen.second; ++n) {
176 Vector::iterator iter = col.begin();
177 for (
int m = 0; m < i; ++m) ++iter;
194 if (_dimen != mat._dimen)
throw invalid_argument(
"Incompatible matrices.");
196 for (
int i = 0; i < _dimen.second; ++i) tbl.push_back(_data[i] + mat._data[i]);
201 if (_dimen != mat._dimen)
throw invalid_argument(
"Incompatible matrices.");
203 for (
int i = 0; i < _dimen.second; ++i) tbl.push_back(_data[i] - mat._data[i]);
208 if (_dimen.second != mat._dimen.first)
throw invalid_argument(
"Incompatible matrices.");
209 Table tbl(mat._dimen.second,
Vector(_dimen.first, 0));
210 for (
int i = 0; i < mat._dimen.second; ++i)
211 for (
int j = 0; j < _dimen.first; ++j)
212 tbl[i] = tbl[i] + mat._data[i][j] * _data[j];
217 return _dimen == mat._dimen && _data == mat._data;
221 return !(_dimen == mat._dimen && _data == mat._data);
227 for (
int i = 0; i < n; ++i) tbl.push_back(scalar * mat.
getCol(i));
252 for (
int i = 0; i < n; ++i) mat(i, i) = 1;
258 int n = mat.
dimen().second;
259 for (
int i = 0; i < n; ++i) vec.push_back(mat(i, i));
266 for (
int i = 0; i < n; ++i) mat(i, i) = vec[i];
271 if (!mat.
square())
throw invalid_argument(
"Square matrix needed.");
273 if (n == 1)
return eps(mat(0, 0)) ? 0 : mat(0, 0);
275 for (
int i = 0; i < n; ++i) {
276 if (!eps(mat(0, i))) {
277 double resid = mat(0, i) *
det(mat.
resid(0, i));
278 if (i % 2 != 0) resid = -resid;
286 if (!mat.
square())
throw invalid_argument(
"Square matrix needed.");
291 for (
int i = 0; i < n; ++i) {
293 if (result.empty())
throw invalid_argument(
"Singular matrix.");
294 tbl.push_back(result);
300 if (!mat.
square())
throw invalid_argument(
"Square matrix needed.");
303 for (
int i = 0; i < n; ++i)
304 if (!eps(mat(i, i))) val += mat(i, i);
309 if (!mat.
symmetric())
throw invalid_argument(
"Symmetric matrix needed.");
311 if (lup.empty())
return false;
314 for (
int i = 0; i < n; ++i)
315 if (U(i, i) < 0)
return false;
320 Table tbl(dimen.second);
323 for (
int i = 0; i < dimen.first; ++i)
324 for (
int j = 0; j < dimen.second; ++j)
325 tbl[j].push_back(data[k++]);
329 for (
int i = 0; i < dimen.second; ++i)
330 for (
int j = 0; j < dimen.first; ++j)
331 tbl[i].push_back(data[k++]);
345 if (mat1.
nrow() != mat2.
nrow())
throw invalid_argument(
"Incompatible matrices.");
347 for (
int i = 0; i < mat1.
ncol(); ++i) tbl.push_back(mat1.
getCol(i));
348 for (
int i = 0; i < mat2.
ncol(); ++i) tbl.push_back(mat2.
getCol(i));
354 if (n != mat2.
ncol())
throw invalid_argument(
"Incompatible matrices.");
356 for (
int i = 0; i < mat1.
nrow(); ++i) {
358 for (
int j = 0; j < n; ++j) tbl[j].push_back(row[j]);
360 for (
int i = 0; i < mat2.
nrow(); ++i) {
362 for (
int j = 0; j < n; ++j) tbl[j].push_back(row[j]);
368 if (!mat.
square())
throw invalid_argument(
"Square matrix needed.");
372 for (
int i = 0; i < n; ++i) {
374 bool exchanged =
false;
375 for (
int k = i + 1; k < n; ++k) {
379 for (
int j = 0; j < i; ++j) {
380 double temp = L(i, j);
390 for (
int j = i + 1; j < n; ++j) {
392 double s = U(j, i) / U(i, i);
409 for (
int i = 0; i < n; ++i) {
411 for (
int j = 0; j < i; ++j) vec = vec - tbl[j] *
dot_prod(tbl[j], vec);
412 tbl.push_back(vec /
norm(vec));
416 for (
int i = 0; i < n; ++i)
417 for (
int j = i; j < n; ++j)
430 if (lup.empty())
return Vector();
431 const Matrix &L = lup[0], &U = lup[1], &P = lup[2];
432 unsigned int n = U.
nrow();
433 if (n != vec.size())
throw invalid_argument(
"Incompatible arguments.");
435 if (P !=
eye(n)) b = P * b;
437 for (
unsigned int i = 0; i < n; ++i) {
439 for (
unsigned int j = 0; j < i; ++j)
440 if (!eps(L(i, j)))val += L(i, j) * result[j];
441 result.push_back(b(i, 0) - val);
444 for (
int i = n - 1; i >= 0; --i) {
446 for (
int j = n - 1; j > i; --j)
447 if (!eps(U(i, j))) val += U(i, j) * result[j];
448 result[i] = (b(i, 0) - val) / U(i, i);
454 if (qrmat.empty())
return Vector();
455 const Matrix &Q = qrmat[0], &R = qrmat[1];
456 unsigned int m = Q.
nrow();
457 if (vec.size() != m)
throw invalid_argument(
"Incompatible arguments.");
462 for (
int i = n - 1; i >= 0; --i) {
464 for (
int j = n - 1; j > i; --j)
465 if (!eps(R(i, j))) val += R(i, j) * result[j];
466 result[i] = (b(i, 0) - val) / R(i, i);
475 for (
int i = 0; i < m; ++i) oss <<
to_str(mat.
getRow(i));
479 for (
int i = 0; i < n; ++i) oss <<
to_str(mat.
getCol(i));
481 string result = oss.str();
482 if (delim.size() == 2) result = delim[0] + result + delim[1];