00001
00002 #define WANT_STREAM
00003
00004 #include "include.h"
00005
00006 #include "newmatap.h"
00007 #include "newmatio.h"
00008
00009 #include "tmt.h"
00010
00011 #ifdef use_namespace
00012 using namespace NEWMAT;
00013 #endif
00014
00015
00016 static inline int my_min(int x, int y) { return x < y ? x : y; }
00017 static inline int my_max(int x, int y) { return x > y ? x : y; }
00018
00019
00020 #ifdef SETUP_C_SUBSCRIPTS
00021
00022 void trymatk()
00023 {
00024 Tracer et("Twentieth test of Matrix package");
00025 Tracer::PrintTrace();
00026
00027
00028
00029 int i,j; Matrix X, Y;
00030
00031 cout << "Matrix\n";
00032 Matrix A(15,35), B(15, 35);
00033 for (i=0; i<15; i++) for (j=0; j<35; j++)
00034 { A[i][j] = i+100*j; B(i+1,j+1) = i+100*j; }
00035 X = A - B; Print(X); Y = X;
00036 for (i=0; i<15; i++) for (j=0; j<35; j++)
00037 {
00038 X.element(i,j) = A.element(i,j) - B[i][j];
00039 Y.element(i,j) = ((const Matrix&)A)[i][j] - B[i][j];
00040 }
00041 Print(X); Print(Y);
00042 A.CleanUp(); B.CleanUp();
00043
00044 cout << "UpperTriangularMatrix\n";
00045 UpperTriangularMatrix A1(15), B1(15);
00046 for (i=0; i<15; i++) for (j=i; j<15; j++)
00047 { A1[i][j] = i+100*j; B1(i+1,j+1) = i+100*j; }
00048 X = A1 - B1; Print(X); Y = X;
00049 for (i=0; i<15; i++) for (j=i; j<15; j++)
00050 {
00051 X.element(i,j) = A1.element(i,j) - B1[i][j];
00052 Y.element(i,j) = ((const UpperTriangularMatrix&)A1)[i][j] - B1[i][j];
00053 }
00054 Print(X); Print(Y);
00055 A1.CleanUp(); B1.CleanUp();
00056
00057 cout << "LowerTriangularMatrix\n";
00058 LowerTriangularMatrix A2(35), B2(35);
00059 for (i=0; i<35; i++) for (j=0; j<=i; j++)
00060 { A2[i][j] = i+100*j; B2(i+1,j+1) = i+100*j; }
00061 X = A2 - B2; Print(X); Y = X;
00062 for (i=0; i<35; i++) for (j=0; j<=i; j++)
00063 {
00064 X.element(i,j) = A2.element(i,j) - B2[i][j];
00065 Y.element(i,j) = ((const LowerTriangularMatrix&)A2)[i][j] - B2[i][j];
00066 }
00067 Print(X); Print(Y);
00068 A2.CleanUp(); B2.CleanUp();
00069
00070 cout << "SymmetricMatrix\n";
00071 SymmetricMatrix A3(10), B3(10);
00072 for (i=0; i<10; i++) for (j=0; j<=i; j++)
00073 { A3[i][j] = i+100*j; B3(i+1,j+1) = i+100*j; }
00074 X = A3 - B3; Print(X); Y = X;
00075 for (i=0; i<10; i++) for (j=0; j<=i; j++)
00076 {
00077 X.element(i,j) = A3.element(i,j) - B3[i][j];
00078 Y.element(i,j) = ((const SymmetricMatrix&)A3)[i][j] - B3[i][j];
00079 }
00080 Print(X); Print(Y);
00081 A3.CleanUp(); B3.CleanUp();
00082
00083 cout << "DiagonalMatrix\n";
00084 DiagonalMatrix A4(10), B4(10);
00085 for (i=0; i<10; i++)
00086 { A4[i] = i+100; B4(i+1) = i+100; }
00087 X = A4 - B4; Print(X); Y = X;
00088 for (i=0; i<10; i++)
00089 {
00090 X.element(i,i) = A4.element(i) - B4[i];
00091 Y.element(i,i) = ((const DiagonalMatrix&)A4)[i] - B4[i];
00092 }
00093 Print(X); Print(Y);
00094 A4.CleanUp(); B4.CleanUp();
00095
00096 cout << "RowVector\n";
00097 RowVector A5(10), B5(10);
00098 for (i=0; i<10; i++)
00099 { A5[i] = i+100; B5(i+1) = i+100; }
00100 X = A5 - B5; Print(X); Y = X;
00101 for (i=0; i<10; i++)
00102 {
00103 X.element(0,i) = A5.element(i) - B5[i];
00104 Y.element(0,i) = ((const RowVector&)A5)[i] - B5[i];
00105 }
00106 Print(X); Print(Y);
00107 A5.CleanUp(); B5.CleanUp();
00108
00109 cout << "ColumnVector\n";
00110 ColumnVector A6(10), B6(10);
00111 for (i=0; i<10; i++)
00112 { A6[i] = i+100; B6(i+1) = i+100; }
00113 X = A6 - B6; Print(X); Y = X;
00114 for (i=0; i<10; i++)
00115 {
00116 X.element(i,0) = A6.element(i) - B6[i];
00117 Y.element(i,0) = ((const ColumnVector&)A6)[i] - B6[i];
00118 }
00119 Print(X); Print(Y);
00120 A6.CleanUp(); B6.CleanUp();
00121
00122 cout << "BandMatrix\n";
00123 BandMatrix A7(55,10, 5), B7(55, 10, 5);
00124 for (i=0; i<55; i++) for (j=my_max(0,i-10); j<=my_min(54,i+5); j++)
00125 { A7[i][j] = i+100*j; B7(i+1,j+1) = i+100*j; }
00126 X = A7 - B7; Print(X); Y = X;
00127 for (i=0; i<55; i++) for (j=my_max(0,i-10); j<=my_min(54,i+5); j++)
00128 {
00129 X.element(i,j) = A7.element(i,j) - B7[i][j];
00130 Y.element(i,j) = ((const BandMatrix&)A7)[i][j] - B7[i][j];
00131 }
00132 Print(X); Print(Y);
00133 A7.CleanUp(); B7.CleanUp();
00134
00135 cout << "UpperBandMatrix\n";
00136 UpperBandMatrix A8(80,15), B8(80,15);
00137 for (i=0; i<80; i++) for (j=i; j<=my_min(79,i+15); j++)
00138 { A8[i][j] = i+100*j; B8(i+1,j+1) = i+100*j; }
00139 X = A8 - B8; Print(X); Y = X;
00140 for (i=0; i<80; i++) for (j=i; j<=my_min(79,i+15); j++)
00141 {
00142 X.element(i,j) = A8.element(i,j) - B8[i][j];
00143 Y.element(i,j) = ((const UpperBandMatrix&)A8)[i][j] - B8[i][j];
00144 }
00145 Print(X); Print(Y);
00146 A8.CleanUp(); B8.CleanUp();
00147
00148 cout << "LowerBandMatrix\n";
00149 LowerBandMatrix A9(75,27), B9(75,27);
00150 for (i=0; i<75; i++) for (j=my_max(0,i-27); j<=i; j++)
00151 { A9[i][j] = i+100*j; B9(i+1,j+1) = i+100*j; }
00152 X = A9 - B9; Print(X); Y = X;
00153 for (i=0; i<75; i++) for (j=my_max(0,i-27); j<=i; j++)
00154 {
00155 X.element(i,j) = A9.element(i,j) - B9[i][j];
00156 Y.element(i,j) = ((const LowerBandMatrix&)A9)[i][j] - B9[i][j];
00157 }
00158 Print(X); Print(Y);
00159 A9.CleanUp(); B9.CleanUp();
00160
00161 cout << "SymmetricBandMatrix\n";
00162 SymmetricBandMatrix Aa(69,15), Ba(69,15);
00163 for (i=0; i<69; i++) for (j=my_max(0,i-15); j<=i; j++)
00164 { Aa[i][j] = i+100*j; Ba(i+1,j+1) = i+100*j; }
00165 X = Aa - Ba; Print(X); Y = X;
00166 for (i=0; i<69; i++) for (j=my_max(0,i-15); j<=i; j++)
00167 {
00168 X.element(i,j) = Aa.element(i,j) - Ba[i][j];
00169 Y.element(i,j) = ((const SymmetricBandMatrix&)Aa)[i][j] - Ba[i][j];
00170 }
00171 Print(X); Print(Y);
00172 Aa.CleanUp(); Ba.CleanUp();
00173
00174
00175 }
00176
00177
00178 #else
00179
00180 void trymatk()
00181 {
00182 Tracer et("Twentieth test of Matrix package");
00183 Tracer::PrintTrace();
00184
00185 cout << "C subscripts not enabled, not tested\n\n";
00186 }
00187
00188 #endif