00001
00002
00003
00004 #include "include.h"
00005 #include "newmatap.h"
00006
00007 #include "tmt.h"
00008
00009 #ifdef use_namespace
00010 using namespace NEWMAT;
00011 #endif
00012
00013 ReturnMatrix Inverter(const CroutMatrix& X)
00014 {
00015 Matrix Y = X.i();
00016 Y.Release();
00017 return Y.ForReturn();
00018 }
00019
00020
00021 void trymatd()
00022 {
00023 Tracer et("Thirteenth test of Matrix package");
00024 Tracer::PrintTrace();
00025 Matrix X(5,20);
00026 int i,j;
00027 for (j=1;j<=20;j++) X(1,j) = j+1;
00028 for (i=2;i<=5;i++) for (j=1;j<=20; j++) X(i,j) = (long)X(i-1,j) * j % 1001;
00029 SymmetricMatrix S; S << X * X.t();
00030 Matrix SM = X * X.t() - S;
00031 Print(SM);
00032 LowerTriangularMatrix L = Cholesky(S);
00033 Matrix Diff = L*L.t()-S; Clean(Diff, 0.000000001);
00034 Print(Diff);
00035 {
00036 Tracer et1("Stage 1");
00037 LowerTriangularMatrix L1(5);
00038 Matrix Xt = X.t(); Matrix Xt2 = Xt;
00039 QRZT(X,L1);
00040 Diff = L - L1; Clean(Diff,0.000000001); Print(Diff);
00041 UpperTriangularMatrix Ut(5);
00042 QRZ(Xt,Ut);
00043 Diff = L - Ut.t(); Clean(Diff,0.000000001); Print(Diff);
00044 Matrix Y(3,20);
00045 for (j=1;j<=20;j++) Y(1,j) = 22-j;
00046 for (i=2;i<=3;i++) for (j=1;j<=20; j++)
00047 Y(i,j) = (long)Y(i-1,j) * j % 101;
00048 Matrix Yt = Y.t(); Matrix M,Mt; Matrix Y2=Y;
00049 QRZT(X,Y,M); QRZ(Xt,Yt,Mt);
00050 Diff = Xt - X.t(); Clean(Diff,0.000000001); Print(Diff);
00051 Diff = Yt - Y.t(); Clean(Diff,0.000000001); Print(Diff);
00052 Diff = Mt - M.t(); Clean(Diff,0.000000001); Print(Diff);
00053 Diff = Y2 * Xt2 * S.i() - M * L.i();
00054 Clean(Diff,0.000000001); Print(Diff);
00055 }
00056
00057 ColumnVector C1(5);
00058 {
00059 Tracer et1("Stage 2");
00060 X.ReSize(5,5);
00061 for (j=1;j<=5;j++) X(1,j) = j+1;
00062 for (i=2;i<=5;i++) for (j=1;j<=5; j++)
00063 X(i,j) = (long)X(i-1,j) * j % 1001;
00064 for (i=1;i<=5;i++) C1(i) = i*i;
00065 CroutMatrix A = X;
00066 ColumnVector C2 = A.i() * C1; C1 = X.i() * C1;
00067 X = C1 - C2; Clean(X,0.000000001); Print(X);
00068 }
00069
00070 {
00071 Tracer et1("Stage 3");
00072 X.ReSize(7,7);
00073 for (j=1;j<=7;j++) X(1,j) = j+1;
00074 for (i=2;i<=7;i++) for (j=1;j<=7; j++)
00075 X(i,j) = (long)X(i-1,j) * j % 1001;
00076 C1.ReSize(7);
00077 for (i=1;i<=7;i++) C1(i) = i*i;
00078 RowVector R1 = C1.t();
00079 Diff = R1 * X.i() - ( X.t().i() * R1.t() ).t(); Clean(Diff,0.000000001);
00080 Print(Diff);
00081 }
00082
00083 {
00084 Tracer et1("Stage 4");
00085 X.ReSize(5,5);
00086 for (j=1;j<=5;j++) X(1,j) = j+1;
00087 for (i=2;i<=5;i++) for (j=1;j<=5; j++)
00088 X(i,j) = (long)X(i-1,j) * j % 1001;
00089 C1.ReSize(5);
00090 for (i=1;i<=5;i++) C1(i) = i*i;
00091 CroutMatrix A1 = X*X;
00092 ColumnVector C2 = A1.i() * C1; C1 = X.i() * C1; C1 = X.i() * C1;
00093 X = C1 - C2; Clean(X,0.000000001); Print(X);
00094 }
00095
00096
00097 {
00098 Tracer et1("Stage 5");
00099 int n = 40;
00100 SymmetricBandMatrix B(n,2); B = 0.0;
00101 for (i=1; i<=n; i++)
00102 {
00103 B(i,i) = 6;
00104 if (i<=n-1) B(i,i+1) = -4;
00105 if (i<=n-2) B(i,i+2) = 1;
00106 }
00107 B(1,1) = 5; B(n,n) = 5;
00108 SymmetricMatrix A = B;
00109 ColumnVector X(n);
00110 X(1) = 429;
00111 for (i=2;i<=n;i++) X(i) = (long)X(i-1) * 31 % 1001;
00112 X = X / 100000L;
00113
00114
00115
00116 ColumnVector Y1 = A.i() * X;
00117 LowerTriangularMatrix C1 = Cholesky(A);
00118 ColumnVector Y2 = C1.t().i() * (C1.i() * X) - Y1;
00119 Clean(Y2, 0.000000001); Print(Y2);
00120 UpperTriangularMatrix CU = C1.t().i();
00121 LowerTriangularMatrix CL = C1.i();
00122 Y2 = CU * (CL * X) - Y1;
00123 Clean(Y2, 0.000000001); Print(Y2);
00124 Y2 = B.i() * X - Y1; Clean(Y2, 0.000000001); Print(Y2);
00125
00126 LowerBandMatrix C2 = Cholesky(B);
00127 Matrix M = C2 - C1; Clean(M, 0.000000001); Print(M);
00128 ColumnVector Y3 = C2.t().i() * (C2.i() * X) - Y1;
00129 Clean(Y3, 0.000000001); Print(Y3);
00130 CU = C1.t().i();
00131 CL = C1.i();
00132 Y3 = CU * (CL * X) - Y1;
00133 Clean(Y3, 0.000000001); Print(Y3);
00134
00135 Y3 = B.i() * X - Y1; Clean(Y3, 0.000000001); Print(Y3);
00136
00137 SymmetricMatrix AI = A.i();
00138 Y2 = AI*X - Y1; Clean(Y2, 0.000000001); Print(Y2);
00139 SymmetricMatrix BI = B.i();
00140 BandMatrix C = B; Matrix CI = C.i();
00141 M = A.i() - CI; Clean(M, 0.000000001); Print(M);
00142 M = B.i() - CI; Clean(M, 0.000000001); Print(M);
00143 M = AI-BI; Clean(M, 0.000000001); Print(M);
00144 M = AI-CI; Clean(M, 0.000000001); Print(M);
00145
00146 M = A; AI << M; M = AI-A; Clean(M, 0.000000001); Print(M);
00147 C = B; BI << C; M = BI-B; Clean(M, 0.000000001); Print(M);
00148
00149
00150 }
00151
00152 {
00153 Tracer et1("Stage 5");
00154 SymmetricMatrix A(4), B(4);
00155 A << 5
00156 << 1 << 4
00157 << 2 << 1 << 6
00158 << 1 << 0 << 1 << 7;
00159 B << 8
00160 << 1 << 5
00161 << 1 << 0 << 9
00162 << 2 << 1 << 0 << 6;
00163 LowerTriangularMatrix AB = Cholesky(A) * Cholesky(B);
00164 Matrix M = Cholesky(A) * B * Cholesky(A).t() - AB*AB.t();
00165 Clean(M, 0.000000001); Print(M);
00166 M = A * Cholesky(B); M = M * M.t() - A * B * A;
00167 Clean(M, 0.000000001); Print(M);
00168 }
00169 {
00170 Tracer et1("Stage 6");
00171 int N=49;
00172 int i;
00173 SymmetricBandMatrix S(N,1);
00174 Matrix B(N,N+1); B=0;
00175 for (i=1;i<=N;i++) { S(i,i)=1; B(i,i)=1; B(i,i+1)=-1; }
00176 for (i=1;i<N; i++) S(i,i+1)=-.5;
00177 DiagonalMatrix D(N+1); D = 1;
00178 B = B.t()*S.i()*B - (D-1.0/(N+1))*2.0;
00179 Clean(B, 0.000000001); Print(B);
00180 }
00181 {
00182 Tracer et1("Stage 7");
00183
00184 Matrix A(4,4);
00185 A.Row(1) << 3 << 2 << -1 << 4;
00186 A.Row(2) << -8 << 7 << 2 << 0;
00187 A.Row(3) << 2 << -2 << 3 << 1;
00188 A.Row(4) << -1 << 5 << 2 << 2;
00189 CroutMatrix B = A;
00190 Matrix C = A * Inverter(B) - IdentityMatrix(4);
00191 Clean(C, 0.000000001); Print(C);
00192 }
00193
00194
00195
00196 }