Для квадратной матрицы 2-го порядка определены: заполнение матрицы указанным значением, получение обратной матрицы, вычисление детерминанта и оператор умножения. template <class T> class Matrix2 { public: T aa, ab, ba, bb; Matrix2 & set ( const T & v ) { aa = ab = ba = bb = v; return *this; } Def<Matrix2> operator ~ () const { SLU2<T, Set2<T>> slu; slu.base() = *this; slu.ac.a = 1; slu.bc.a = 0; slu.ac.b = 0; slu.bc.b = 1; Def<Matrix2> m; Set2<T> sa, sb; if ( ! slu.gauss ( sa, sb ) ) return m; m.isDef = true; m.aa = sa.a; m.ab = sa.b; m.ba = sb.a; m.bb = sb.b; return m; } T determinant() const { return aa * bb - ab * ba; } }; template <class T> Matrix2<T> operator * ( const Matrix2<T> & l, const Matrix2<T> & r ) { Matrix2<T> m; m.aa = l.aa * r.aa + l.ab * r.ba; m.ab = l.aa * r.ab + l.ab * r.bb; m.ba = l.ba * r.aa + l.bb * r.ba; m.bb = l.ba * r.ab + l.bb * r.bb; return m; }Для квадратной матрицы 3-го порядка определены: заполнение матрицы указанным значением, получение обратной матрицы, вычисление детерминанта и оператор умножения. template <class T> class Matrix3 { public: T aa, ab, ac, ba, bb, bc, ca, cb, cc; Matrix3 & set ( const T & v ) { aa = ab = ac = ba = bb = bc = ca = cb = cc = v; return *this; } Def<Matrix3> operator ~ () const { SLU3<T, Set3<T>> slu; slu.base() = *this; slu.ad.a = 1; slu.ad.b = 0; slu.ad.c = 0; slu.bd.a = 0; slu.bd.b = 1; slu.bd.c = 0; slu.cd.a = 0; slu.cd.b = 0; slu.cd.c = 1; Def<Matrix3> m; Set3<T> sa, sb, sc; if ( ! slu.gauss ( sa, sb, sc ) ) return m; m.isDef = true; m.aa = sa.a; m.ab = sa.b; m.ac = sa.c; m.ba = sb.a; m.bb = sb.b; m.bc = sb.c; m.ca = sc.a; m.cb = sc.b; m.cc = sc.c; return m; } T determinant() const { return ( ab * bc - ac * bb ) * ca + ( ac * ba - aa * bc ) * cb + ( aa * bb - ab * ba ) * cc; } }; template <class T> Matrix3<T> operator * ( const Matrix3<T> & l, const Matrix3<T> & r ) { Matrix3<T> m; m.aa = l.aa * r.aa + l.ab * r.ba + l.ac * r.ca; m.ab = l.aa * r.ab + l.ab * r.bb + l.ac * r.cb; m.ac = l.aa * r.ac + l.ab * r.bc + l.ac * r.cc; m.ba = l.ba * r.aa + l.bb * r.ba + l.bc * r.ca; m.bb = l.ba * r.ab + l.bb * r.bb + l.bc * r.cb; m.bc = l.ba * r.ac + l.bb * r.bc + l.bc * r.cc; m.ca = l.ca * r.aa + l.cb * r.ba + l.cc * r.ca; m.cb = l.ca * r.ab + l.cb * r.bb + l.cc * r.cb; m.cc = l.ca * r.ac + l.cb * r.bc + l.cc * r.cc; return m; }Для симметричной квадратной матрицы 3-го порядка определены: заполнение матрицы указанным значением, получение обратной матрицы, вычисление детерминанта и оператор умножения. template <class T> class SymMatrix3 { public: T aa, ab, ac, bb, bc, cc; SymMatrix3 & set ( const T & v ) { aa = ab = ac = bb = bc = cc = v; return *this; } Def<SymMatrix3> operator ~ () const { SymSLU3<T, Set3<T>> slu; slu.base() = *this; slu.ad.a = 1; slu.ad.b = 0; slu.ad.c = 0; slu.bd.a = 0; slu.bd.b = 1; slu.bd.c = 0; slu.cd.a = 0; slu.cd.b = 0; slu.cd.c = 1; Def<SymMatrix3> m; Set3<T> sa, sb, sc; if ( ! slu.gauss ( sa, sb, sc ) ) return m; m.isDef = true; m.aa = sa.a; m.ab = sa.b; m.ac = sa.c; m.bb = sb.b; m.bc = sb.c; m.cc = sc.c; return m; } T determinant() const { return ( ab * bc - ac * bb ) * ac + ( ac * ab - aa * bc ) * bc + ( aa * bb - ab * ab ) * cc; } }; template <class T> Matrix3<T> operator * ( const SymMatrix3<T> & l, const SymMatrix3<T> & r ) { Matrix3<T> m; m.aa = l.aa * r.aa + l.ab * r.ab + l.ac * r.ac; m.ab = l.aa * r.ab + l.ab * r.bb + l.ac * r.bc; m.ac = l.aa * r.ac + l.ab * r.bc + l.ac * r.cc; m.ba = l.ab * r.aa + l.bb * r.ab + l.bc * r.ac; m.bb = l.ab * r.ab + l.bb * r.bb + l.bc * r.bc; m.bc = l.ab * r.ac + l.bb * r.bc + l.bc * r.cc; m.ca = l.ac * r.aa + l.bc * r.ab + l.cc * r.ac; m.cb = l.ac * r.ab + l.bc * r.bb + l.cc * r.bc; m.cc = l.ac * r.ac + l.bc * r.bc + l.cc * r.cc; return m; }Для квадратной матрицы 4-го порядка определены: заполнение матрицы указанным значением и вычисление детерминанта. template <class T> class Matrix4 { public: T aa, ab, ac, ad, ba, bb, bc, bd, ca, cb, cc, cd, da, db, dc, dd; Matrix4 & set ( const T & v ) { aa = ab = ac = ad = ba = bb = bc = bd = ca = cb = cc = cd = da = db = dc = dd = v; return *this; } T determinant() const { Matrix3<T> m; T a, b, c; const double ma = qmod ( ad ); const double mb = qmod ( bd ); const double mc = qmod ( cd ); const double md = qmod ( dd ); if ( md >= mc && md >= mb && md >= ma ) { if ( dd == 0 ) return 0; a = da / dd; b = db / dd; c = dc / dd; m.aa = aa - a * ad; m.ab = ab - b * ad; m.ac = ac - c * ad; m.ba = ba - a * bd; m.bb = bb - b * bd; m.bc = bc - c * bd; m.ca = ca - a * cd; m.cb = cb - b * cd; m.cc = cc - c * cd; return dd * m.determinant(); } if ( mc >= mb && mc >= ma ) { a = ca / cd; b = cb / cd; c = cc / cd; m.aa = aa - a * ad; m.ab = ab - b * ad; m.ac = ac - c * ad; m.ba = ba - a * bd; m.bb = bb - b * bd; m.bc = bc - c * bd; m.ca = da - a * dd; m.cb = db - b * dd; m.cc = dc - c * dd; return -cd * m.determinant(); } if ( mb >= ma ) { a = ba / bd; b = bb / bd; c = bc / bd; m.aa = aa - a * ad; m.ab = ab - b * ad; m.ac = ac - c * ad; m.ba = ca - a * cd; m.bb = cb - b * cd; m.bc = cc - c * cd; m.ca = da - a * dd; m.cb = db - b * dd; m.cc = dc - c * dd; return bd * m.determinant(); } else { a = aa / ad; b = ab / ad; c = ac / ad; m.aa = ba - a * bd; m.ab = bb - b * bd; m.ac = bc - c * bd; m.ba = ca - a * cd; m.bb = cb - b * cd; m.bc = cc - c * cd; m.ca = da - a * dd; m.cb = db - b * dd; m.cc = dc - c * dd; return -ad * m.determinant(); } } }; Исходники находятся в файле LinAlg.h. Наверх |