📄 tmth.cpp
字号:
/// \ingroup newmat
///@{
/// \file tmth.cpp
/// Part of matrix library test program.
//#define WANT_STREAM
#include "include.h"
#include "newmatap.h"
//#include "newmatio.h"
#include "tmt.h"
#ifdef use_namespace
using namespace NEWMAT;
#endif
static int my_max(int p, int q) { return (p > q) ? p : q; }
static int my_min(int p, int q) { return (p < q) ? p : q; }
void BandFunctions(int l1, int u1, int l2, int u2)
{
int i, j;
BandMatrix BM1(20, l1, u1); BM1 = 999999.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
if (i - j <= l1 && i - j >= -u1) BM1(i, j) = 100 * i + j;
BandMatrix BM2 = BM1; Matrix M = BM2 - BM1; Print(M);
BM2.ReSize(20, l2, u2); BM2 = 777777.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
if (i - j <= l2 && i - j >= -u2) BM2(i, j) = (100 * i + j) * 11;
BandMatrix BMA = BM1 + BM2, BMS = BM1 - BM2, BMSP = SP(BM1, BM2),
BMM = BM1 * BM2, BMN = -BM1;
Matrix M1(20,20); M1 = 0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
if (i - j <= l1 && i - j >= -u1) M1(i, j) = 100 * i + j;
Matrix M2(20,20); M2 = 0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
if (i - j <= l2 && i - j >= -u2) M2(i, j) = (100 * i + j) * 11;
Matrix MA = M1 + M2, MS = M1 - M2, MSP = SP(M1, M2), MM = M1 * M2, MN = -M1;
MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
Print(MA); Print(MS); Print(MSP); Print(MM); Print(MN);
Matrix Test(7, 2);
Test(1,1) = BM1.BandWidth().Lower() - l1;
Test(1,2) = BM1.BandWidth().Upper() - u1;
Test(2,1) = BM2.BandWidth().Lower() - l2;
Test(2,2) = BM2.BandWidth().Upper() - u2;
Test(3,1) = BMA.BandWidth().Lower() - my_max(l1,l2);
Test(3,2) = BMA.BandWidth().Upper() - my_max(u1,u2);
Test(4,1) = BMS.BandWidth().Lower() - my_max(l1,l2);
Test(4,2) = BMS.BandWidth().Upper() - my_max(u1,u2);
Test(5,1) = BMSP.BandWidth().Lower() - my_min(l1,l2);
Test(5,2) = BMSP.BandWidth().Upper() - my_min(u1,u2);
Test(6,1) = BMM.BandWidth().Lower() - (l1 + l2);
Test(6,2) = BMM.BandWidth().Upper() - (u1 + u2);
Test(7,1) = BMN.BandWidth().Lower() - l1;
Test(7,2) = BMN.BandWidth().Upper() - u1;
Print(Test);
// test element function
BandMatrix BM1E(20, l1, u1); BM1E = 999999.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
if (i - j <= l1 && i - j >= -u1) BM1E.element(i-1, j-1) = 100 * i + j;
BandMatrix BM2E = BM1E; BM2E.ReSize(BM2); BM2E = 777777.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
if (i - j <= l2 && i - j >= -u2)
BM2E.element(i-1, j-1) = (100 * i + j) * 11;
M1 = BM1E - BM1; Print(M1);
M2 = BM2E - BM2; Print(M2);
// test element function with constant
BM1E = 444444.04; BM2E = 555555.0;
const BandMatrix BM1C = BM1, BM2C = BM2;
for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
if (i - j <= l1 && i - j >= -u1)
BM1E.element(i-1, j-1) = BM1C.element(i-1, j-1);
for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
if (i - j <= l2 && i - j >= -u2)
BM2E.element(i-1, j-1) = BM2C.element(i-1, j-1);
M1 = BM1E - BM1; Print(M1);
M2 = BM2E - BM2; Print(M2);
// test subscript function with constant
BM1E = 444444.04; BM2E = 555555.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
if (i - j <= l1 && i - j >= -u1) BM1E(i, j) = BM1C(i, j);
for (i = 1; i <= 20; ++i) for (j = 1; j <= 20; ++j)
if (i - j <= l2 && i - j >= -u2) BM2E(i, j) = BM2C(i, j);
M1 = BM1E - BM1; Print(M1);
M2 = BM2E - BM2; Print(M2);
}
void LowerBandFunctions(int l1, int l2)
{
int i, j;
LowerBandMatrix BM1(20, l1); BM1 = 999999.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l1) BM1(i, j) = 100 * i + j;
LowerBandMatrix BM2 = BM1; Matrix M = BM2 - BM1; Print(M);
BM2.ReSize(20, l2); BM2 = 777777.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l2) BM2(i, j) = (100 * i + j) * 11;
LowerBandMatrix BMA = BM1 + BM2, BMS = BM1 - BM2, BMSP = SP(BM1, BM2),
BMM = BM1 * BM2, BMN = -BM1;
Matrix M1(20,20); M1 = 0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l1) M1(i, j) = 100 * i + j;
Matrix M2(20,20); M2 = 0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l2) M2(i, j) = (100 * i + j) * 11;
Matrix MA = M1 + M2, MS = M1 - M2, MSP = SP(M1, M2), MM = M1 * M2, MN = -M1;
MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
Print(MA); Print(MS); Print(MSP); Print(MM); Print(MN);
Matrix Test(7, 2);
Test(1,1) = BM1.BandWidth().Lower() - l1;
Test(1,2) = BM1.BandWidth().Upper();
Test(2,1) = BM2.BandWidth().Lower() - l2;
Test(2,2) = BM2.BandWidth().Upper();
Test(3,1) = BMA.BandWidth().Lower() - my_max(l1,l2);
Test(3,2) = BMA.BandWidth().Upper();
Test(4,1) = BMS.BandWidth().Lower() - my_max(l1,l2);
Test(4,2) = BMS.BandWidth().Upper();
Test(5,1) = BMSP.BandWidth().Lower() - my_min(l1,l2);
Test(5,2) = BMSP.BandWidth().Upper();
Test(6,1) = BMM.BandWidth().Lower() - (l1 + l2);
Test(6,2) = BMM.BandWidth().Upper();
Test(7,1) = BMN.BandWidth().Lower() - l1;
Test(7,2) = BMN.BandWidth().Upper();
Print(Test);
// test element function
LowerBandMatrix BM1E(20, l1); BM1E = 999999.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l1) BM1E.element(i-1, j-1) = 100 * i + j;
LowerBandMatrix BM2E = BM1E; BM2E.ReSize(BM2); BM2E = 777777.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l2) BM2E.element(i-1, j-1) = (100 * i + j) * 11;
M1 = BM1E - BM1; Print(M1);
M2 = BM2E - BM2; Print(M2);
// test element function with constant
BM1E = 444444.04; BM2E = 555555.0;
const LowerBandMatrix BM1C = BM1, BM2C = BM2;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l1) BM1E.element(i-1, j-1) = BM1C.element(i-1, j-1);
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l2) BM2E.element(i-1, j-1) = BM2C.element(i-1, j-1);
M1 = BM1E - BM1; Print(M1);
M2 = BM2E - BM2; Print(M2);
// test subscript function with constant
BM1E = 444444.04; BM2E = 555555.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l1) BM1E(i, j) = BM1C(i, j);
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l2) BM2E(i, j) = BM2C(i, j);
M1 = BM1E - BM1; Print(M1);
M2 = BM2E - BM2; Print(M2);
}
void UpperBandFunctions(int u1, int u2)
{
int i, j;
UpperBandMatrix BM1(20, u1); BM1 = 999999.0;
for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
if (i - j >= -u1) BM1(i, j) = 100 * i + j;
UpperBandMatrix BM2 = BM1; Matrix M = BM2 - BM1; Print(M);
BM2.ReSize(20, u2); BM2 = 777777.0;
for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
if (i - j >= -u2) BM2(i, j) = (100 * i + j) * 11;
UpperBandMatrix BMA = BM1 + BM2, BMS = BM1 - BM2, BMSP = SP(BM1, BM2),
BMM = BM1 * BM2, BMN = -BM1;
Matrix M1(20,20); M1 = 0;
for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
if (i - j >= -u1) M1(i, j) = 100 * i + j;
Matrix M2(20,20); M2 = 0;
for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
if (i - j >= -u2) M2(i, j) = (100 * i + j) * 11;
Matrix MA = M1 + M2, MS = M1 - M2, MSP = SP(M1, M2), MM = M1 * M2, MN = -M1;
MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
Print(MA); Print(MS); Print(MSP); Print(MM); Print(MN);
Matrix Test(7, 2);
Test(1,1) = BM1.BandWidth().Lower();
Test(1,2) = BM1.BandWidth().Upper() - u1;
Test(2,1) = BM2.BandWidth().Lower();
Test(2,2) = BM2.BandWidth().Upper() - u2;
Test(3,1) = BMA.BandWidth().Lower();
Test(3,2) = BMA.BandWidth().Upper() - my_max(u1,u2);
Test(4,1) = BMS.BandWidth().Lower();
Test(4,2) = BMS.BandWidth().Upper() - my_max(u1,u2);
Test(5,1) = BMSP.BandWidth().Lower();
Test(5,2) = BMSP.BandWidth().Upper() - my_min(u1,u2);
Test(6,1) = BMM.BandWidth().Lower();
Test(6,2) = BMM.BandWidth().Upper() - (u1 + u2);
Test(7,1) = BMN.BandWidth().Lower();
Test(7,2) = BMN.BandWidth().Upper() - u1;
Print(Test);
// test element function
UpperBandMatrix BM1E(20, u1); BM1E = 999999.0;
for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
if (i - j >= -u1) BM1E.element(i-1, j-1) = 100 * i + j;
UpperBandMatrix BM2E = BM1E; BM2E.ReSize(BM2); BM2E = 777777.0;
for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
if (i - j >= -u2) BM2E.element(i-1, j-1) = (100 * i + j) * 11;
M1 = BM1E - BM1; Print(M1);
M2 = BM2E - BM2; Print(M2);
// test element function with constant
BM1E = 444444.04; BM2E = 555555.0;
const UpperBandMatrix BM1C = BM1, BM2C = BM2;
for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
if (i - j >= -u1) BM1E.element(i-1, j-1) = BM1C.element(i-1, j-1);
for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
if (i - j >= -u2) BM2E.element(i-1, j-1) = BM2C.element(i-1, j-1);
M1 = BM1E - BM1; Print(M1);
M2 = BM2E - BM2; Print(M2);
// test subscript function with constant
BM1E = 444444.04; BM2E = 555555.0;
for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
if (i - j >= -u1) BM1E(i, j) = BM1C(i, j);
for (i = 1; i <= 20; ++i) for (j = i; j <= 20; ++j)
if (i - j >= -u2) BM2E(i, j) = BM2C(i, j);
M1 = BM1E - BM1; Print(M1);
M2 = BM2E - BM2; Print(M2);
}
void SymmetricBandFunctions(int l1, int l2)
{
int i, j;
SymmetricBandMatrix BM1(20, l1); BM1 = 999999.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l1) BM1(i, j) = 100 * i + j;
SymmetricBandMatrix BM2 = BM1; Matrix M = BM2 - BM1; Print(M);
BM2.ReSize(20, l2); BM2 = 777777.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l2) BM2(i, j) = (100 * i + j) * 11;
SymmetricBandMatrix BMA = BM1 + BM2, BMS = BM1 - BM2, BMSP = SP(BM1, BM2),
BMN = -BM1;
BandMatrix BMM = BM1 * BM2;
SymmetricMatrix M1(20); M1 = 0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l1) M1(i, j) = 100 * i + j;
SymmetricMatrix M2(20); M2 = 0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l2) M2(i, j) = (100 * i + j) * 11;
SymmetricMatrix MA = M1 + M2, MS = M1 - M2, MSP = SP(M1, M2), MN = -M1;
Matrix MM = M1 * M2;
MA -= BMA; MS -= BMS; MSP -= BMSP; MM -= BMM; MN -= BMN;
Print(MA); Print(MS); Print(MSP); Print(MM); Print(MN);
Matrix Test(7, 2);
Test(1,1) = BM1.BandWidth().Lower() - l1;
Test(1,2) = BM1.BandWidth().Upper() - l1;
Test(2,1) = BM2.BandWidth().Lower() - l2;
Test(2,2) = BM2.BandWidth().Upper() - l2;
Test(3,1) = BMA.BandWidth().Lower() - my_max(l1,l2);
Test(3,2) = BMA.BandWidth().Upper() - my_max(l1,l2);
Test(4,1) = BMS.BandWidth().Lower() - my_max(l1,l2);
Test(4,2) = BMS.BandWidth().Upper() - my_max(l1,l2);
Test(5,1) = BMSP.BandWidth().Lower() - my_min(l1,l2);
Test(5,2) = BMSP.BandWidth().Upper() - my_min(l1,l2);
Test(6,1) = BMM.BandWidth().Lower() - (l1 + l2);
Test(6,2) = BMM.BandWidth().Upper() - (l1 + l2);
Test(7,1) = BMN.BandWidth().Lower() - l1;
Test(7,2) = BMN.BandWidth().Upper() - l1;
Print(Test);
// test element function
SymmetricBandMatrix BM1E(20, l1); BM1E = 999999.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l1) BM1E.element(i-1, j-1) = 100 * i + j;
SymmetricBandMatrix BM2E = BM1E; BM2E.ReSize(BM2); BM2E = 777777.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l2) BM2E.element(i-1, j-1) = (100 * i + j) * 11;
M1 = BM1E - BM1; Print(M1);
M2 = BM2E - BM2; Print(M2);
// reverse subscripts
BM1E = 999999.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l1) BM1E.element(j-1, i-1) = 100 * i + j;
BM2E = 777777.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l2) BM2E.element(j-1, i-1) = (100 * i + j) * 11;
M1 = BM1E - BM1; Print(M1);
M2 = BM2E - BM2; Print(M2);
// test element function with constant
BM1E = 444444.04; BM2E = 555555.0;
const SymmetricBandMatrix BM1C = BM1, BM2C = BM2;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l1) BM1E.element(i-1, j-1) = BM1C.element(i-1, j-1);
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l2) BM2E.element(i-1, j-1) = BM2C.element(i-1, j-1);
M1 = BM1E - BM1; Print(M1);
M2 = BM2E - BM2; Print(M2);
// reverse subscripts
BM1E = 444444.0; BM2E = 555555.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l1) BM1E.element(j-1, i-1) = BM1C.element(j-1, i-1);
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l2) BM2E.element(j-1, i-1) = BM2C.element(j-1, i-1);
M1 = BM1E - BM1; Print(M1);
M2 = BM2E - BM2; Print(M2);
// test subscript function with constant
BM1E = 444444.0; BM2E = 555555.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l1) BM1E(i, j) = BM1C(i, j);
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l2) BM2E(i, j) = BM2C(i, j);
M1 = BM1E - BM1; Print(M1);
M2 = BM2E - BM2; Print(M2);
// reverse subscripts
BM1E = 444444.0; BM2E = 555555.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l1) BM1E(j, i) = BM1C(j, i);
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l2) BM2E(j, i) = BM2C(j, i);
M1 = BM1E - BM1; Print(M1);
M2 = BM2E - BM2; Print(M2);
// partly reverse subscripts
BM1E = 444444.0; BM2E = 555555.0;
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l1) BM1E(j, i) = BM1C(i, j);
for (i = 1; i <= 20; ++i) for (j = 1; j <= i; ++j)
if (i - j <= l2) BM2E(j, i) = BM2C(i, j);
M1 = BM1E - BM1; Print(M1);
M2 = BM2E - BM2; Print(M2);
}
void trymath()
{
// cout << "\nSeventeenth test of Matrix package\n";
// cout << "\n";
Tracer et("Seventeenth test of Matrix package");
Tracer::PrintTrace();
{
Tracer et1("Stage 1");
int i, j;
BandMatrix B(8,3,1);
for (i=1; i<=8; i++) for (j=-3; j<=1; j++)
{ int k = i+j; if (k>0 && k<=8) B(i,k) = i + k/64.0; }
IdentityMatrix I(8); BandMatrix B1; B1 = I;
UpperTriangularMatrix UT = I; Print(Matrix(B1-UT));
Print(Matrix(B * B - B * 2 + I - (B - I) * (B - I)));
Matrix A = B; BandMatrix C; C = B;
Print(Matrix(B * A - C * 2 + I - (A - I) * (B - I)));
C.ReSize(8,2,2); C = 0.0; C.Inject(B);
Matrix X = A.t();
B1.ReSize(8,2,2); B1.Inject(X); Print(Matrix(C.t()-B1));
Matrix BI = B.i(); A = A.i()-BI; Clean(A,0.000000001); Print(A);
BandLUMatrix BLU = B.t();
BI = BLU.i(); A = X.i()-BI; Clean(A,0.000000001); Print(A);
X.ReSize(1,1);
X(1,1) = BLU.LogDeterminant().Value()-B.LogDeterminant().Value();
Clean(X,0.000000001); Print(X);
UpperBandMatrix U; U << B; LowerBandMatrix L; L << B;
DiagonalMatrix D; D << B;
Print( Matrix(L + (U - D - B)) );
for (i=1; i<=8; i++) A.Column(i) << B.Column(i);
Print(Matrix(A-B));
}
{
Tracer et1("Stage 2");
BandMatrix A(7,2,2);
int i,j;
for (i=1; i<=7; i++) for (j=1; j<=7; j++)
{
int k=i-j; if (k<0) k = -k;
if (k==0) A(i,j)=6;
else if (k==1) A(i,j) = -4;
else if (k==2) A(i,j) = 1;
A(1,1) = A(7,7) = 5;
}
DiagonalMatrix D(7); D = 0.0; A = A - D;
BandLUMatrix B(A); Matrix M = A;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -