Здесь представлены шаблоны классов SLU2, SLU3 и SLU4, предназначенные для решения систем линейных уравнений второго, третьего и четвёртого порядка:
inline float qmod ( const float & x ) { return x * x; }
inline double qmod ( const double & x ) { return x * x; }
// Решение систем линейных уравнений 2-го порядка методом Гаусса
// Выбор ведущего элемента по столбцам
// aa * x + ab * y = ac
// ba * x + bb * y = bc
template <class T1, class T2 = T1> class SLU2 : public Matrix2<T1>
{
public:
T2 ac, bc;
// Решение системы линейных уравнений 2-го порядка методом Гаусса
bool gauss ( T2 & x, T2 & y ) const
{
T2 v1, v2;
const double ma = qmod ( aa );
const double mb = qmod ( ba );
if ( ma >= mb )
{
if ( ! ma ) return false;
const T1 c = ab / aa;
const T1 b = bb - c * ba;
if ( ! b ) return false;
( v1 = ac ) /= aa;
( v2 = v1 ) *= ba;
( ( y = bc ) -= v2 ) /= b;
( v2 = y ) *= c;
}
else
{
const T1 c = bb / ba;
const T1 a = ab - c * aa;
if ( a == 0 ) return false;
( v1 = bc ) /= ba;
( v2 = v1 ) *= aa;
( ( y = ac ) -= v2 ) /= a;
( v2 = y ) *= c;
}
( x = v1 ) -= v2;
return true;
}
SLU2 & fill ( const T1 & v1, const T2 & v2 )
{
set ( v1 );
ac = bc = v2;
return *this;
}
SLU2 & operator += ( const SLU2 & slu )
{
aa += slu.aa; ab += slu.ab; ac += slu.ac;
ba += slu.ba; bb += slu.bb; bc += slu.bc;
return *this;
}
};
// Решение систем линейных уравнений 3-го порядка методом Гаусса
// Выбор ведущего элемента по столбцам
// aa * x + ab * y + ac * z = ad
// ba * x + bb * y + bc * z = bd
// ca * x + cb * y + cc * z = cd
template <class T1, class T2 = T1> class SLU3 : public Matrix3<T1>
{
public:
T2 ad, bd, cd;
// Решение системы линейных уравнений 3-го порядка методом Гаусса
bool gauss ( T2 & x, T2 & y, T2 & z ) const
{
SLU2<T1, T2> slu;
T1 a, b;
T2 d, e;
const double ma = qmod ( ac );
const double mb = qmod ( bc );
const double mc = qmod ( cc );
if ( mc >= mb && mc >= ma )
{
if ( ! mc ) return false;
a = ca / cc;
b = cb / cc;
( d = cd ) /= cc;
slu.aa = aa - a * ac; slu.ab = ab - b * ac; ( e = d ) *= ac; ( slu.ac = ad ) -= e;
slu.ba = ba - a * bc; slu.bb = bb - b * bc; ( e = d ) *= bc; ( slu.bc = bd ) -= e;
}
else
if ( mb >= ma )
{
a = ba / bc;
b = bb / bc;
( d = bd ) /= bc;
slu.aa = aa - a * ac; slu.ab = ab - b * ac; ( e = d ) *= ac; ( slu.ac = ad ) -= e;
slu.ba = ca - a * cc; slu.bb = cb - b * cc; ( e = d ) *= cc; ( slu.bc = cd ) -= e;
}
else
{
a = aa / ac;
b = ab / ac;
( d = ad ) /= ac;
slu.aa = ba - a * bc; slu.ab = bb - b * bc; ( e = d ) *= bc; ( slu.ac = bd ) -= e;
slu.ba = ca - a * cc; slu.bb = cb - b * cc; ( e = d ) *= cc; ( slu.bc = cd ) -= e;
}
if ( ! slu.gauss ( x, y ) ) return false;
z = d;
( e = x ) *= a; z -= e;
( e = y ) *= b; z -= e;
return true;
}
SLU3 & fill ( const T1 & v1, const T2 & v2 )
{
set ( v1 );
ad = bd = cd = v2;
return *this;
}
};
// Решение симметричной системы линейных уравнений 3-го порядка методом Гаусса
// aa * x + ab * y + ac * z = ad
// ab * x + bb * y + bc * z = bd
// ac * x + bc * y + cc * z = cd
template <class T1, class T2=T1> class SymSLU3 : public SymMatrix3<T1>
{
public:
T2 ad, bd, cd;
SymSLU3 & operator += ( const SymSLU3 & slu )
{
aa += slu.aa; ab += slu.ab; ac += slu.ac; ad += slu.ad;
bb += slu.bb; bc += slu.bc; bd += slu.bd;
cc += slu.cc; cd += slu.cd;
return *this;
}
bool gauss ( T2 & x, T2 & y, T2 & z ) const
{
SLU2<T1, T2> slu;
T1 a, b;
T2 d, e;
const double ma = qmod ( ac );
const double mb = qmod ( bc );
const double mc = qmod ( cc );
if ( mc >= mb && mc >= ma )
{
if ( ! mc ) return false;
a = ac / cc;
b = bc / cc;
( d = cd ) /= cc;
slu.aa = aa - a * ac; slu.ab = ab - b * ac; ( e = d ) *= ac; ( slu.ac = ad ) -= e;
slu.ba = ab - a * bc; slu.bb = bb - b * bc; ( e = d ) *= bc; ( slu.bc = bd ) -= e;
}
else
if ( mb >= ma )
{
a = ab / bc;
b = bb / bc;
( d = bd ) /= bc;
slu.aa = aa - a * ac; slu.ab = ab - b * ac; ( e = d ) *= ac; ( slu.ac = ad ) -= e;
slu.ba = ac - a * cc; slu.bb = bc - b * cc; ( e = d ) *= cc; ( slu.bc = cd ) -= e;
}
else
{
a = aa / ac;
b = ab / ac;
( d = ad ) /= ac;
slu.aa = ab - a * bc; slu.ab = bb - b * bc; ( e = d ) *= bc; ( slu.ac = bd ) -= e;
slu.ba = ac - a * cc; slu.bb = bc - b * cc; ( e = d ) *= cc; ( slu.bc = cd ) -= e;
}
if ( ! slu.gauss ( x, y ) ) return false;
z = d;
( e = x ) *= a; z -= e;
( e = y ) *= b; z -= e;
return true;
}
SymSLU3 & fill ( const T1 & v1, const T2 & v2 )
{
set ( v1 );
ad = bd = cd = v2;
return *this;
}
};
// Решение систем линейных уравнений 4-го порядка методом Гаусса
// Выбор ведущего элемента по столбцам
// aa * xa + ab * xb + ac * xc + ad * xd = ae
// ba * xa + bb * xb + bc * xc + bd * xd = be
// ca * xa + cb * xb + cc * xc + cd * xd = ce
// da * xa + db * xb + dc * xc + dd * xd = de
template <class T1, class T2 = T1> class SLU4 : public Matrix4<T1>
{
public:
T2 ae, be, ce, de;
// Решение системы линейных уравнений 4-го порядка методом Гаусса
bool gauss ( T2 & xa, T2 & xb, T2 & xc, T2 & xd ) const
{
SLU3<T1, T2> slu;
T1 a, b, c;
T2 e, f;
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 false;
a = da / dd;
b = db / dd;
c = dc / dd;
e = de; e /= dd;
f = e; f *= ad;
slu.aa = aa - a * ad; slu.ab = ab - b * ad; slu.ac = ac - c * ad; slu.ad = ae; slu.ad -= f;
f = e; f *= bd;
slu.ba = ba - a * bd; slu.bb = bb - b * bd; slu.bc = bc - c * bd; slu.bd = be; slu.bd -= f;
f = e; f *= cd;
slu.ca = ca - a * cd; slu.cb = cb - b * cd; slu.cc = cc - c * cd; slu.cd = ce; slu.cd -= f;
}
else
if ( mc >= mb && mc >= ma )
{
a = ca / cd;
b = cb / cd;
c = cc / cd;
e = ce; e /= cd;
f = e; f *= ad;
slu.aa = aa - a * ad; slu.ab = ab - b * ad; slu.ac = ac - c * ad; slu.ad = ae; slu.ad -= f;
f = e; f *= bd;
slu.ba = ba - a * bd; slu.bb = bb - b * bd; slu.bc = bc - c * bd; slu.bd = be; slu.bd -= f;
f = e; f *= dd;
slu.ca = da - a * dd; slu.cb = db - b * dd; slu.cc = dc - c * dd; slu.cd = de; slu.cd -= f;
}
else
if ( mb >= ma )
{
a = ba / bd;
b = bb / bd;
c = bc / bd;
e = be; e /= bd;
f = e; f *= ad;
slu.aa = aa - a * ad; slu.ab = ab - b * ad; slu.ac = ac - c * ad; slu.ad = ae; slu.ad -= f;
f = e; f *= cd;
slu.ba = ca - a * cd; slu.bb = cb - b * cd; slu.bc = cc - c * cd; slu.bd = ce; slu.bd -= f;
f = e; f *= dd;
slu.ca = da - a * dd; slu.cb = db - b * dd; slu.cc = dc - c * dd; slu.cd = de; slu.cd -= f;
}
else
{
a = aa / ad;
b = ab / ad;
c = ac / ad;
e = ae; e /= ad;
f = e; f *= bd;
slu.aa = ba - a * bd; slu.ab = bb - b * bd; slu.ac = bc - c * bd; slu.ad = be; slu.ad -= f;
f = e; f *= cd;
slu.ba = ca - a * cd; slu.bb = cb - b * cd; slu.bc = cc - c * cd; slu.bd = ce; slu.bd -= f;
f = e; f *= dd;
slu.ca = da - a * dd; slu.cb = db - b * dd; slu.cc = dc - c * dd; slu.cd = de; slu.cd -= f;
}
if ( ! slu.gauss ( xa, xb, xc ) ) return false;
xd = e;
f = xa; f *= a; xd -= f;
f = xb; f *= b; xd -= f;
f = xc; f *= c; xd -= f;
return true;
}
SLU4 & fill ( const T1 & v1, const T2 & v2 )
{
set ( v1 );
ae = be = ce = de = v2;
return *this;
}
};
В качестве типа Т1 можно задавать действительные ( float, double ) или комплексные числа ( Complex ).
В качестве типа Т2, кроме одиночных чисел, можно задавать наборы чисел, такие как: Vector2d,
Vector3d, Vector4d, FixArray и т. д.
Исходники находятся в файле LinAlg.h.
|