8 using std::ostringstream;
9 using std::invalid_argument;
17 _dimen =
Dimen(data[0].size(), data.size());
21 _data =
Table(dimen.second,
Vector(dimen.first, val));
31 _dimen =
Dimen(data[0].size(), data.size());
35 _data =
Table(dimen.second,
Vector(dimen.first, val));
41 _data =
Table(dimen.second);
44 for (
int i = 0; i < dimen.first; ++i)
45 for (
int j = 0; j < dimen.second; ++j)
46 _data[j].push_back(data[k++]);
50 for (
int i = 0; i < dimen.second; ++i)
51 for (
int j = 0; j < dimen.first; ++j)
52 _data[i].push_back(data[k++]);
65 if (!_data.empty())
return false;
80 int m = mat._dimen.first, n = mat._dimen.second;
81 for (
int i = 0; i < m; ++i)
82 for (
int j = 0; j < n; ++j)
83 mat(i, j) =
prec(mat(i, j), ndec);
88 if (_dimen != mat._dimen)
return false;
89 int m = _dimen.first, n = _dimen.second;
90 for (
int i = 0; i < m; ++i)
91 for (
int j = 0; j < n; ++j)
92 if (!eps(_data[j][i] - mat(i, j)))
return false;
97 return _dimen.first == _dimen.second;
101 if (!
square())
return false;
106 if (!
square())
return false;
107 int n = _dimen.second;
108 for (
int i = 0; i < n; ++i)
109 for (
int j = i + 1; j < n; ++j)
110 if (!eps(_data[i][j]))
return false;
115 if (!
square())
return false;
116 int n = _dimen.first;
117 for (
int i = 0; i < n; ++i)
118 for (
int j = i + 1; j < n; ++j)
119 if (!eps(_data[j][i]))
return false;
128 Vector row(_dimen.second);
129 for (
int i = 0; i < _dimen.second; ++i) row[i] = _data[i][m];
138 for (
int i = 0; i < _dimen.second; ++i) _data[i][m] = row[i];
161 for (
int i = 0; i < _dimen.first; ++i) tbl.push_back(
getRow(i));
167 int m = mat._dimen.first;
168 int n = mat._dimen.second;
169 int limit = m < n ? m : n;
171 for (
int i = 0; i < n; ++i) {
172 if (eps(mat(r, i))) {
173 bool exchanged =
false;
174 for (
int j = r + 1; j < m; ++j) {
175 if (!eps(mat(j, i))) {
181 if (!exchanged)
continue;
183 for (
int k = r + 1; k < m; ++k) {
184 if (!eps(mat(k, i))) {
185 double s = mat(k, i) / mat(r, i);
190 if (r == limit)
break;
196 int m = lrow - frow + 1, n = lcol - fcol + 1;
198 for (
int i = 0; i < m; ++i)
199 for (
int j = 0; j < n; ++j)
200 tbl[j][i] = _data[fcol + j][frow + i];
206 for (
int n = 0; n < _dimen.second; ++n) {
209 Vector::iterator iter = col.begin();
210 for (
int m = 0; m < i; ++m) ++iter;
227 if (_dimen != mat._dimen)
throw invalid_argument(
"Incompatible matrices.");
229 for (
int i = 0; i < _dimen.second; ++i) tbl.push_back(_data[i] + mat._data[i]);
234 if (_dimen != mat._dimen)
throw invalid_argument(
"Incompatible matrices.");
236 for (
int i = 0; i < _dimen.second; ++i) tbl.push_back(_data[i] - mat._data[i]);
241 if (_dimen.second != mat._dimen.first)
throw invalid_argument(
"Incompatible matrices.");
242 Table tbl(mat._dimen.second,
Vector(_dimen.first, 0));
243 for (
int i = 0; i < mat._dimen.second; ++i)
244 for (
int j = 0; j < _dimen.second; ++j)
245 tbl[i] = tbl[i] + mat._data[i][j] * _data[j];
250 return _dimen == mat._dimen && _data == mat._data;
254 return !(_dimen == mat._dimen && _data == mat._data);
260 for (
int i = 0; i < n; ++i) tbl.push_back(scalar * mat.
getCol(i));
285 for (
int i = 0; i < n; ++i) mat(i, i) = 1;
291 int n = mat.
dimen().second;
292 for (
int i = 0; i < n; ++i) vec.push_back(mat(i, i));
299 for (
int i = 0; i < n; ++i) mat(i, i) = vec[i];
304 if (!mat.
square())
throw invalid_argument(
"Square matrix needed.");
306 if (n == 1)
return eps(mat(0, 0)) ? 0 : mat(0, 0);
308 for (
int i = 0; i < n; ++i) {
309 if (!eps(mat(0, i))) {
310 double resid = mat(0, i) *
det(mat.
resid(0, i));
311 if (i % 2 != 0) resid = -resid;
319 if (!mat.
square())
throw invalid_argument(
"Square matrix needed.");
324 for (
int i = 0; i < n; ++i) {
326 if (result.empty())
throw invalid_argument(
"Singular matrix.");
327 tbl.push_back(result);
333 if (!mat.
square())
throw invalid_argument(
"Square matrix needed.");
336 for (
int i = 0; i < n; ++i)
337 if (!eps(mat(i, i))) val += mat(i, i);
342 if (!mat.
symmetric(eps))
throw invalid_argument(
"Symmetric matrix needed.");
344 if (lup.empty())
return false;
347 for (
int i = 0; i < n; ++i)
348 if (U(i, i) < 0)
return false;
353 Table tbl(dimen.second);
356 for (
int i = 0; i < dimen.first; ++i)
357 for (
int j = 0; j < dimen.second; ++j)
358 tbl[j].push_back(data[k++]);
362 for (
int i = 0; i < dimen.second; ++i)
363 for (
int j = 0; j < dimen.first; ++j)
364 tbl[i].push_back(data[k++]);
378 if (mat1.
nrow() != mat2.
nrow())
throw invalid_argument(
"Incompatible matrices.");
380 for (
int i = 0; i < mat1.
ncol(); ++i) tbl.push_back(mat1.
getCol(i));
381 for (
int i = 0; i < mat2.
ncol(); ++i) tbl.push_back(mat2.
getCol(i));
387 if (n != mat2.
ncol())
throw invalid_argument(
"Incompatible matrices.");
389 for (
int i = 0; i < mat1.
nrow(); ++i) {
391 for (
int j = 0; j < n; ++j) tbl[j].push_back(row[j]);
393 for (
int i = 0; i < mat2.
nrow(); ++i) {
395 for (
int j = 0; j < n; ++j) tbl[j].push_back(row[j]);
401 if (!mat.
square())
throw invalid_argument(
"Square matrix needed.");
404 Matrix m = qrmat[1] * qrmat[0];
407 m = qrmat[1] * qrmat[0];
413 if (!mat.
symmetric(eps))
throw invalid_argument(
"Symmetric matrix needed.");
423 qs.push_back(qrmat[0]);
424 Matrix m = qrmat[1] * qrmat[0];
425 while (!m.
diag(eps)) {
427 qs.push_back(qrmat[0]);
428 m = qrmat[1] * qrmat[0];
431 Matrices::const_iterator citer;
432 for (citer = qs.begin(); citer != qs.end(); ++citer)
435 for (
int i = 0; i < n; ++i) {
437 if (
zero(col, eps)) {
438 zerocol.push_back(i);
444 if (!zerocol.empty()) {
446 int i = zerocol.size();
447 int j = nullmat.
ncol();
449 while (k < i && k < j) {
459 Matrix aat = mat *
t(mat), ata =
t(mat) * mat;
462 Vector::iterator iter;
464 for (iter = eigata.values.begin(); iter != eigata.values.end(); ++iter) *iter = sqrt(*iter);
468 for (iter = eigaat.
values.begin(); iter != eigaat.
values.end(); ++iter) *iter = sqrt(*iter);
471 int k = m < n ? m : n;
472 for (
int i = 0; i < k; ++i) {
473 Vector col = (mat *
colmat(eigata.vectors.getCol(i))).getCol(0);
477 result.push_back(eigaat.
vectors);
478 result.push_back(eigata.vectors);
485 int m = coef_mat.
nrow();
486 int n = augment ? coef_mat.
ncol() - 1 : coef_mat.
ncol();
487 int limit = m < n ? m : n;
490 for (
int i = 0; i < n; ++i) {
491 if (eps(coef_mat(r, i))) {
492 bool exchanged =
false;
493 for (
int j = r + 1; j < m; ++j) {
494 if (!eps(coef_mat(j, i))) {
500 if (!exchanged)
continue;
502 for (
int k = r + 1; k < m; ++k) {
503 if (!eps(coef_mat(k, i))) {
504 double s = coef_mat(k, i) / coef_mat(r, i);
510 if (r == limit)
break;
512 for (
int i = 0; i < n; ++i)
523 int n = augment ? mat.
ncol() - 1 : mat.
ncol();
525 for (
int i = 0; i < m; ++i)
526 for (
int j = 0; j < n; ++j)
532 for (
int i = 0; i < k; ++i) {
535 for (
int j = 0; j < m; ++j)
536 if (!eps(pivmat(j, elm.
freecol[i]))) b[j] = -pivmat(j, elm.
freecol[i]);
538 for (
int j = 0; j < m; ++j) basis(elm.
pivotcol[j], i) = pivar[j];
540 basis =
qr(basis, eps)[0];
545 if (!mat.
square())
throw invalid_argument(
"Square matrix needed.");
549 for (
int i = 0; i < n; ++i) {
551 bool exchanged =
false;
552 for (
int k = i + 1; k < n; ++k) {
556 for (
int j = 0; j < i; ++j) {
557 double temp = L(i, j);
567 for (
int j = i + 1; j < n; ++j) {
569 double s = U(j, i) / U(i, i);
585 for (
int i = 0; i < n; ++i) {
587 for (
int j = 0; j < i; ++j) vec = vec - tbl[j] *
dot_prod(tbl[j], vec);
588 if (
zero(vec, eps)) {
589 tbl.push_back(
Vector(vec.size(), 0));
592 tbl.push_back(vec /
norm(vec));
596 for (
int i = 0; i < n; ++i)
597 for (
int j = i; j < n; ++j)
610 if (lup.empty())
return Vector();
611 const Matrix &L = lup[0], &U = lup[1], &P = lup[2];
612 unsigned int n = U.
nrow();
613 if (n != vec.size())
throw invalid_argument(
"Incompatible arguments.");
615 if (P !=
eye(n)) b = P * b;
617 for (
unsigned int i = 0; i < n; ++i) {
619 for (
unsigned int j = 0; j < i; ++j)
620 if (!eps(L(i, j)))val += L(i, j) * result[j];
621 result.push_back(b(i, 0) - val);
624 for (
int i = n - 1; i >= 0; --i) {
626 for (
int j = n - 1; j > i; --j)
627 if (!eps(U(i, j))) val += U(i, j) * result[j];
628 result[i] = (b(i, 0) - val) / U(i, i);
634 const Matrix &Q = qrmat[0], &R = qrmat[1];
635 unsigned int k = Q.
ncol();
636 for (
unsigned int i = 0; i < k; ++i)
if (
zero(Q.
getCol(i), eps))
return Vector();
637 unsigned int m = Q.
nrow();
638 if (vec.size() != m)
throw invalid_argument(
"Incompatible arguments.");
643 for (
int i = n - 1; i >= 0; --i) {
645 for (
int j = n - 1; j > i; --j)
646 if (!eps(R(i, j))) val += R(i, j) * result[j];
647 result[i] = (b(i, 0) - val) / R(i, i);
653 if (!upmat.
uptri(eps))
throw invalid_argument(
"The coefficient matrix should be upper triangular.");
654 unsigned int n = upmat.
nrow();
655 if (n != vec.size())
throw invalid_argument(
"Incompatible arguments.");
658 for (
int i = n - 1; i >= 0; --i) {
660 for (
int j = n - 1; j > i; --j)
661 if (!eps(upmat(i, j))) val += upmat(i, j) * result[j];
662 result[i] = (b(i, 0) - val) / upmat(i, i);
671 for (
int i = 0; i < m; ++i) oss <<
to_str(mat.
getRow(i));
675 for (
int i = 0; i < n; ++i) oss <<
to_str(mat.
getCol(i));
677 string result = oss.str();
678 if (delim.size() == 2) result = delim[0] + result + delim[1];