mathkit  2.0
 All Classes Namespaces Functions Variables Typedefs
Matrix.cpp
1 #include <sstream>
2 #include <stdexcept>
3 #include <algorithm>
4 #include "Matrix.hpp"
5 
6 using std::abs;
7 using std::sqrt;
8 using std::ostringstream;
9 using std::invalid_argument;
10 
11 namespace mathkit {
12 
13  Matrix::Matrix() : _data(), _dimen(0, 0) {}
14 
15  Matrix::Matrix(const Table & data) {
16  _data = data;
17  _dimen = Dimen(data[0].size(), data.size());
18  }
19 
20  Matrix::Matrix(Dimen dimen, double val) {
21  _data = Table(dimen.second, Vector(dimen.first, val));
22  _dimen = dimen;
23  }
24 
25  Matrix::Matrix(const Vector & data, Dimen dimen, bool byrow) {
26  setData(data, dimen, byrow);
27  }
28 
29  void Matrix::setData(const Table & data) {
30  _data = data;
31  _dimen = Dimen(data[0].size(), data.size());
32  }
33 
34  void Matrix::setData(Dimen dimen, double val) {
35  _data = Table(dimen.second, Vector(dimen.first, val));
36  _dimen = dimen;
37  }
38 
39  void Matrix::setData(const Vector & data, Dimen dimen, bool byrow) {
40  _dimen = dimen;
41  _data = Table(dimen.second);
42  if (byrow) {
43  int k = 0;
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++]);
47  }
48  else {
49  int k = 0;
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++]);
53  }
54  }
55 
57  return _data;
58  }
59 
60  Dimen Matrix::dimen() const {
61  return _dimen;
62  }
63 
64  bool Matrix::dimen(Dimen dimen) {
65  if (!_data.empty()) return false;
66  _dimen = dimen;
67  return true;
68  }
69 
70  int Matrix::nrow() const {
71  return _dimen.first;
72  }
73 
74  int Matrix::ncol() const {
75  return _dimen.second;
76  }
77 
78  Matrix Matrix::round(int ndec) const {
79  Matrix mat = *this;
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);
84  return mat;
85  }
86 
87  bool Matrix::eq(const Matrix & mat, const Epsilon & eps) const {
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;
93  return true;
94  }
95 
96  bool Matrix::square() const {
97  return _dimen.first == _dimen.second;
98  }
99 
100  bool Matrix::symmetric(const Epsilon & eps) const {
101  if (!square()) return false;
102  return eq(transpose(), eps);
103  }
104 
105  bool Matrix::uptri(const Epsilon & eps) const {
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;
111  return true;
112  }
113 
114  bool Matrix::lowtri(const Epsilon & eps) const {
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;
120  return true;
121  }
122 
123  bool Matrix::diag(const Epsilon & eps) const {
124  return uptri(eps) && lowtri(eps);
125  }
126 
127  Vector Matrix::getRow(int m) const {
128  Vector row(_dimen.second);
129  for (int i = 0; i < _dimen.second; ++i) row[i] = _data[i][m];
130  return row;
131  }
132 
133  Vector Matrix::getCol(int n) const {
134  return _data[n];
135  }
136 
137  void Matrix::setRow(int m, const Vector & row) {
138  for (int i = 0; i < _dimen.second; ++i) _data[i][m] = row[i];
139  }
140 
141  void Matrix::setCol(int n, const Vector & col) {
142  _data[n] = col;
143  }
144 
145  void Matrix::exchange(int i, int j, bool row) {
146  Vector temp;
147  if (row) {
148  temp = getRow(i);
149  setRow(i, getRow(j));
150  setRow(j, temp);
151  }
152  else {
153  temp = getCol(i);
154  setCol(i, getCol(j));
155  setCol(j, temp);
156  }
157  }
158 
160  Table tbl;
161  for (int i = 0; i < _dimen.first; ++i) tbl.push_back(getRow(i));
162  return Matrix(tbl);
163  }
164 
165  int Matrix::rank(const Epsilon & eps) const {
166  Matrix mat = *this;
167  int m = mat._dimen.first;
168  int n = mat._dimen.second;
169  int limit = m < n ? m : n;
170  int r = 0;
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))) {
176  mat.exchange(r, j);
177  exchanged = true;
178  break;
179  }
180  }
181  if (!exchanged) continue;
182  }
183  for (int k = r + 1; k < m; ++k) {
184  if (!eps(mat(k, i))) {
185  double s = mat(k, i) / mat(r, i);
186  mat.setRow(k, -s * mat.getRow(r) + mat.getRow(k));
187  }
188  }
189  r += 1;
190  if (r == limit) break;
191  }
192  return r;
193  }
194 
195  Matrix Matrix::submat(int frow, int lrow, int fcol, int lcol) {
196  int m = lrow - frow + 1, n = lcol - fcol + 1;
197  Table tbl(n, Vector(m));
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];
201  return Matrix(tbl);
202  }
203 
204  Matrix Matrix::resid(int i, int j) const {
205  Table tbl;
206  for (int n = 0; n < _dimen.second; ++n) {
207  if (n != j) {
208  Vector col = _data[n];
209  Vector::iterator iter = col.begin();
210  for (int m = 0; m < i; ++m) ++iter;
211  col.erase(iter);
212  tbl.push_back(col);
213  }
214  }
215  return Matrix(tbl);
216  }
217 
218  double & Matrix::operator()(int m, int n) {
219  return _data[n][m];
220  }
221 
222  double Matrix::operator()(int m, int n) const {
223  return _data[n][m];
224  }
225 
226  Matrix Matrix::operator+(const Matrix & mat) const {
227  if (_dimen != mat._dimen) throw invalid_argument("Incompatible matrices.");
228  Table tbl;
229  for (int i = 0; i < _dimen.second; ++i) tbl.push_back(_data[i] + mat._data[i]);
230  return Matrix(tbl);
231  }
232 
233  Matrix Matrix::operator-(const Matrix & mat) const {
234  if (_dimen != mat._dimen) throw invalid_argument("Incompatible matrices.");
235  Table tbl;
236  for (int i = 0; i < _dimen.second; ++i) tbl.push_back(_data[i] - mat._data[i]);
237  return Matrix(tbl);
238  }
239 
240  Matrix Matrix::operator*(const Matrix & mat) const {
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];
246  return Matrix(tbl);
247  }
248 
249  bool Matrix::operator==(const Matrix & mat) const {
250  return _dimen == mat._dimen && _data == mat._data;
251  }
252 
253  bool Matrix::operator!=(const Matrix & mat) const {
254  return !(_dimen == mat._dimen && _data == mat._data);
255  }
256 
257  Matrix operator*(const Matrix & mat, double scalar) {
258  Table tbl;
259  int n = mat.ncol();
260  for (int i = 0; i < n; ++i) tbl.push_back(scalar * mat.getCol(i));
261  return Matrix(tbl);
262  }
263 
264  Matrix operator*(double scalar, const Matrix & mat) {
265  return mat * scalar;
266  }
267 
268  ostream & operator<<(ostream & outs, const Matrix & mat) {
269  save(outs, mat.getData());
270  return outs;
271  }
272 
273  istream & operator>>(istream & ins, Matrix & mat) {
274  if (mat.dimen() != Dimen(0, 0))
275  mat.setData(load(ins, mat.nrow(), mat.ncol()));
276  return ins;
277  }
278 
279  Matrix t(const Matrix & mat) {
280  return mat.transpose();
281  }
282 
283  Matrix eye(int n) {
284  Matrix mat(Dimen(n, n));
285  for (int i = 0; i < n; ++i) mat(i, i) = 1;
286  return mat;
287  }
288 
289  Vector diag(const Matrix & mat) {
290  Vector vec;
291  int n = mat.dimen().second;
292  for (int i = 0; i < n; ++i) vec.push_back(mat(i, i));
293  return vec;
294  }
295 
296  Matrix diag(const Vector & vec) {
297  int n = vec.size();
298  Matrix mat(Dimen(n, n));
299  for (int i = 0; i < n; ++i) mat(i, i) = vec[i];
300  return mat;
301  }
302 
303  double det(const Matrix & mat, const Epsilon & eps) {
304  if (!mat.square()) throw invalid_argument("Square matrix needed.");
305  int n = mat.ncol();
306  if (n == 1) return eps(mat(0, 0)) ? 0 : mat(0, 0);
307  double result = 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;
312  result += resid;
313  }
314  }
315  return result;
316  }
317 
318  Matrix inv(const Matrix & mat, const Epsilon & eps) {
319  if (!mat.square()) throw invalid_argument("Square matrix needed.");
320  int n = mat.ncol();
321  Matrix I = eye(n);
322  Table tbl;
323  Matrices lup = lu(mat, eps);
324  for (int i = 0; i < n; ++i) {
325  Vector result = lusolve(lup, I.getCol(i), eps);
326  if (result.empty()) throw invalid_argument("Singular matrix.");
327  tbl.push_back(result);
328  }
329  return Matrix(tbl);
330  }
331 
332  double trace(const Matrix & mat, const Epsilon & eps) {
333  if (!mat.square()) throw invalid_argument("Square matrix needed.");
334  int n = mat.ncol();
335  double val = 0;
336  for (int i = 0; i < n; ++i)
337  if (!eps(mat(i, i))) val += mat(i, i);
338  return val;
339  }
340 
341  bool positive(const Matrix & mat, const Epsilon & eps) {
342  if (!mat.symmetric(eps)) throw invalid_argument("Symmetric matrix needed.");
343  Matrices lup = lu(mat, eps);
344  if (lup.empty()) return false;
345  const Matrix &U = lup[1];
346  int n = U.ncol();
347  for (int i = 0; i < n; ++i)
348  if (U(i, i) < 0) return false;
349  return true;
350  }
351 
352  Table table(const Vector & data, Dimen dimen, bool byrow) {
353  Table tbl(dimen.second);
354  if (byrow) {
355  int k = 0;
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++]);
359  }
360  else {
361  int k = 0;
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++]);
365  }
366  return tbl;
367  }
368 
369  Matrix rowmat(const Vector & data) {
370  return Matrix(data, Dimen(1, data.size()));
371  }
372 
373  Matrix colmat(const Vector & data) {
374  return Matrix(data, Dimen(data.size(), 1));
375  }
376 
377  Matrix cbind(const Matrix & mat1, const Matrix & mat2) {
378  if (mat1.nrow() != mat2.nrow()) throw invalid_argument("Incompatible matrices.");
379  Table tbl;
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));
382  return Matrix(tbl);
383  }
384 
385  Matrix rbind(const Matrix & mat1, const Matrix & mat2) {
386  int n = mat1.ncol();
387  if (n != mat2.ncol()) throw invalid_argument("Incompatible matrices.");
388  Table tbl(n);
389  for (int i = 0; i < mat1.nrow(); ++i) {
390  Vector row = mat1.getRow(i);
391  for (int j = 0; j < n; ++j) tbl[j].push_back(row[j]);
392  }
393  for (int i = 0; i < mat2.nrow(); ++i) {
394  Vector row = mat2.getRow(i);
395  for (int j = 0; j < n; ++j) tbl[j].push_back(row[j]);
396  }
397  return Matrix(tbl);
398  }
399 
400  Vector eigval(const Matrix & mat, const Epsilon & eps) {
401  if (!mat.square()) throw invalid_argument("Square matrix needed.");
402  if (mat.uptri(eps) || mat.lowtri(eps)) return diag(mat);
403  Matrices qrmat = qr(mat, eps);
404  Matrix m = qrmat[1] * qrmat[0];
405  while (!(m.uptri(eps) || m.lowtri(eps))) {
406  qrmat = qr(m, eps);
407  m = qrmat[1] * qrmat[0];
408  }
409  return diag(m);
410  }
411 
412  Eigen eigen(const Matrix & mat, const Epsilon & eps) {
413  if (!mat.symmetric(eps)) throw invalid_argument("Symmetric matrix needed.");
414  Eigen eig;
415  int n = mat.ncol();
416  eig.vectors = eye(n);
417  if (mat.diag(eps)) {
418  eig.values = diag(mat);
419  return eig;
420  }
421  Matrices qs;
422  Matrices qrmat = qr(mat, eps);
423  qs.push_back(qrmat[0]);
424  Matrix m = qrmat[1] * qrmat[0];
425  while (!m.diag(eps)) {
426  qrmat = qr(m, eps);
427  qs.push_back(qrmat[0]);
428  m = qrmat[1] * qrmat[0];
429  }
430  eig.values = diag(m);
431  Matrices::const_iterator citer;
432  for (citer = qs.begin(); citer != qs.end(); ++citer)
433  eig.vectors = eig.vectors * *citer;
434  vector<int> zerocol;
435  for (int i = 0; i < n; ++i) {
436  Vector col = eig.vectors.getCol(i);
437  if (zero(col, eps)) {
438  zerocol.push_back(i);
439  eig.vectors.setCol(i, Vector(col.size(), 0));
440  continue;
441  }
442  eig.vectors.setCol(i, col / norm(col));
443  }
444  if (!zerocol.empty()) {
445  Matrix nullmat = nullbasis(mat, false, eps);
446  int i = zerocol.size();
447  int j = nullmat.ncol();
448  int k = 0;
449  while (k < i && k < j) {
450  eig.vectors.setCol(zerocol[k], nullmat.getCol(k));
451  k += 1;
452  }
453  }
454  return eig;
455  }
456 
457  Matrices svd(const Matrix & mat, const Epsilon & eps) {
458  int m = mat.nrow(), n = mat.ncol();
459  Matrix aat = mat * t(mat), ata = t(mat) * mat;
460  Eigen eigaat = eigen(aat, eps), eigata = eigen(ata, eps);
461  Matrix W;
462  Vector::iterator iter;
463  if (m > n) {
464  for (iter = eigata.values.begin(); iter != eigata.values.end(); ++iter) *iter = sqrt(*iter);
465  W = rbind(diag(eigata.values), Matrix(Dimen(m - n, n), 0));
466  }
467  else {
468  for (iter = eigaat.values.begin(); iter != eigaat.values.end(); ++iter) *iter = sqrt(*iter);
469  if (m < n) W = cbind(diag(eigaat.values), Matrix(Dimen(m, n - m), 0));
470  }
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);
474  eigaat.vectors.setCol(i, col / norm(col));
475  }
476  Matrices result;
477  result.push_back(eigaat.vectors);
478  result.push_back(eigata.vectors);
479  result.push_back(W);
480  return result;
481  }
482 
483  Elim elim(const Matrix & mat, bool augment, const Epsilon & eps) {
484  Matrix coef_mat = mat;
485  int m = coef_mat.nrow();
486  int n = augment ? coef_mat.ncol() - 1 : coef_mat.ncol();
487  int limit = m < n ? m : n;
488  Elim result;
489  int r = 0;
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))) {
495  coef_mat.exchange(r, j);
496  exchanged = true;
497  break;
498  }
499  }
500  if (!exchanged) continue;
501  }
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);
505  coef_mat.setRow(k, -s * coef_mat.getRow(r) + coef_mat.getRow(k));
506  }
507  }
508  result.pivotcol.push_back(i);
509  r += 1;
510  if (r == limit) break;
511  }
512  for (int i = 0; i < n; ++i)
513  if (!binary_search(result.pivotcol.begin(), result.pivotcol.end(), i))
514  result.freecol.push_back(i);
515  result.coef_mat = coef_mat;
516  return result;
517  }
518 
519  Matrix nullbasis(const Matrix & mat, bool augment, const Epsilon & eps) {
520  Elim elm = elim(mat, augment, eps);
521  if (elm.freecol.empty()) return Matrix();
522  int m = elm.pivotcol.size();
523  int n = augment ? mat.ncol() - 1 : mat.ncol();
524  Matrix pivmat(Dimen(m, n), 0);
525  for (int i = 0; i < m; ++i)
526  for (int j = 0; j < n; ++j)
527  pivmat(i, j) = elm.coef_mat(i, j);
528  Matrix U(Dimen(m, m), 0);
529  for (int i = 0; i < m; ++i) U.setCol(i, pivmat.getCol(elm.pivotcol[i]));
530  int k = elm.freecol.size();
531  Matrix basis(Dimen(n, k), 0);
532  for (int i = 0; i < k; ++i) {
533  basis(elm.freecol[i], i) = 1;
534  Vector b(m, 0);
535  for (int j = 0; j < m; ++j)
536  if (!eps(pivmat(j, elm.freecol[i]))) b[j] = -pivmat(j, elm.freecol[i]);
537  Vector pivar = upsolve(U, b, eps);
538  for (int j = 0; j < m; ++j) basis(elm.pivotcol[j], i) = pivar[j];
539  }
540  basis = qr(basis, eps)[0];
541  return basis;
542  }
543 
544  Matrices lu(const Matrix & mat, const Epsilon & eps) {
545  if (!mat.square()) throw invalid_argument("Square matrix needed.");
546  int n = mat.ncol();
547  Matrix L = eye(n), P = eye(n);
548  Matrix U = mat;
549  for (int i = 0; i < n; ++i) {
550  if (eps(U(i, i))) {
551  bool exchanged = false;
552  for (int k = i + 1; k < n; ++k) {
553  if (!eps(U(k, i))) {
554  U.exchange(i, k);
555  P.exchange(i, k);
556  for (int j = 0; j < i; ++j) {
557  double temp = L(i, j);
558  L(i, j) = L(k, j);
559  L(k, j) = temp;
560  }
561  exchanged = true;
562  break;
563  }
564  }
565  if (!exchanged) return Matrices();
566  }
567  for (int j = i + 1; j < n; ++j) {
568  if (!eps(U(j, i))) {
569  double s = U(j, i) / U(i, i);
570  L(j, i) = s;
571  U.setRow(j, -s * U.getRow(i) + U.getRow(j));
572  }
573  }
574  }
575  Matrices result;
576  result.push_back(L);
577  result.push_back(U);
578  result.push_back(P);
579  return result;
580  }
581 
582  Matrices qr(const Matrix & mat, const Epsilon & eps) {
583  int n = mat.ncol();
584  Table tbl;
585  for (int i = 0; i < n; ++i) {
586  Vector vec = mat.getCol(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));
590  continue;
591  }
592  tbl.push_back(vec / norm(vec));
593  }
594  Matrix Q(tbl);
595  Matrix R(Dimen(n, n));
596  for (int i = 0; i < n; ++i)
597  for (int j = i; j < n; ++j)
598  R(i, j) = dot_prod(tbl[i], mat.getCol(j));
599  Matrices result;
600  result.push_back(Q);
601  result.push_back(R);
602  return result;
603  }
604 
605  Vector linsolve(const Matrix & mat, const Vector & vec, const Epsilon & eps) {
606  return lusolve(lu(mat, eps), vec, eps);
607  }
608 
609  Vector lusolve(const Matrices & lup, const Vector & vec, const Epsilon & eps) {
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.");
614  Matrix b = colmat(vec);
615  if (P != eye(n)) b = P * b;
616  Vector result;
617  for (unsigned int i = 0; i < n; ++i) {
618  double val = 0;
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);
622  }
623  b = colmat(result);
624  for (int i = n - 1; i >= 0; --i) {
625  double val = 0;
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);
629  }
630  return result;
631  }
632 
633  Vector qrsolve(const Matrices & qrmat, const Vector & vec, const Epsilon & eps) {
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.");
639  Matrix b = colmat(vec);
640  b = t(Q) * b;
641  int n = R.nrow();
642  Vector result(n);
643  for (int i = n - 1; i >= 0; --i) {
644  double val = 0;
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);
648  }
649  return result;
650  }
651 
652  Vector upsolve(const Matrix & upmat, const Vector & vec, const Epsilon & eps) {
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.");
656  Vector result(n, 0);
657  Matrix b = colmat(vec);
658  for (int i = n - 1; i >= 0; --i) {
659  double val = 0;
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);
663  }
664  return result;
665  }
666 
667  string to_str(const Matrix & mat, bool byrow, string delim) {
668  ostringstream oss;
669  if (byrow) {
670  int m = mat.nrow();
671  for (int i = 0; i < m; ++i) oss << to_str(mat.getRow(i));
672  }
673  else {
674  int n = mat.ncol();
675  for (int i = 0; i < n; ++i) oss << to_str(mat.getCol(i));
676  }
677  string result = oss.str();
678  if (delim.size() == 2) result = delim[0] + result + delim[1];
679  return result;
680  }
681 }