00001
00002
00003
00006
00007
00008 #define WANT_STREAM
00009
00010
00011
00012 #include "include.h"
00013
00014 #include "newmat.h"
00015
00016 #include "tmt.h"
00017
00018 #ifdef use_namespace
00019 using namespace NEWMAT;
00020 #endif
00021
00022
00023
00024
00025
00026 void trymat1()
00027 {
00028
00029 Tracer et("First test of Matrix package");
00030 Tracer::PrintTrace();
00031 {
00032 Tracer et1("Stage 1");
00033 int i,j;
00034
00035 LowerTriangularMatrix L(10);
00036 for (i=1;i<=10;i++) for (j=1;j<=i;j++) L(i,j)=2.0+i*i+j;
00037 SymmetricMatrix S(10);
00038 for (i=1;i<=10;i++) for (j=1;j<=i;j++) S(i,j)=i*j+1.0;
00039 SymmetricMatrix S1 = S / 2.0;
00040 S = S1 * 2.0;
00041 UpperTriangularMatrix U=L.t()*2.0;
00042 Print(LowerTriangularMatrix(L-U.t()*0.5));
00043 DiagonalMatrix D(10);
00044 for (i=1;i<=10;i++) D(i,i)=(i-4)*(i-5)*(i-6);
00045 Matrix M=(S+U-D+L)*(L+U-D+S);
00046 DiagonalMatrix DD=D*D;
00047 LowerTriangularMatrix LD=L*D;
00048
00049 Matrix M1 = S*L + U*L - D*L + L*L + 10.0;
00050 { M1 = M1 + S*U + U*U - D*U + L*U - S*D; }
00051 { M1 = M1 - U*D + DD - LD + S*S; }
00052 { M1 = M1 + U*S - D*S + L*S - 10.0; }
00053 M=M1-M;
00054 Print(M);
00055 }
00056 {
00057 Tracer et1("Stage 2");
00058 int i,j;
00059
00060 LowerTriangularMatrix L(9);
00061 for (i=1;i<=9;i++) for (j=1;j<=i;j++) L(i,j)=1.0+j;
00062 UpperTriangularMatrix U1(9);
00063 for (j=1;j<=9;j++) for (i=1;i<=j;i++) U1(i,j)=1.0+i;
00064 LowerTriangularMatrix LX(9);
00065 for (i=1;i<=9;i++) for (j=1;j<=i;j++) LX(i,j)=1.0+i*i;
00066 UpperTriangularMatrix UX(9);
00067 for (j=1;j<=9;j++) for (i=1;i<=j;i++) UX(i,j)=1.0+j*j;
00068 {
00069 L=L+LX/0.5; L=L-LX*3.0; L=LX*2.0+L;
00070 U1=U1+UX*2.0; U1=U1-UX*3.0; U1=UX*2.0+U1;
00071 }
00072
00073
00074 SymmetricMatrix S(9);
00075 for (i=1;i<=9;i++) for (j=1;j<=i;j++) S(i,j)=i*i+j;
00076 {
00077 SymmetricMatrix S1 = S;
00078 S=S1+5.0;
00079 S=S-3.0;
00080 }
00081
00082 DiagonalMatrix D(9);
00083 for (i=1;i<=9;i++) D(i,i)=S(i,i);
00084 UpperTriangularMatrix U=L.t()*2.0;
00085 {
00086 U1=U1*2.0 - U; Print(U1);
00087 L=L*2.0-D; U=U-D;
00088 }
00089 Matrix M=U+L; S=S*2.0; M=S-M; Print(M);
00090 }
00091 {
00092 Tracer et1("Stage 3");
00093 int i,j;
00094 Matrix M(10,3), N(10,3);
00095 for (i = 1; i<=10; i++) for (j = 1; j<=3; j++)
00096 { M(i,j) = 2*i-j; N(i,j) = i*j + 20; }
00097 Matrix MN = M + N, M1;
00098
00099 M1 = M; M1 += N; M1 -= MN; Print(M1);
00100 M1 = M; M1 += M1; M1 = M1 - M * 2; Print(M1);
00101 M1 = M; M1 += N * 2; M1 -= (MN + N); Print(M1);
00102 M1 = M; M1 -= M1; Print(M1);
00103 M1 = M; M1 -= MN + M1; M1 += N + M; Print(M1);
00104 M1 = M; M1 -= 5; M1 -= M; M1 *= 0.2; M1 = M1 + 1; Print(M1);
00105 Matrix NT = N.t();
00106 M1 = M; M1 *= NT; M1 -= M * N.t(); Print(M1);
00107 M = M * M.t();
00108 DiagonalMatrix D(10); D = 2;
00109 M1 = M; M1 += D; M1 -= M; M1 = M1 - D; Print(M1);
00110 M1 = M; M1 -= D; M1 -= M; M1 = M1 + D; Print(M1);
00111 M1 = M; M1 *= D; M1 /= 2; M1 -= M; Print(M1);
00112 SymmetricMatrix SM; SM << M;
00113
00114 SM += 10; M1 = SM - M; M1 /=10; M1 = M1 - 1; Print(M1);
00115 }
00116 {
00117 Tracer et1("Stage 4");
00118 int i,j;
00119 Matrix M(10,3), N(10,5);
00120 for (i = 1; i<=10; i++) for (j = 1; j<=3; j++) M(i,j) = 2*i-j;
00121 for (i = 1; i<=10; i++) for (j = 1; j<=5; j++) N(i,j) = i*j + 20;
00122 Matrix M1;
00123 M1 = M; M1 |= N; M1 &= N | M;
00124 M1 -= (M | N) & (N | M); Print(M1);
00125 M1 = M; M1 |= M1; M1 &= M1;
00126 M1 -= (M | M) & (M | M); Print(M1);
00127
00128 }
00129 {
00130 Tracer et1("Stage 5");
00131 int i,j;
00132 BandMatrix BM1(10,2,3), BM2(10,4,1); Matrix M1(10,10), M2(10,10);
00133 for (i=1;i<=10;i++) for (j=1;j<=10;j++)
00134 { M1(i,j) = 0.5*i+j*j-50; M2(i,j) = (i*101 + j*103) % 13; }
00135 BM1.Inject(M1); BM2.Inject(M2);
00136 BandMatrix BM = BM1; BM += BM2;
00137 Matrix M1X = BM1; Matrix M2X = BM2; Matrix MX = BM;
00138 MX -= M1X + M2X; Print(MX);
00139 MX = BM1; MX += BM2; MX -= M1X; MX -= M2X; Print(MX);
00140 SymmetricBandMatrix SM1; SM1 << BM1 * BM1.t();
00141 SymmetricBandMatrix SM2; SM2 << BM2 * BM2.t();
00142 SM1 *= 5.5;
00143 M1X *= M1X.t(); M1X *= 5.5; M2X *= M2X.t();
00144 SM1 -= SM2; M1 = SM1 - M1X + M2X; Print(M1);
00145 M1 = BM1; BM1 *= SM1; M1 = M1 * SM1 - BM1; Print(M1);
00146 M1 = BM1; BM1 -= SM1; M1 = M1 - SM1 - BM1; Print(M1);
00147 M1 = BM1; BM1 += SM1; M1 = M1 + SM1 - BM1; Print(M1);
00148
00149 }
00150 {
00151 Tracer et1("Stage 6");
00152 int i,j;
00153 Matrix M(10,10), N(10,10);
00154 for (i = 1; i<=10; i++) for (j = 1; j<=10; j++)
00155 { M(i,j) = 2*i-j; N(i,j) = i*j + 20; }
00156 GenericMatrix GM = M;
00157 GM += N; Matrix M1 = GM - N - M; Print(M1);
00158 DiagonalMatrix D(10); D = 3;
00159 GM = D; GM += N; GM += M; GM += D;
00160 M1 = D*2 - GM + M + N; Print(M1);
00161 GM = D; GM *= 4; GM += 16; GM /= 8; GM -= 2;
00162 GM -= D / 2; M1 = GM; Print(M1);
00163 GM = D; GM *= M; GM *= N; GM /= 3; M1 = M*N - GM; Print(M1);
00164 GM = D; GM |= M; GM &= N | D; M1 = GM - ((D | M) & (N | D));
00165 Print(M1);
00166 GM = M; M1 = M; GM += 5; GM *= 3; M *= 3; M += 15; M1 = GM - M;
00167 Print(M1);
00168 D.ReSize(10); for (i = 1; i<=10; i++) D(i) = i;
00169 M1 = D + 10; GM = D; GM += 10; M1 -= GM; Print(M1);
00170 GM = M; GM -= D; M1 = GM; GM = D; GM -= M; M1 += GM; Print(M1);
00171 GM = M; GM *= D; M1 = GM; GM = D; GM *= M.t();
00172 M1 -= GM.t(); Print(M1);
00173 GM = M; GM += 2 * GM; GM -= 3 * M; M1 = GM; Print(M1);
00174 GM = M; GM |= GM; GM -= (M | M); M1 = GM; Print(M1);
00175 GM = M; GM &= GM; GM -= (M & M); M1 = GM; Print(M1);
00176 M1 = M; M1 = (M1.t() & M.t()) - (M | M).t(); Print(M1);
00177 M1 = M; M1 = (M1.t() | M.t()) - (M & M).t(); Print(M1);
00178
00179 }
00180
00181 {
00182 Tracer et1("Stage 7");
00183
00184 int n = 3;
00185 int k; int j;
00186 Matrix A(n,n), B(n,n);
00187
00188
00189 for (k=1; k<=n; k++)
00190 {
00191 for (j = 1; j <= n; j++) A(k,j) = ((k-1) * (2*j-1));
00192 }
00193
00194
00195 for (k=1; k<=n; k++)
00196 {
00197 const int k1 = k-1;
00198 for (j = 1; j <= n; j++) B(k,j) = (k1 * (2*j-1));
00199 }
00200
00201 if (A != B)
00202 {
00203 cout << "\nVisual C++ version 5 compiler error?";
00204 cout << "\nTurn off optimisation";
00205 }
00206
00207 A -= B; Print(A);
00208
00209 }
00210
00211 {
00212 Tracer et1("Stage 8");
00213
00214 ColumnVector i(3); i << 1 << 0 << 0;
00215 ColumnVector j(3); j << 0 << 1 << 0;
00216 ColumnVector k(3); k << 0 << 0 << 1;
00217 ColumnVector X;
00218 X = CrossProduct(i,j) - k; Print(X);
00219 X = CrossProduct(j,k) - i; Print(X);
00220 X = CrossProduct(k,i) - j; Print(X);
00221 X = CrossProduct(j,i) + k; Print(X);
00222 X = CrossProduct(k,j) + i; Print(X);
00223 X = CrossProduct(i,k) + j; Print(X);
00224 X = CrossProduct(i,i); Print(X);
00225 X = CrossProduct(j,j); Print(X);
00226 X = CrossProduct(k,k); Print(X);
00227
00228 ColumnVector A(3); A << 2.25 << 1.75 << -7.5;
00229 ColumnVector B(3); B << -0.5 << 4.75 << 3.25;
00230 ColumnVector C(3); C << 9.25 << 3.5 << 1.25;
00231
00232 Real d0 = (A | B | C).Determinant();
00233 Real d1 = DotProduct(CrossProduct(A, B), C);
00234 Real d2 = DotProduct(CrossProduct(B, C), A);
00235 Real d3 = DotProduct(CrossProduct(C, A), B);
00236 X << (d1 - d0) << (d2 - d0) << (d3 - d0);
00237 Clean(X, 0.000000001); Print(X);
00238
00239 X = CrossProduct(A, CrossProduct(B, C))
00240 + CrossProduct(B, CrossProduct(C, A))
00241 + CrossProduct(C, CrossProduct(A, B));
00242 Print(X);
00243
00244 RowVector XT = CrossProduct(A.AsRow(), B.AsRow());
00245 XT -= CrossProduct(A, B).AsRow();
00246 Print(XT);
00247 }
00248
00249 {
00250 Tracer et1("Stage 9");
00251
00252 int i, j;
00253 Matrix M(10,3), N(10,3);
00254 for (i = 1; i<=10; i++) for (j = 1; j<=3; j++)
00255 { M(i,j) = 2*i-j; N(i,j) = i*j + 20; }
00256
00257 Matrix CP1 = CrossProductRows(M, N);
00258 Matrix CP2(10,3);
00259 for (i = 1; i<=10; i++)
00260 CP2.Row(i) = CrossProduct(M.Row(i), N.Row(i));
00261 CP2 -= CP1; Print(CP2);
00262
00263 CP2 = CrossProductColumns(M.t(), N.t());
00264 CP2 -= CP1.t(); Print(CP2);
00265 }
00266
00267 {
00268 Tracer et1("Stage 10");
00269
00270 MultWithCarry mwc;
00271 ColumnVector cv(10);
00272 for (int i = 1; i <= 10; ++i) cv(i) = mwc.Next();
00273 cv *= 100.0;
00274 cv(1) -= 6.27874;
00275 cv(2) -= 42.1718;
00276 cv(3) -= 80.2854;
00277 cv(4) -= 12.961;
00278 cv(5) -= 17.7499;
00279 cv(6) -= 13.2657;
00280 cv(7) -= 50.4923;
00281 cv(8) -= 26.095;
00282 cv(9) -= 57.9147;
00283 cv(10) -= 30.1778;
00284 Clean(cv, 0.0001); Print(cv);
00285 }
00286
00287
00288
00289 }
00290
00291