Для всех матриц определены: заполнение матрицы указанным значением, получение обратной матрицы, вычисление детерминанта и оператор умножения. Матрицы 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;
}
Симметричные матрицы 2-го порядка:
template <class T> class SymMatrix2
{
public:
T aa, ab, bb;
SymMatrix2 & set ( const T & v )
{
aa = ab = bb = v;
return *this;
}
Def<SymMatrix2> operator ~ () const
{
SymSLU2<T, Set2
Матрицы 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;
}
Def<Matrix4> operator ~ () const
{
SLU4<T, Set4<T> > slu;
slu.base() = *this;
slu.ae.a = 1; slu.ae.b = 0; slu.ae.c = 0; slu.ae.d = 0;
slu.be.a = 0; slu.be.b = 1; slu.be.c = 0; slu.be.d = 0;
slu.ce.a = 0; slu.ce.b = 0; slu.ce.c = 1; slu.ce.d = 0;
slu.de.a = 0; slu.de.b = 0; slu.de.c = 0; slu.de.d = 1;
Def<Matrix4> m;
Set4<T> sa, sb, sc, sd;
if ( ! slu.gauss ( sa, sb, sc, sd ) ) return m;
m.isDef = true;
m.aa = sa.a; m.ab = sa.b; m.ac = sa.c; m.ad = sa.d;
m.ba = sb.a; m.bb = sb.b; m.bc = sb.c; m.bd = sb.d;
m.ca = sc.a; m.cb = sc.b; m.cc = sc.c; m.cd = sc.d;
m.da = sd.a; m.db = sd.b; m.dc = sd.c; m.dd = sd.d;
return m;
}
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();
}
}
};
template <class T> Matrix4<T> operator * ( const Matrix4<T> & l, const Matrix4<T> & r )
{
Matrix4<T> m;
m.aa = l.aa * r.aa + l.ab * r.ba + l.ac * r.ca + l.ad * r.da; m.ab = l.aa * r.ab + l.ab * r.bb + l.ac * r.cb + l.ad * r.db;
m.ac = l.aa * r.ac + l.ab * r.bc + l.ac * r.cc + l.ad * r.dc; m.ad = l.aa * r.ad + l.ab * r.bd + l.ac * r.cd + l.ad * r.dd;
m.ba = l.ba * r.aa + l.bb * r.ba + l.bc * r.ca + l.bd * r.da; m.bb = l.ba * r.ab + l.bb * r.bb + l.bc * r.cb + l.bd * r.db;
m.bc = l.ba * r.ac + l.bb * r.bc + l.bc * r.cc + l.bd * r.dc; m.bd = l.ba * r.ad + l.bb * r.bd + l.bc * r.cd + l.bd * r.dd;
m.ca = l.ca * r.aa + l.cb * r.ba + l.cc * r.ca + l.cd * r.da; m.cb = l.ca * r.ab + l.cb * r.bb + l.cc * r.cb + l.cd * r.db;
m.cc = l.ca * r.ac + l.cb * r.bc + l.cc * r.cc + l.cd * r.dc; m.cd = l.ca * r.ad + l.cb * r.bd + l.cc * r.cd + l.cd * r.dd;
m.da = l.da * r.aa + l.db * r.ba + l.dc * r.ca + l.dd * r.da; m.db = l.da * r.ab + l.db * r.bb + l.dc * r.cb + l.dd * r.db;
m.dc = l.da * r.ac + l.db * r.bc + l.dc * r.cc + l.dd * r.dc; m.dd = l.da * r.ad + l.db * r.bd + l.dc * r.cd + l.dd * r.dd;
return m;
}
Исходники находятся в файле LinAlg.h.
|