1 #ifndef __MATRIX_UTILITIES_H__
2 #define __MATRIX_UTILITIES_H__
44 typedef int64_t msize_t;
50 template <
typename T1,
typename T2,
typename T3,
typename A>
51 static void multiply(
const vector<vector<T1> > &left,
const vector<vector<T2> > &right, vector<vector<T3> > &result);
56 template <
typename T1,
typename T2,
typename T3>
57 static void multiply(
const vector<vector<T1> > &left,
const T2 right, vector<vector<T3> > &result);
63 static void rref(vector<vector<T> > &inout);
69 static void inverse(
const vector<vector<T> > &in, vector<vector<T> > &result);
74 template <
typename T1,
typename T2,
typename T3>
75 static void add(
const vector<vector<T1> > &left,
const vector<vector<T2> > &right, vector<vector<T3> > &result);
80 template <
typename T1,
typename T2,
typename T3>
81 static void add(
const vector<vector<T1> > &left,
const T2 right, vector<vector<T3> > &result);
86 template <
typename T1,
typename T2,
typename T3>
87 static void subtract(
const vector<vector<T1> > &left,
const vector<vector<T2> > &right, vector<vector<T3> > &result);
93 static void transpose(
const vector<vector<T> > &in, vector<vector<T> > &result);
99 static bool checkDim(
const vector<vector<T> > &in);
104 template <
typename T>
105 static void resize(
const msize_t rows,
const msize_t columns, vector<vector<T> > &result,
bool destructive =
false);
110 template <
typename T>
111 static void zeros(
const msize_t rows,
const msize_t columns, vector<vector<T> > &result);
116 template <
typename T>
117 static void ones(
const msize_t rows,
const msize_t columns, vector<vector<T> > &result);
122 template <
typename T>
123 static void identity(
const msize_t size, vector<vector<T> > &result);
128 template <
typename T1,
typename T2,
typename T3>
129 static void horizCat(
const vector<vector<T1> > &left,
const vector<vector<T2> > &right, vector<vector<T3> > &result);
134 template <
typename T1,
typename T2,
typename T3>
135 static void vertCat(
const vector<vector<T1> > &top,
const vector<vector<T2> > &bottom, vector<vector<T3> > &result);
140 template <
typename T>
141 static void getChunk(
const msize_t firstrow,
const msize_t lastrow,
const msize_t firstcol,
const msize_t lastcol,
const vector<vector<T> > &in, vector<vector<T> > &result);
147 template <
typename T>
148 static void rref_big(vector<vector<T> > &inout);
151 template <
typename T1,
typename T2,
typename T3,
typename A>
152 void MatrixFunctions::multiply(
const vector<vector<T1> >& left,
const vector<vector<T2> >& right, vector<vector<T3> >& result)
154 msize_t leftrows = (msize_t)left.size(), rightrows = (msize_t)right.size(), leftcols, rightcols;
155 vector<vector<T3> > tempstorage, *tresult = &result;
156 bool copyout =
false;
157 if (&left == &result || &right == &result)
160 tresult = &tempstorage;
162 if (leftrows && rightrows)
164 leftcols = (msize_t)left[0].size();
165 rightcols = (msize_t)right[0].size();
166 if (leftcols && rightcols && (rightrows == leftcols))
168 resize(leftrows, rightcols, (*tresult),
true);
170 for (i = 0; i < leftrows; ++i)
172 for (j = 0; j < rightcols; ++j)
175 for (k = 0; k < leftcols; ++k)
177 accum += left[i][k] * right[k][j];
179 (*tresult)[i][j] = accum;
192 result = tempstorage;
196 template <
typename T1,
typename T2,
typename T3>
197 void MatrixFunctions::multiply(
const vector<vector<T1> > &left,
const T2 right, vector<vector<T3> > &result)
199 msize_t leftrows = (msize_t)left.size(), leftcols;
200 bool doresize =
true;
201 if (&left == &result)
207 leftcols = (msize_t)left[0].size();
210 if (doresize) resize(leftrows, leftcols, result,
true);
212 for (i = 0; i < leftrows; ++i)
214 for (j = 0; j < leftcols; ++j)
216 result[i][j] = left[i][j] * right;
230 void MatrixFunctions::rref_big(vector<vector<T> > &inout)
232 msize_t rows = (msize_t)inout.size(), cols;
235 cols = (msize_t)inout[0].size();
238 vector<msize_t> pivots(rows, -1), missingPivots;
239 msize_t i, j, k, myrow = 0;
242 for (i = 0; i < cols; ++i)
244 if (myrow >= rows)
break;
247 for (j = myrow; j < rows; ++j)
249 if (abs(inout[j][i]) > tempval)
251 pivotrow = (msize_t)j;
252 tempval = abs(inout[j][i]);
257 missingPivots.push_back(i);
260 inout[pivotrow].swap(inout[myrow]);
262 tempval = inout[myrow][i];
263 inout[myrow][i] = (T)1;
264 for (j = i + 1; j < cols; ++j)
266 inout[myrow][j] /= tempval;
268 for (j = myrow + 1; j < rows; ++j)
270 tempval = inout[j][i];
272 for (k = i + 1; k < cols; ++k)
274 inout[j][k] -= tempval * inout[myrow][k];
279 msize_t numMissing = (msize_t)missingPivots.size();
282 msize_t lastPivotCol = pivots[myrow - 1];
283 for (i = myrow - 1; i > 0; --i)
285 msize_t pivotCol = pivots[i];
286 for (j = i - 1; j >= 0; --j)
288 tempval = inout[j][pivotCol];
289 inout[j][pivotCol] = (T)0;
290 for (k = numMissing - 1; k >= 0; --k)
292 msize_t missingCol = missingPivots[k];
293 if (missingCol <= pivotCol)
break;
294 inout[j][missingCol] -= tempval * inout[i][missingCol];
296 for (k = lastPivotCol + 1; k < cols; ++k)
298 inout[j][k] -= tempval * inout[i][k];
314 void MatrixFunctions::rref(vector<vector<T> > &inout)
316 msize_t rows = (msize_t)inout.size(), cols;
319 cols = (msize_t)inout[0].size();
322 if (rows > 7 || cols > 7)
327 msize_t i, j, k, myrow = 0;
330 for (i = 0; i < cols; ++i)
332 if (myrow >= rows)
break;
335 for (j = myrow; j < rows; ++j)
337 if (abs(inout[j][i]) > tempval)
339 pivotrow = (msize_t)j;
340 tempval = abs(inout[j][i]);
347 inout[pivotrow].swap(inout[myrow]);
348 tempval = inout[myrow][i];
350 for (j = i + 1; j < cols; ++j)
352 inout[myrow][j] /= tempval;
354 for (j = 0; j < myrow; ++j)
356 tempval = inout[j][i];
358 for (k = i + 1; k < cols; ++k)
360 inout[j][k] -= tempval * inout[myrow][k];
363 for (j = myrow + 1; j < rows; ++j)
365 tempval = inout[j][i];
367 for (k = i + 1; k < cols; ++k)
369 inout[j][k] -= tempval * inout[myrow][k];
385 void MatrixFunctions::inverse(
const vector<vector<T> > &in, vector<vector<T> > &result)
387 msize_t inrows = (msize_t)in.size(), incols;
390 incols = (msize_t)in[0].size();
391 if (incols == inrows)
393 vector<vector<T> > inter, inter2;
394 identity(incols, inter2);
395 horizCat(in, inter2, inter);
397 getChunk(0, inrows, incols, incols * 2, inter, result);
408 template <
typename T1,
typename T2,
typename T3>
409 void MatrixFunctions::add(
const vector<vector<T1> >& left,
const vector<vector<T2> >& right, vector<vector<T3> >& result)
411 msize_t inrows = (msize_t)left.size(), incols;
412 bool doresize =
true;
413 if (&left == &result || &right == &result)
419 incols = (msize_t)left[0].size();
420 if (inrows == (msize_t)right.size() && incols == (msize_t)right[0].size())
422 if (doresize) resize(inrows, incols, result,
true);
423 for (msize_t i = 0; i < inrows; ++i)
425 for (msize_t j = 0; j < incols; ++j)
427 result[i][j] = left[i][j] + right[i][j];
440 template <
typename T1,
typename T2,
typename T3>
441 void MatrixFunctions::add(
const vector<vector<T1> >& left,
const T2 right, vector<vector<T3> >& result)
443 msize_t inrows = (msize_t)left.size(), incols;
444 bool doresize =
true;
445 if (&left == &result)
451 incols = (msize_t)left[0].size();
452 if (doresize) resize(inrows, incols, result,
true);
453 for (msize_t i = 0; i < inrows; ++i)
455 for (msize_t j = 0; j < incols; ++j)
457 result[i][j] = left[i][j] + right;
466 template <
typename T1,
typename T2,
typename T3>
467 void MatrixFunctions::subtract(
const vector<vector<T1> >& left,
const vector<vector<T2> >& right, vector<vector<T3> >& result)
469 msize_t inrows = (msize_t)left.size(), incols;
470 bool doresize =
true;
471 if (&left == &result || &right == &result)
477 incols = (msize_t)left[0].size();
478 if (inrows == (msize_t)right.size() && incols == (msize_t)right[0].size())
480 if (doresize) resize(inrows, incols, result,
true);
481 for (msize_t i = 0; i < inrows; ++i)
483 for (msize_t j = 0; j < incols; ++j)
485 result[i][j] = left[i][j] - right[i][j];
499 void MatrixFunctions::transpose(
const vector<vector<T> > &in, vector<vector<T> > &result)
501 msize_t inrows = (msize_t)in.size(), incols;
502 vector<vector<T> > tempstorage, *tresult = &result;
503 bool copyout =
false;
507 tresult = &tempstorage;
511 incols = (msize_t)in[0].size();
512 resize(incols, inrows, (*tresult),
true);
513 for (msize_t i = 0; i < inrows; ++i)
515 for (msize_t j = 0; j < incols; ++j)
517 (*tresult)[j][i] = in[i][j];
525 result = tempstorage;
530 bool MatrixFunctions::checkDim(
const vector<vector<T> > &in)
533 msize_t rows = (msize_t)in.size(), columns;
536 columns = (msize_t)in[0].size();
537 for (msize_t i = 1; i < rows; ++i)
539 if (in[i].size() != columns)
549 void MatrixFunctions::resize(
const msize_t rows,
const msize_t columns, vector<vector<T> >& result,
bool destructive)
551 if (destructive && result.size() && ((msize_t)result.capacity() < rows || (msize_t)result[0].capacity() < columns))
556 for (msize_t i = 0; i < (
const msize_t)rows; ++i)
558 result[i].resize(columns);
563 void MatrixFunctions::zeros(
const msize_t rows,
const msize_t columns, vector<vector<T> >& result)
565 resize(rows, columns, result,
true);
566 for (msize_t i = 0; i < rows; ++i)
568 for (msize_t j = 0; j < columns; ++j)
576 void MatrixFunctions::ones(
const msize_t rows,
const msize_t columns, vector<vector<T> >& result)
578 resize(rows, columns, result,
true);
579 for (msize_t i = 0; i < rows; ++i)
581 for (msize_t j = 0; j < columns; ++j)
589 void MatrixFunctions::identity(
const msize_t size, vector<vector<T> >& result)
591 resize(size, size, result,
true);
592 for (msize_t i = 0; i < (
const msize_t)size; ++i)
594 for (msize_t j = 0; j < (
const msize_t)size; ++j)
596 result[i][j] = ((i == j) ? 1 : 0);
601 template <
typename T1,
typename T2,
typename T3>
602 void MatrixFunctions::horizCat(
const vector<vector<T1> >& left,
const vector<vector<T2> >& right, vector<vector<T3> >& result)
604 msize_t inrows = (msize_t)left.size(), leftcols, rightcols;
605 vector<vector<T3> > tempstorage, *tresult = &result;
606 bool copyout =
false;
607 if (&left == &result || &right == &result)
610 tresult = &tempstorage;
612 if (inrows && inrows == (msize_t)right.size())
614 leftcols = (msize_t)left[0].size();
615 rightcols = (msize_t)right[0].size();
617 resize(inrows, leftcols + rightcols, (*tresult));
618 for (msize_t i = 0; i < inrows; ++i)
620 for (msize_t j = 0; j < rightcols; ++j)
622 (*tresult)[i][j + leftcols] = right[i][j];
631 result = tempstorage;
635 template <
typename T1,
typename T2,
typename T3>
636 void MatrixFunctions::vertCat(
const vector<vector<T1> >& top,
const vector<vector<T2> >& bottom, vector<vector<T3> >& result)
638 msize_t toprows = (msize_t)top.size(), botrows = (msize_t)bottom.size(), incols;
639 vector<vector<T3> > tempstorage, *tresult = &result;
640 bool copyout =
false;
641 if (&top == &result || &bottom == &result)
644 tresult = &tempstorage;
646 if (toprows && botrows)
648 incols = (msize_t)top[0].size();
649 if (incols == (msize_t)bottom[0].size())
652 resize(toprows + botrows, incols, (*tresult));
653 for (msize_t i = 0; i < botrows; ++i)
655 for (msize_t j = 0; j < incols; ++j)
657 (*tresult)[i + toprows][j] = bottom[i][j];
670 result = tempstorage;
675 void MatrixFunctions::getChunk(
const msize_t firstrow,
const msize_t lastrow,
const msize_t firstcol,
const msize_t lastcol,
const vector<vector<T> >& in, vector<vector<T> >& result)
677 msize_t outrows = lastrow - firstrow;
678 msize_t outcols = lastcol - firstcol;
679 if (lastrow <= firstrow || lastcol <= firstcol || firstrow < 0 || firstcol < 0 || lastrow > (msize_t)in.size() || lastcol > (msize_t)in[0].size())
684 vector<vector<T> > tempstorage, *tresult = &result;
685 bool copyout =
false;
689 tresult = &tempstorage;
691 resize(outrows, outcols, (*tresult),
true);
692 for (msize_t i = 0; i < outrows; ++i)
694 for (msize_t j = 0; j < outcols; ++j)
696 (*tresult)[i][j] = in[i + firstrow][j + firstcol];
701 result = tempstorage;