tmt7.cpp
Go to the documentation of this file.
00001 
00002 
00003 
00006 
00007 
00008 //#define WANT_STREAM
00009 
00010 #include "include.h"
00011 
00012 #include "newmat.h"
00013 
00014 #include "tmt.h"
00015 
00016 #ifdef use_namespace
00017 using namespace NEWMAT;
00018 #endif
00019 
00020 // C matrix mulitply - for testing RealStarStar
00021 
00022 void c_matrix_multiply(int p, int q, int r,
00023    const Real** a, const Real** b, Real** c)
00024 {
00025    for (int i = 0; i < p; ++i) for (int k = 0; k < r; ++k)
00026    {
00027       Real sum = 0.0;
00028       for (int j = 0; j < q; ++j) sum += a[i][j] * b[j][k];
00029       c[i][k] = sum;
00030    }
00031 }
00032    
00033 
00034 
00035 void trymat7()
00036 {
00037 //   cout << "\nSeventh test of Matrix package\n";
00038    Tracer et("Seventh test of Matrix package");
00039    Tracer::PrintTrace();
00040 
00041    int i,j;
00042 
00043 
00044    DiagonalMatrix D(6);
00045    UpperTriangularMatrix U(6);
00046    for (i=1;i<=6;i++) { for (j=i;j<=6;j++) U(i,j)=i*i*j-50; D(i,i)=i*i+i-10; }
00047    LowerTriangularMatrix L=(U*3.0).t();
00048    SymmetricMatrix S(6);
00049    for (i=1;i<=6;i++) for (j=i;j<=6;j++) S(i,j)=i*i+2.0+j;
00050    Matrix MD=D; Matrix ML=L; Matrix MU=U;
00051    Matrix MS=S;
00052    Matrix M(6,6);
00053    for (i=1;i<=6;i++) for (j=1;j<=6;j++) M(i,j)=i*j+i*i-10.0;  
00054    {
00055       Tracer et1("Stage 1");
00056       Print(Matrix((S-M)-(MS-M)));
00057       Print(Matrix((-M-S)+(MS+M)));
00058       Print(Matrix((U-M)-(MU-M)));
00059    }
00060    {
00061       Tracer et1("Stage 2");
00062       Print(Matrix((L-M)+(M-ML)));
00063       Print(Matrix((D-M)+(M-MD)));
00064       Print(Matrix((D-S)+(MS-MD)));
00065       Print(Matrix((D-L)+(ML-MD)));
00066    }
00067 
00068    { M=MU.t(); }
00069    LowerTriangularMatrix LY=D.i()*U.t();
00070    {
00071       Tracer et1("Stage 3");
00072       MS=D*LY-M; Clean(MS,0.00000001); Print(MS);
00073       L=U.t();
00074       LY=D.i()*L; MS=D*LY-M; Clean(MS,0.00000001); Print(MS);
00075    }
00076    {
00077       Tracer et1("Stage 4");
00078       UpperTriangularMatrix UT(11);
00079       int i, j;
00080       for (i=1;i<=11;i++) for (j=i;j<=11;j++) UT(i,j)=i*i+j*3;
00081       GenericMatrix GM; Matrix X;
00082       UpperBandMatrix UB(11,3); UB.Inject(UT); UT = UB;
00083       UpperBandMatrix UB2 = UB / 8;
00084       GM = UB2-UT/8; X = GM; Print(X);
00085       SymmetricBandMatrix SB(11,4); SB << (UB + UB.t());
00086       X = SB - UT - UT.t(); Print(X);
00087       BandMatrix B = UB + UB.t()*2;
00088       DiagonalMatrix D; D << B;
00089       X.ReSize(1,1); X(1,1) = Trace(B)-Sum(D); Print(X);
00090       X = SB + 5; Matrix X1=X; X = SP(UB,X); Matrix X2 =UB;
00091       X1 = (X1.AsDiagonal() * X2.AsDiagonal()).AsRow()-X.AsColumn().t();
00092       Print(X1);
00093       X1=SB.t(); X2 = B.t(); X = SB.i() * B - X1.i() * X2.t();
00094       Clean(X,0.00000001); Print(X);
00095       X = SB.i(); X = X * B - X1.i() * X2.t();
00096       Clean(X,0.00000001); Print(X);
00097       D = 1; X = SB.i() * SB - D; Clean(X,0.00000001); Print(X);
00098       ColumnVector CV(11);
00099       CV << 2 << 6 <<3 << 8 << -4 << 17.5 << 2 << 1 << -2 << 5 << 3.75;
00100       D << 2 << 6 <<3 << 8 << -4 << 17.5 << 2 << 1 << -2 << 5 << 3.75;
00101       X = CV.AsDiagonal(); X = X-D; Print(X);
00102       SymmetricBandMatrix SB1(11,7); SB1 = 5; 
00103       SymmetricBandMatrix SB2 = SB1 + D;
00104       X.ReSize(11,11); X=0;
00105       for (i=1;i<=11;i++) for (j=1;j<=11;j++)
00106       {
00107          if (abs(i-j)<=7) X(i,j)=5;
00108          if (i==j) X(i,j)+=CV(i);
00109       }
00110       SymmetricMatrix SM; SM.ReSize(11);
00111       SM=SB; SB = SB+SB2; X1 = SM+X-SB; Print(X1);
00112       SB2=0; X2=SB2; X1=SB; Print(X2);
00113       for (i=1;i<=11;i++) SB2.Column(i)<<SB.Column(i);
00114       X1=X1-SB2; Print(X1);
00115       X = SB; SB2.ReSize(11,4); SB2 = SB*5; SB2 = SB + SB2;
00116       X1 = X*6 - SB2; Print(X1);
00117       X1 = SP(SB,SB2/3); X1=X1-SP(X,X*2); Print(X1);
00118       X1 = SP(SB2/6,X*2); X1=X1-SP(X*2,X); Print(X1);
00119    }
00120 
00121    {
00122       // test the simple integer array class
00123       Tracer et("Stage 5");
00124       ColumnVector Test(10); Test = 0.0;
00125       int i;
00126       SimpleIntArray A(100);
00127       for (i = 0; i < 100; i++) A[i] = i*i+1;
00128       SimpleIntArray B(100), C(50), D;
00129       B = A; A.ReSize(50, true); C = A; A.ReSize(150, true); D = A;
00130       for (i = 0; i < 100; i++) if (B[i] != i*i+1) Test(1)=1;
00131       for (i = 0; i < 50; i++) if (C[i] != i*i+1) Test(2)=1;
00132       for (i = 0; i < 50; i++) if (D[i] != i*i+1) Test(3)=1;
00133       for (i = 50; i < 150; i++) if (D[i] != 0) Test(3)=1;
00134       A.resize(75); A = A.size();
00135       for (i = 0; i < 75; i++) if (A[i] != 75) Test(4)=1;
00136       A.resize(25); A = A.size();
00137       for (i = 0; i < 25; i++) if (A[i] != 25) Test(5)=1;
00138       A.ReSize(25); A = 23;
00139       for (i = 0; i < 25; i++) if (A[i] != 23) Test(6)=1;
00140       A.ReSize(0); A.ReSize(15); A = A.Size();
00141       for (i = 0; i < 15; i++) if (A[i] != 15) Test(7)=1;
00142       const SimpleIntArray E = B;
00143       for (i = 0; i < 100; i++) if (E[i] != i*i+1) Test(8)=1;
00144       SimpleIntArray F; F.resize_keep(5);
00145       for (i = 0; i < 5; i++) if (F[i] != 0) Test(9)=1;
00146       Print(Test);
00147    }
00148    
00149    {
00150       // testing RealStarStar
00151       Tracer et("Stage 6");
00152       MultWithCarry MWC;
00153       
00154       Matrix A(10, 12), B(12, 15), C(10, 15);
00155       FillWithValues(MWC, A); FillWithValues(MWC, B);
00156       ConstRealStarStar a(A);
00157       ConstRealStarStar b(B);
00158       RealStarStar c(C);
00159       c_matrix_multiply(10,12,15,a,b,c);
00160       Matrix X = C - A * B; Clean(X,0.00000001); Print(X);
00161       A.ReSize(11, 10); B.ReSize(10,8); C.ReSize(11,8);
00162       FillWithValues(MWC, A); FillWithValues(MWC, B);
00163       C = -1;
00164       c_matrix_multiply(11,10,8,
00165          ConstRealStarStar(A),ConstRealStarStar(B),RealStarStar(C));
00166       X = C - A * B; Clean(X,0.00000001); Print(X);
00167    }
00168    
00169    {
00170       // testing resize_keep
00171       Tracer et("Stage 7");
00172       Matrix X, Y;
00173       MultWithCarry MWC;
00174       
00175       X.resize(20,35); FillWithValues(MWC, X);
00176       Matrix M(20,35); M = X;
00177       X = M.submatrix(1,15,1,25);
00178       M.resize_keep(15,25); Y  = X - M; Print(Y);
00179       M.resize_keep(15,25); Y  = X - M; Print(Y);
00180       Y.resize(29,27); Y = 0; Y.submatrix(1,15,1,25) = X;
00181       M.resize_keep(29,27); Y -= M; Print(Y);
00182       M.resize_keep(0,5);  M.resize_keep(10,10); Print(M);
00183       M.resize_keep(15,0); M.resize_keep(10,10); Print(M);
00184       
00185       X.resize(20,35); FillWithValues(MWC, X);
00186       M = X;
00187       M.resize_keep(38,17);
00188       Y.resize(38,17); Y = 0;
00189       Y.submatrix(1,20,1,17) = X.submatrix(1,20,1,17);
00190       Y -= M; Print(Y);
00191             
00192       X.resize(40,12); FillWithValues(MWC, X);
00193       M = X;
00194       M.resize_keep(38,17);
00195       Y.resize(38,17); Y = 0;
00196       Y.submatrix(1,38,1,12) = X.submatrix(1,38,1,12);
00197       Y -= M; Print(Y);
00198 
00199 #ifndef DONT_DO_NRIC
00200 
00201       X.resize(20,35); FillWithValues(MWC, X);
00202       nricMatrix nM(20,35); nM = X;
00203       X = nM.submatrix(1,15,1,25);
00204       nM.resize_keep(15,25); Y  = X - nM; Print(Y);
00205       nM.resize_keep(15,25); Y  = X - nM; Print(Y);
00206       Y.resize(29,27); Y = 0; Y.submatrix(1,15,1,25) = X;
00207       nM.resize_keep(29,27); Y -= nM; Print(Y);
00208       nM.resize_keep(0,5);  nM.resize_keep(10,10); Print(nM);
00209       nM.resize_keep(15,0); nM.resize_keep(10,10); Print(nM);
00210       
00211       X.resize(20,35); FillWithValues(MWC, X);
00212       nM = X;
00213       nM.resize_keep(38,17);
00214       Y.resize(38,17); Y = 0;
00215       Y.submatrix(1,20,1,17) = X.submatrix(1,20,1,17);
00216       Y -= nM; Print(Y);
00217             
00218       X.resize(40,12); FillWithValues(MWC, X);
00219       nM = X;
00220       nM.resize_keep(38,17);
00221       Y.resize(38,17); Y = 0;
00222       Y.submatrix(1,38,1,12) = X.submatrix(1,38,1,12);
00223       Y -= nM; Print(Y);      
00224 
00225 #endif      
00226       
00227       X.resize(20,20); FillWithValues(MWC, X);
00228       SquareMatrix SQM(20); SQM << X;
00229       X = SQM.sym_submatrix(1,13);
00230       SQM.resize_keep(13); Y  = X - SQM; Print(Y);
00231       SQM.resize_keep(13); Y  = X - SQM; Print(Y);
00232       Y.resize(23,23); Y = 0; Y.sym_submatrix(1,13) = X;
00233       SQM.resize_keep(23,23); Y -= SQM; Print(Y);
00234       SQM.resize_keep(0); SQM.resize_keep(50); Print(SQM);
00235       
00236       X.resize(20,20); FillWithValues(MWC, X);
00237       SymmetricMatrix SM(20); SM << X;
00238       X = SM.sym_submatrix(1,13);
00239       SM.resize_keep(13); Y  = X - SM; Print(Y);
00240       SM.resize_keep(13); Y  = X - SM; Print(Y);
00241       Y.resize(23,23); Y = 0; Y.sym_submatrix(1,13) = X;
00242       SM.resize_keep(23); Y -= SM; Print(Y);
00243       SM.resize_keep(0); SM.resize_keep(50); Print(SM);
00244       
00245       X.resize(20,20); FillWithValues(MWC, X);
00246       LowerTriangularMatrix LT(20); LT << X;
00247       X = LT.sym_submatrix(1,13);
00248       LT.resize_keep(13); Y  = X - LT; Print(Y);
00249       LT.resize_keep(13); Y  = X - LT; Print(Y);
00250       Y.resize(23,23); Y = 0; Y.sym_submatrix(1,13) = X;
00251       LT.resize_keep(23); Y -= LT; Print(Y);
00252       LT.resize_keep(0); LT.resize_keep(50); Print(LT);
00253       
00254       X.resize(20,20); FillWithValues(MWC, X);
00255       UpperTriangularMatrix UT(20); UT << X;
00256       X = UT.sym_submatrix(1,13);
00257       UT.resize_keep(13); Y  = X - UT; Print(Y);
00258       UT.resize_keep(13); Y  = X - UT; Print(Y);
00259       Y.resize(23,23); Y = 0; Y.sym_submatrix(1,13) = X;
00260       UT.resize_keep(23); Y -= UT; Print(Y);
00261       UT.resize_keep(0); UT.resize_keep(50); Print(UT);
00262       
00263       X.resize(20,20); FillWithValues(MWC, X);
00264       DiagonalMatrix DM(20); DM << X;
00265       X = DM.sym_submatrix(1,13);
00266       DM.resize_keep(13); Y  = X - DM; Print(Y);
00267       DM.resize_keep(13); Y  = X - DM; Print(Y);
00268       Y.resize(23,23); Y = 0; Y.sym_submatrix(1,13) = X;
00269       DM.resize_keep(23); Y -= DM; Print(Y);
00270       DM.resize_keep(0); DM.resize_keep(50); Print(DM);
00271       
00272       X.resize(1,20); FillWithValues(MWC, X);
00273       RowVector RV(20); RV << X;
00274       X = RV.columns(1,13);
00275       RV.resize_keep(13); Y  = X - RV; Print(Y);
00276       RV.resize_keep(13); Y  = X - RV; Print(Y);
00277       Y.resize(1,23); Y = 0; Y.columns(1,13) = X;
00278       RV.resize_keep(1,23); Y -= RV; Print(Y);
00279       RV.resize_keep(0); RV.resize_keep(50); Print(RV);
00280       
00281       X.resize(20,1); FillWithValues(MWC, X);
00282       ColumnVector CV(20); CV << X;
00283       X = CV.rows(1,13);
00284       CV.resize_keep(13); Y  = X - CV; Print(Y);
00285       CV.resize_keep(13); Y  = X - CV; Print(Y);
00286       Y.resize(23,1); Y = 0; Y.rows(1,13) = X;
00287       CV.resize_keep(23,1); Y -= CV; Print(Y);
00288       CV.resize_keep(0); CV.resize_keep(50); Print(CV);
00289       
00290       
00291    }    
00292 
00293 
00294 //   cout << "\nEnd of seventh test\n";
00295 }
00296 
00297 
00298 


kni
Author(s): Martin Günther
autogenerated on Mon Aug 14 2017 02:44:13