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