home *** CD-ROM | disk | FTP | other *** search
Text File | 1992-12-06 | 59.0 KB | 1,992 lines |
- Newsgroups: comp.sources.misc
- From: robertd@kauri.vuw.ac.nz (Robert Davies)
- Subject: v34i013: newmat06 - A matrix package in C++, Part07/07
- Message-ID: <1992Dec6.045901.4846@sparky.imd.sterling.com>
- X-Md4-Signature: e4653de0d2693bd454a2450abcae8eef
- Date: Sun, 6 Dec 1992 04:59:01 GMT
- Approved: kent@sparky.imd.sterling.com
-
- Submitted-by: robertd@kauri.vuw.ac.nz (Robert Davies)
- Posting-number: Volume 34, Issue 13
- Archive-name: newmat06/part07
- Environment: C++
- Supersedes: newmat04: Volume 26, Issue 87-91
-
- #! /bin/sh
- # This is a shell archive. Remove anything before this line, then unpack
- # it by saving it into a file and typing "sh file". To overwrite existing
- # files, type "sh file -c". You can also feed this as standard input via
- # unshar, or by typing "sh <file", e.g.. If this archive is complete, you
- # will see the following message at the end:
- # "End of archive 7 (of 7)."
- # Contents: tmt1.cxx tmt2.cxx tmt3.cxx tmt4.cxx tmt5.cxx tmt6.cxx
- # tmt7.cxx tmt8.cxx tmt9.cxx tmta.cxx tmtb.cxx tmtc.cxx tmtd.cxx
- # tmte.cxx tmtf.cxx tmtg.cxx tmth.cxx tmti.cxx
- # Wrapped by robert@kea on Thu Dec 3 22:37:51 1992
- PATH=/bin:/usr/bin:/usr/ucb ; export PATH
- if test -f 'tmt1.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmt1.cxx'\"
- else
- echo shar: Extracting \"'tmt1.cxx'\" \(2304 characters\)
- sed "s/^X//" >'tmt1.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X
- X
- X
- X#include "include.h"
- X
- X#include "newmat.h"
- X
- X/**************************** test program ******************************/
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- X
- Xvoid trymat1()
- X{
- X// cout << "\nFirst test of Matrix package\n\n";
- X Tracer et("First test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X {
- X Tracer et1("Stage 1");
- X int i,j;
- X
- X LowerTriangularMatrix L(10);
- X for (i=1;i<=10;i++) for (j=1;j<=i;j++) L(i,j)=2.0+i*i+j;
- X SymmetricMatrix S(10);
- X for (i=1;i<=10;i++) for (j=1;j<=i;j++) S(i,j)=i*j+1.0;
- X SymmetricMatrix S1 = S / 2.0;
- X S = S1 * 2.0;
- X UpperTriangularMatrix U=L.t()*2.0;
- X Print(LowerTriangularMatrix(L-U.t()*0.5));
- X DiagonalMatrix D(10);
- X for (i=1;i<=10;i++) D(i,i)=(i-4)*(i-5)*(i-6);
- X Matrix M=(S+U-D+L)*(L+U-D+S);
- X DiagonalMatrix DD=D*D;
- X LowerTriangularMatrix LD=L*D;
- X // expressions split for Turbo C
- X Matrix M1 = S*L + U*L - D*L + L*L + 10.0;
- X { M1 = M1 + S*U + U*U - D*U + L*U - S*D; }
- X { M1 = M1 - U*D + DD - LD + S*S; }
- X { M1 = M1 + U*S - D*S + L*S - 10.0; }
- X M=M1-M;
- X Print(M);
- X }
- X {
- X Tracer et1("Stage 2");
- X int i,j;
- X
- X LowerTriangularMatrix L(9);
- X for (i=1;i<=9;i++) for (j=1;j<=i;j++) L(i,j)=1.0+j;
- X UpperTriangularMatrix U1(9);
- X for (j=1;j<=9;j++) for (i=1;i<=j;i++) U1(i,j)=1.0+i;
- X LowerTriangularMatrix LX(9);
- X for (i=1;i<=9;i++) for (j=1;j<=i;j++) LX(i,j)=1.0+i*i;
- X UpperTriangularMatrix UX(9);
- X for (j=1;j<=9;j++) for (i=1;i<=j;i++) UX(i,j)=1.0+j*j;
- X {
- X L=L+LX/0.5; L=L-LX*3.0; L=LX*2.0+L;
- X U1=U1+UX*2.0; U1=U1-UX*3.0; U1=UX*2.0+U1;
- X }
- X
- X
- X SymmetricMatrix S(9);
- X for (i=1;i<=9;i++) for (j=1;j<=i;j++) S(i,j)=i*i+j;
- X {
- X SymmetricMatrix S1 = S;
- X S=S1+5.0;
- X S=S-3.0;
- X }
- X
- X DiagonalMatrix D(9);
- X for (i=1;i<=9;i++) D(i,i)=S(i,i);
- X UpperTriangularMatrix U=L.t()*2.0;
- X {
- X U1=U1*2.0 - U; Print(U1);
- X L=L*2.0-D; U=U-D;
- X }
- X Matrix M=U+L; S=S*2.0; M=S-M; Print(M);
- X }
- X// cout << "\nEnd of first test\n";
- X}
- X
- END_OF_FILE
- if test 2304 -ne `wc -c <'tmt1.cxx'`; then
- echo shar: \"'tmt1.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmt1.cxx'
- fi
- if test -f 'tmt2.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmt2.cxx'\"
- else
- echo shar: Extracting \"'tmt2.cxx'\" \(1761 characters\)
- sed "s/^X//" >'tmt2.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X
- X
- X#include "include.h"
- X
- X#include "newmat.h"
- X
- X/**************************** test program ******************************/
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- X
- Xvoid trymat2()
- X{
- X// cout << "\nSecond test of Matrix package\n\n";
- X Tracer et("Second test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X
- X int i,j;
- X Matrix M(3,5);
- X for (i=1; i<=3; i++) for (j=1; j<=5; j++) M(i,j) = 100*i + j;
- X Matrix X(8,10);
- X for (i=1; i<=8; i++) for (j=1; j<=10; j++) X(i,j) = 1000*i + 10*j;
- X Matrix Y = X; Matrix Z = X;
- X { X.SubMatrix(2,4,3,7) << M; }
- X for (i=1; i<=3; i++) for (j=1; j<=5; j++) Y(i+1,j+2) = 100*i + j;
- X Print(Matrix(X-Y));
- X
- X
- X Real a[15]; Real* r = a;
- X for (i=1; i<=3; i++) for (j=1; j<=5; j++) *r++ = 100*i + j;
- X { Z.SubMatrix(2,4,3,7) << a; }
- X Print(Matrix(Z-Y));
- X
- X { M=33; X.SubMatrix(2,4,3,7) << M; }
- X { Z.SubMatrix(2,4,3,7) << 33; }
- X Print(Matrix(Z-X));
- X
- X for (i=1; i<=8; i++) for (j=1; j<=10; j++) X(i,j) = 1000*i + 10*j;
- X Y = X;
- X UpperTriangularMatrix U(5);
- X for (i=1; i<=5; i++) for (j=i; j<=5; j++) U(i,j) = 100*i + j;
- X { X.SubMatrix(3,7,5,9) << U; }
- X for (i=1; i<=5; i++) for (j=i; j<=5; j++) Y(i+2,j+4) = 100*i + j;
- X for (i=1; i<=5; i++) for (j=1; j<i; j++) Y(i+2,j+4) = 0.0;
- X Print(Matrix(X-Y));
- X for (i=1; i<=8; i++) for (j=1; j<=10; j++) X(i,j) = 1000*i + 10*j;
- X Y = X;
- X for (i=1; i<=5; i++) for (j=i; j<=5; j++) U(i,j) = 100*i + j;
- X { X.SubMatrix(3,7,5,9).Inject(U); }
- X for (i=1; i<=5; i++) for (j=i; j<=5; j++) Y(i+2,j+4) = 100*i + j;
- X Print(Matrix(X-Y));
- X
- X
- X// cout << "\nEnd of second test\n";
- X}
- END_OF_FILE
- if test 1761 -ne `wc -c <'tmt2.cxx'`; then
- echo shar: \"'tmt2.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmt2.cxx'
- fi
- if test -f 'tmt3.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmt3.cxx'\"
- else
- echo shar: Extracting \"'tmt3.cxx'\" \(2192 characters\)
- sed "s/^X//" >'tmt3.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X
- X#include "include.h"
- X
- X#include "newmat.h"
- X
- X
- X/**************************** test program ******************************/
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- Xvoid trymat3()
- X{
- X// cout << "\nThird test of Matrix package\n";
- X Tracer et("Third test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X
- X int i,j;
- X
- X {
- X Tracer et1("Stage 1");
- X SymmetricMatrix S(7);
- X for (i=1;i<=7;i++) for (j=1;j<=i;j++) S(i,j)=i*i+j;
- X S=-S+2.0;
- X
- X DiagonalMatrix D(7);
- X for (i=1;i<=7;i++) D(i,i)=S(i,i);
- X
- X Matrix M4(7,7); { M4=D+(D+4.0); M4=M4-D*2.0; M4=M4-4.0; Print(M4); }
- X SymmetricMatrix S2=D; Matrix M2=S2; { M2=-D+M2; Print(M2); }
- X UpperTriangularMatrix U2=D; { M2=U2; M2=D-M2; Print(M2); }
- X LowerTriangularMatrix L2=D; { M2=L2; M2=D-M2; Print(M2); }
- X M2=D; M2=M2-D; Print(M2);
- X for (i=1;i<=7;i++) for (j=1;j<=i;j++) L2(i,j)=2.0-i*i-j;
- X U2=L2.t(); D=D.t(); S=S.t();
- X M4=(L2-1.0)+(U2+1.0)-D-S; Print(M4);
- X M4=(-L2+1.0)+(-U2-1.0)+D+S; Print(M4);
- X }
- X
- X {
- X Tracer et1("Stage 2");
- X DiagonalMatrix D(6);
- X for (i=1;i<=6;i++) D(i,i)=i*3.0+i*i+2.0;
- X UpperTriangularMatrix U2(7); LowerTriangularMatrix L2(7);
- X for (i=1;i<=7;i++) for (j=1;j<=i;j++) L2(i,j)=2.0-i*i+j;
- X { U2=L2.t(); }
- X DiagonalMatrix D1(7); for (i=1;i<=7;i++) D1(i,i)=(i-2)*(i-4);
- X Matrix M2(6,7);
- X for (i=1;i<=6;i++) for (j=1;j<=7;j++) M2(i,j)=2.0+i*j+i*i+2.0*j*j;
- X Matrix MD=D; SymmetricMatrix MD1(1); MD1=D1;
- X Matrix MX=MD*M2*MD1 - D*(M2*D1); Print(MX);
- X MX=MD*M2*MD1 - (D*M2)*D1; Print(MX);
- X {
- X D.ReDimension(7); for (i=1;i<=7;i++) D(i,i)=i*3.0+i*i+2.0;
- X LowerTriangularMatrix LD(1); LD=D;
- X UpperTriangularMatrix UD(1); UD=D;
- X M2=U2; M2=LD*M2*MD1 - D*(U2*D1); Print(M2);
- X M2=U2; M2=UD*M2*MD1 - (D*U2)*D1; Print(M2);
- X M2=L2; M2=LD*M2*MD1 - D*(L2*D1); Print(M2);
- X M2=L2; M2=UD*M2*MD1 - (D*L2)*D1; Print(M2);
- X }
- X }
- X
- X// cout << "\nEnd of third test\n";
- X}
- X
- END_OF_FILE
- if test 2192 -ne `wc -c <'tmt3.cxx'`; then
- echo shar: \"'tmt3.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmt3.cxx'
- fi
- if test -f 'tmt4.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmt4.cxx'\"
- else
- echo shar: Extracting \"'tmt4.cxx'\" \(2436 characters\)
- sed "s/^X//" >'tmt4.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X
- X
- X#include "include.h"
- X
- X#include "newmat.h"
- X
- X
- X/**************************** test program ******************************/
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- Xvoid trymat4()
- X{
- X// cout << "\nFourth test of Matrix package\n";
- X Tracer et("Fourth test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X
- X int i,j;
- X
- X {
- X Tracer et1("Stage 1");
- X Matrix M(10,10);
- X UpperTriangularMatrix U(10);
- X for (i=1;i<=10;i++) for (j=1;j<=10;j++) M(i,j) = 100*i+j;
- X U << -M;
- X Matrix X1 = M.Rows(2,4);
- X Matrix Y1 = U.t().Rows(2,4);
- X Matrix X = U; { Print(Matrix(X.Columns(2,4).t()-Y1)); }
- X RowVector RV = M.Row(5);
- X {
- X X.ReDimension(3,10);
- X X.Row(1) << M.Row(2); X.Row(2) << M.Row(3); X.Row(3) << M.Row(4);
- X Print(Matrix(X-X1));
- X }
- X {
- X UpperTriangularMatrix V = U.SymSubMatrix(3,5);
- X Matrix MV = U.SubMatrix(3,5,3,5); { Print(Matrix(MV-V)); }
- X Matrix X2 = M.t().Columns(2,4); { Print(Matrix(X2-X1.t())); }
- X Matrix Y2 = U.Columns(2,4); { Print(Matrix(Y2-Y1.t())); }
- X ColumnVector CV = M.t().Column(5); { Print(ColumnVector(CV-RV.t())); }
- X X.ReDimension(10,3); M = M.t();
- X X.Column(1) << M.Column(2); X.Column(2) << M.Column(3);
- X X.Column(3) << M.Column(4);
- X Print(Matrix(X-X2));
- X }
- X }
- X
- X {
- X Tracer et1("Stage 2");
- X Matrix M; Matrix X; M.ReDimension(5,8);
- X for (i=1;i<=5;i++) for (j=1;j<=8;j++) M(i,j) = 100*i+j;
- X {
- X X = M.Columns(5,8); M.Columns(5,8) << M.Columns(1,4);
- X M.Columns(1,4) << X;
- X X = M.Columns(3,4); M.Columns(3,4) << M.Columns(1,2);
- X M.Columns(1,2) << X;
- X X = M.Columns(7,8); M.Columns(7,8) << M.Columns(5,6);
- X M.Columns(5,6) << X;
- X }
- X {
- X X = M.Column(2); M.Column(2) << M.Column(1); M.Column(1) << X;
- X X = M.Column(4); M.Column(4) << M.Column(3); M.Column(3) << X;
- X X = M.Column(6); M.Column(6) << M.Column(5); M.Column(5) << X;
- X X = M.Column(8); M.Column(8) << M.Column(7); M.Column(7) << X;
- X X.ReDimension(5,8);
- X }
- X for (i=1;i<=5;i++) for (j=1;j<=8;j++) X(i,9-j) = 100*i+j;
- X Print(Matrix(X-M));
- X }
- X
- X// cout << "\nEnd of fourth test\n";
- X}
- X
- END_OF_FILE
- if test 2436 -ne `wc -c <'tmt4.cxx'`; then
- echo shar: \"'tmt4.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmt4.cxx'
- fi
- if test -f 'tmt5.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmt5.cxx'\"
- else
- echo shar: Extracting \"'tmt5.cxx'\" \(1255 characters\)
- sed "s/^X//" >'tmt5.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X
- X
- X#include "include.h"
- X
- X#include "newmat.h"
- X
- X
- X/**************************** test program ******************************/
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- Xvoid trymat5()
- X{
- X// cout << "\nFifth test of Matrix package\n";
- X Tracer et("Fifth test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X
- X int i,j;
- X
- X Matrix A(5,6);
- X for (i=1;i<=5;i++) for (j=1;j<=6;j++) A(i,j)=1+i*j+i*i+j*j;
- X ColumnVector CV(6);
- X for (i=1;i<=6;i++) CV(i)=i*i+3;
- X ColumnVector CV2(5); for (i=1;i<=5;i++) CV2(i)=1.0;
- X ColumnVector CV1=CV;
- X
- X {
- X CV=A*CV;
- X RowVector RV=CV.t(); // RowVector RV; RV=CV.t();
- X RV=RV-1.0;
- X CV=(RV*A).t()+A.t()*CV2; CV1=(A.t()*A)*CV1 - CV;
- X Print(CV1);
- X }
- X
- X CV1.ReDimension(6);
- X CV2.ReDimension(6);
- X CV.ReDimension(6);
- X for (i=1;i<=6;i++) { CV1(i)=i*3+1; CV2(i)=10-i; CV(i)=11+i*2; }
- X ColumnVector CX=CV2-CV; { CX=CX+CV1; Print(CX); }
- X Print(ColumnVector(CV1+CV2-CV));
- X RowVector RV=CV.t(); RowVector RV1=CV1.t();
- X RowVector R=RV-RV1; Print(RowVector(R-CV2.t()));
- X// cout << "\nEnd of fifth test\n";
- X}
- END_OF_FILE
- if test 1255 -ne `wc -c <'tmt5.cxx'`; then
- echo shar: \"'tmt5.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmt5.cxx'
- fi
- if test -f 'tmt6.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmt6.cxx'\"
- else
- echo shar: Extracting \"'tmt6.cxx'\" \(1326 characters\)
- sed "s/^X//" >'tmt6.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X
- X#include "include.h"
- X
- X#include "newmat.h"
- X
- X
- X/**************************** test program ******************************/
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- Xvoid trymat6()
- X{
- X// cout << "\nSixth test of Matrix package\n";
- X Tracer et("Sixth test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X
- X int i,j;
- X
- X
- X DiagonalMatrix D(6);
- X UpperTriangularMatrix U(6);
- X for (i=1;i<=6;i++) { for (j=i;j<=6;j++) U(i,j)=i*i*i-50; D(i,i)=i*i+i-10; }
- X LowerTriangularMatrix L=(U*3.0).t();
- X SymmetricMatrix S(6);
- X for (i=1;i<=6;i++) for (j=i;j<=6;j++) S(i,j)=i*i+2.0+j;
- X Matrix MD=D; Matrix ML=L; Matrix MU=U; Matrix MS=S;
- X Matrix M(6,6);
- X for (i=1;i<=6;i++) for (j=1;j<=6;j++) M(i,j)=i*j+i*i-10.0;
- X {
- X Tracer et1("Stage 1");
- X Print(Matrix(MS+(-MS)));
- X Print(Matrix((S+M)-(MS+M)));
- X Print(Matrix((M+U)-(M+MU)));
- X Print(Matrix((M+L)-(M+ML)));
- X }
- X {
- X Tracer et1("Stage 2");
- X Print(Matrix((M+D)-(M+MD)));
- X Print(Matrix((U+D)-(MU+MD)));
- X Print(Matrix((D+L)-(ML+MD)));
- X Print(Matrix((-U+D)+MU-MD));
- X Print(Matrix((-L+D)+ML-MD));
- X }
- X
- X
- X
- X// cout << "\nEnd of sixth test\n";
- X}
- X
- END_OF_FILE
- if test 1326 -ne `wc -c <'tmt6.cxx'`; then
- echo shar: \"'tmt6.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmt6.cxx'
- fi
- if test -f 'tmt7.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmt7.cxx'\"
- else
- echo shar: Extracting \"'tmt7.cxx'\" \(1520 characters\)
- sed "s/^X//" >'tmt7.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X
- X#include "include.h"
- X
- X#include "newmat.h"
- X
- X
- X/**************************** test program ******************************/
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- Xvoid Clean(Matrix&, Real);
- X
- Xvoid trymat7()
- X{
- X// cout << "\nSeventh test of Matrix package\n";
- X Tracer et("Seventh test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X
- X int i,j;
- X
- X
- X DiagonalMatrix D(6);
- X UpperTriangularMatrix U(6);
- X 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; }
- X LowerTriangularMatrix L=(U*3.0).t();
- X SymmetricMatrix S(6);
- X for (i=1;i<=6;i++) for (j=i;j<=6;j++) S(i,j)=i*i+2.0+j;
- X Matrix MD=D; Matrix ML=L; Matrix MU=U;
- X Matrix MS=S;
- X Matrix M(6,6);
- X for (i=1;i<=6;i++) for (j=1;j<=6;j++) M(i,j)=i*j+i*i-10.0;
- X {
- X Tracer et1("Stage 1");
- X Print(Matrix((S-M)-(MS-M)));
- X Print(Matrix((-M-S)+(MS+M)));
- X Print(Matrix((U-M)-(MU-M)));
- X }
- X {
- X Tracer et1("Stage 2");
- X Print(Matrix((L-M)+(M-ML)));
- X Print(Matrix((D-M)+(M-MD)));
- X Print(Matrix((D-S)+(MS-MD)));
- X Print(Matrix((D-L)+(ML-MD)));
- X }
- X
- X { M=MU.t(); }
- X LowerTriangularMatrix LY=D.i()*U.t();
- X {
- X Tracer et1("Stage 3");
- X MS=D*LY-M; Clean(MS,0.00000001); Print(MS);
- X L=U.t();
- X LY=D.i()*L; MS=D*LY-M; Clean(MS,0.00000001); Print(MS);
- X }
- X
- X// cout << "\nEnd of seventh test\n";
- X}
- END_OF_FILE
- if test 1520 -ne `wc -c <'tmt7.cxx'`; then
- echo shar: \"'tmt7.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmt7.cxx'
- fi
- if test -f 'tmt8.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmt8.cxx'\"
- else
- echo shar: Extracting \"'tmt8.cxx'\" \(1414 characters\)
- sed "s/^X//" >'tmt8.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X
- X
- X#include "include.h"
- X
- X#include "newmat.h"
- X
- X
- X/**************************** test program ******************************/
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- Xvoid trymat8()
- X{
- X// cout << "\nEighth test of Matrix package\n";
- X Tracer et("Eighth test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X
- X int i;
- X
- X
- X DiagonalMatrix D(6);
- X for (i=1;i<=6;i++) D(i,i)=i*i+i-10;
- X DiagonalMatrix D2=D;
- X Matrix MD=D;
- X
- X DiagonalMatrix D1(6); for (i=1;i<=6;i++) D1(i,i)=-100+i*i*i;
- X Matrix MD1=D1;
- X Print(Matrix(D*D1-MD*MD1));
- X Print(Matrix((-D)*D1+MD*MD1));
- X Print(Matrix(D*(-D1)+MD*MD1));
- X DiagonalMatrix DX=D;
- X {
- X Tracer et1("Stage 1");
- X DX=(DX+D1)*DX; Print(Matrix(DX-(MD+MD1)*MD));
- X DX=D;
- X DX=-DX*DX+(DX-(-D1))*((-D1)+DX);
- X#ifdef GXX
- X Matrix MX = Matrix(MD1);
- X MD1=DX+(MX.t())*(MX.t()); Print(MD1);
- X#else
- X MD1=DX+(Matrix(MD1).t())*(Matrix(MD1).t()); Print(MD1);
- X#endif
- X DX=D; DX=DX; DX=D2-DX; Print(DiagonalMatrix(DX));
- X DX=D;
- X }
- X {
- X Tracer et1("Stage 2");
- X D.Release(2);
- X D1=D; D2=D;
- X Print(DiagonalMatrix(D1-DX));
- X Print(DiagonalMatrix(D2-DX));
- X MD1=1.0;
- X Print(Matrix(MD1-1.0));
- X }
- X
- X// cout << "\nEnd of eighth test\n";
- X}
- END_OF_FILE
- if test 1414 -ne `wc -c <'tmt8.cxx'`; then
- echo shar: \"'tmt8.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmt8.cxx'
- fi
- if test -f 'tmt9.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmt9.cxx'\"
- else
- echo shar: Extracting \"'tmt9.cxx'\" \(1577 characters\)
- sed "s/^X//" >'tmt9.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X
- X
- X#include "include.h"
- X#include "newmat.h"
- X
- X
- X/**************************** test program ******************************/
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- Xvoid Clean(Matrix&, Real);
- Xvoid Clean(DiagonalMatrix&, Real);
- X
- Xvoid trymat9()
- X{
- X// cout << "\nNinth test of Matrix package\n";
- X Tracer et("Ninth test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X
- X
- X int i; int j;
- X Matrix A(7,7); Matrix X(7,3);
- X for (i=1;i<=7;i++) for (j=1;j<=7;j++) A(i,j)=i*i+j+((i==j) ? 1 : 0);
- X for (i=1;i<=7;i++) for (j=1;j<=3;j++) X(i,j)=i-j;
- X Matrix B = A.i(); DiagonalMatrix D(7); D=1.0;
- X {
- X Tracer et1("Stage 1");
- X Matrix Q = B*A-D; Clean(Q, 0.000000001); Print(Q);
- X Q=A; Q = Q.i() * X; Q = A*Q - X; Clean(Q, 0.000000001); Print(Q);
- X Q=X; Q = A.i() * Q; Q = A*Q - X; Clean(Q, 0.000000001); Print(Q);
- X }
- X for (i=1;i<=7;i++) D(i,i)=i*i+1;
- X DiagonalMatrix E(3); for (i=1;i<=3;i++) E(i,i)=i+23;
- X {
- X Tracer et1("Stage 2");
- X Matrix DXE = D.i() * X * E;
- X DXE = E.i() * DXE.t() * D - X.t(); Clean(DXE, 0.00000001); Print(DXE);
- X E=D; for (i=1;i<=7;i++) E(i,i)=i*3+1;
- X }
- X DiagonalMatrix F=D;
- X {
- X Tracer et1("Stage 3");
- X F=E.i()*F; F=F*E-D; Clean(F,0.00000001); Print(F);
- X F=E.i()*D; F=F*E-D; Clean(F,0.00000001); Print(F);
- X }
- X { F=E; F=F.i()*D; F=F*E-D; Clean(F,0.00000001); Print(F); }
- X// cout << "\nEnd of ninth test\n";
- X}
- X
- END_OF_FILE
- if test 1577 -ne `wc -c <'tmt9.cxx'`; then
- echo shar: \"'tmt9.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmt9.cxx'
- fi
- if test -f 'tmta.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmta.cxx'\"
- else
- echo shar: Extracting \"'tmta.cxx'\" \(2837 characters\)
- sed "s/^X//" >'tmta.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X
- X
- X#include "include.h"
- X
- X#include "newmatap.h"
- X
- X
- X/**************************** test program ******************************/
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- Xvoid Clean(Matrix&, Real);
- X
- X
- Xstatic void process(const GeneralMatrix& A,
- X const ColumnVector& X1, const ColumnVector& X2)
- X{
- X Matrix B = A;
- X LinearEquationSolver L(A);
- X Matrix Y(4,2);
- X Y.Column(1) << L.i() * X1; Y.Column(2) << L.i() * X2;
- X Matrix Z(4,2); Z.Column(1) << X1; Z.Column(2) << X2;
- X Z = B * Y - Z; Clean(Z,0.00000001); Print(Z);
- X}
- X
- X
- X
- Xvoid trymata()
- X{
- X// cout << "\nTenth test of Matrix package\n";
- X Tracer et("Tenth test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X int i; int j;
- X UpperTriangularMatrix U(8);
- X for (i=1;i<=8;i++) for (j=i;j<=8;j++) U(i,j)=i+j*j+5;
- X Matrix X(8,6);
- X for (i=1;i<=8;i++) for (j=1;j<=6;j++) X(i,j)=i*j+1.0;
- X Matrix Y = U.i()*X; Matrix MU=U;
- X Y=Y-MU.i()*X; Clean(Y,0.00000001); Print(Y);
- X Y = U.t().i()*X; Y=Y-MU.t().i()*X; Clean(Y,0.00000001); Print(Y);
- X UpperTriangularMatrix UX(8);
- X for (i=1;i<=8;i++) for (j=i;j<=8;j++) UX(i,j)=i+j+1;
- X UX(4,4)=0; UX(4,5)=0;
- X UpperTriangularMatrix UY = U.i() * UX;
- X { X=UX; MU=U; Y=UY-MU.i()*X; Clean(Y,0.000000001); Print(Y); }
- X LowerTriangularMatrix LY = U.t().i() * UX.t();
- X { Y=LY-MU.i().t()*X.t(); Clean(Y,0.000000001); Print(Y); }
- X DiagonalMatrix D(8); for (i=1;i<=8;i++) D(i,i)=i+1;
- X { X=D.i()*MU; }
- X { UY=U; UY=D.i()*UY; Y=UY-X; Clean(Y,0.00000001); Print(Y); }
- X { UY=D.i()*U; Y=UY-X; Clean(Y,0.00000001); Print(Y); }
- X// X=MU.t();
- X// LY=D.i()*U.t(); Y=D*LY-X; Clean(Y,0.00000001); Print(Y);
- X// LowerTriangularMatrix L=U.t();
- X// LY=D.i()*L; Y=D*LY-X; Clean(Y,0.00000001); Print(Y);
- X U.ReDimension(8);
- X for (i=1;i<=8;i++) for (j=i;j<=8;j++) U(i,j)=i+j*j+5;
- X MU = U;
- X MU = U.i() - MU.i(); Clean(MU,0.00000001); Print(MU);
- X MU = U.t().i() - U.i().t(); Clean(MU,0.00000001); Print(MU);
- X
- X // test LINEQ
- X {
- X ColumnVector X1(4), X2(4);
- X X1(1)=1; X1(2)=2; X1(3)=3; X1(4)=4;
- X X2(1)=1; X2(2)=10; X2(3)=100; X2(4)=1000;
- X
- X
- X Matrix A(4,4);
- X A(1,1)=1; A(1,2)=3; A(1,3)=0; A(1,4)=0;
- X A(2,1)=3; A(2,2)=2; A(2,3)=5; A(2,4)=0;
- X A(3,1)=0; A(3,2)=5; A(3,3)=4; A(3,4)=1;
- X A(4,1)=0; A(4,2)=0; A(4,3)=1; A(4,4)=3;
- X process(A,X1,X2);
- X
- X BandMatrix B(4,1,1); B.Inject(A);
- X process(B,X1,X2);
- X
- X UpperTriangularMatrix U(4);
- X U(1,1)=1; U(1,2)=2; U(1,3)=3; U(1,4)=4;
- X U(2,2)=8; U(2,3)=7; U(2,4)=6;
- X U(3,3)=2; U(3,4)=7;
- X U(4,4)=1;
- X process(U,X1,X2);
- X
- X LowerTriangularMatrix L = U.t();
- X process(L,X1,X2);
- X }
- X
- X// cout << "\nEnd of tenth test\n";
- X}
- END_OF_FILE
- if test 2837 -ne `wc -c <'tmta.cxx'`; then
- echo shar: \"'tmta.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmta.cxx'
- fi
- if test -f 'tmtb.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmtb.cxx'\"
- else
- echo shar: Extracting \"'tmtb.cxx'\" \(4020 characters\)
- sed "s/^X//" >'tmtb.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X
- X#include "include.h"
- X
- X#include "newmat.h"
- X
- X
- X/**************************** test program ******************************/
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- Xvoid Clean(Matrix&, Real);
- X
- Xvoid trymatb()
- X{
- X// cout << "\nEleventh test of Matrix package\n";
- X Tracer et("Eleventh test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X int i; int j;
- X RowVector RV; RV.ReDimension(10);
- X {
- X Tracer et1("Stage 1");
- X for (i=1;i<=10;i++) RV(i)=i*i-3;
- X Matrix X(1,1); X(1,1) = .25;
- X Print(RowVector(X.i() * RV - RV / .25));
- X// Print(RowVector(X.i() * Matrix(RV) - RV / .25)); // != zortech, AT&T
- X Print(RowVector(X.i() * RV - RV / .25));
- X }
- X LowerTriangularMatrix L(5); UpperTriangularMatrix U(5);
- X for (i=1; i<=5; i++) for (j=1; j<=i; j++)
- X { L(i,j) = i*i + j -2.0; U(j,i) = i*i*j+3; }
- X DiagonalMatrix D(5);
- X for (i=1; i<=5; i++) D(i,i) = i*i + i + 2;
- X Matrix M1 = -L; Matrix M2 = L-U; Matrix M3 = U*3; Matrix M4 = U-L;
- X Matrix M5 = M1 - D; M1 = D * (-3) - M3;
- X {
- X Tracer et1("Stage 2");
- X Print(Matrix((M2-M4*2)+M5*3-M1));
- X M1 = L.t(); Print(Matrix(M1.t()-L));
- X M1 = U.t(); Print(Matrix(M1.t()-U));
- X }
- X {
- X Tracer et1("Stage 3");
- X SymmetricMatrix S(5);
- X for (i=1; i<=5; i++) for (j=1; j<=i; j++) S(i,j) = i*j+i-j+5;
- X M2 = S.i() * M4; M1 = S; M3=M1*M2-M4; Clean(M3,0.00000001); Print(M3);
- X SymmetricMatrix T(5);
- X for (i=1; i<=5; i++) for (j=1; j<=i; j++) T(i,j) = i*i*j*j+i-j+5-i*j;
- X M1 = S.i() * T; M1 = S * M1; M1 = M1-T; Clean(M1,0.00000001); Print(M1);
- X ColumnVector CV(5); for (i=1; i<=5; i++) CV(i) = i*i*i+10;
- X M1 = CV * RV;
- X }
- X {
- X Tracer et1("Stage 4");
- X M4.ReDimension(5,10);
- X for (i=1; i<=5; i++) for (j=1; j<=10; j++) M4(i,j) = (i*i*i+10)*(j*j-3);
- X Print(Matrix(M1-M4));
- X M1 = L.t(); M2 = U.t(); M3 = L+U; Print(Matrix(M1-M3.t()+M2));
- X }
- X// UpperTriangularMatrix U2((const UpperTriangularMatrix&)U); // != zortech
- X UpperTriangularMatrix U2((UpperTriangularMatrix&)U);
- X {
- X Tracer et1("Stage 5");
- X Print(Matrix(U2-U));
- X M2.ReDimension(10,10);
- X for (i=1; i<=10; i++) for (j=1; j<=10; j++) M2(i,j) = (i*i*i+10)*(j*j-3);
- X D << M2; L << M2; U << M2; // check copy into
- X Print( Matrix( (D+M2)-(L+U) ) );
- X }
- X {
- X Tracer et1("Stage 6");
- X M1.ReDimension(6,10);
- X for (i=1; i<=6; i++) for (j=1; j<=10; j++) M1(i,j) = 100*i + j;
- X M2 = M1.SubMatrix(3,5,4,7); M3.ReDimension(3,4);
- X for (i=3; i<=5; i++) for (j=4; j<=7; j++) M3(i-2,j-3) = 100*i + j;
- X Print(Matrix(M2-M3));
- X }
- X int a1,a2,a3,a4;
- X {
- X Tracer et1("Stage 7");
- X int a1,a2,a3,a4;
- X a1=4; a2=9; a3=4; a4=7;
- X U.ReDimension(10);
- X for (i=1; i<=10; i++) for (j=i; j<=10; j++) U(i,j) = 100*i + j;
- X M2 = U.SubMatrix(a1,a2,a3,a4);
- X M3.ReDimension(a2-a1+1,a4-a3+1); M3=0.0;
- X for (i=a1; i<=a2; i++) for (j=(i>a3) ? i : a3; j<=a4; j++)
- X M3(i-a1+1,j-a3+1) = 100*i + j;
- X Print(Matrix(M2-M3));
- X }
- X {
- X Tracer et1("Stage 8");
- X a1=3; a2=9; a3=2; a4=7;
- X U.ReDimension(10);
- X for (i=1; i<=10; i++) for (j=i; j<=10; j++) U(i,j) = 100*i + j;
- X M2 = U.SubMatrix(a1,a2,a3,a4);
- X M3.ReDimension(a2-a1+1,a4-a3+1); M3=0.0;
- X for (i=a1; i<=a2; i++) for (j=(i>a3) ? i : a3; j<=a4; j++)
- X M3(i-a1+1,j-a3+1) = 100*i + j;
- X Print(Matrix(M2-M3));
- X }
- X {
- X Tracer et1("Stage 9");
- X a1=4; a2=6; a3=2; a4=5;
- X U.ReDimension(10);
- X for (i=1; i<=10; i++) for (j=i; j<=10; j++) U(i,j) = 100*i + j;
- X M2 = U.SubMatrix(a1,a2,a3,a4);
- X M3.ReDimension(a2-a1+1,a4-a3+1); M3=0.0;
- X for (i=a1; i<=a2; i++) for (j=(i>a3) ? i : a3; j<=a4; j++)
- X M3(i-a1+1,j-a3+1) = 100*i + j;
- X Print(Matrix(M2-M3));
- X }
- X
- X// cout << "\nEnd of eleventh test\n";
- X}
- END_OF_FILE
- if test 4020 -ne `wc -c <'tmtb.cxx'`; then
- echo shar: \"'tmtb.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmtb.cxx'
- fi
- if test -f 'tmtc.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmtc.cxx'\"
- else
- echo shar: Extracting \"'tmtc.cxx'\" \(2900 characters\)
- sed "s/^X//" >'tmtc.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X
- X
- X#include "include.h"
- X#include "newmat.h"
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- Xvoid Clean(Matrix&, Real);
- X
- X
- X
- Xvoid trymatc()
- X{
- X// cout << "\nTwelfth test of Matrix package\n";
- X Tracer et("Twelfth test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X DiagonalMatrix D(15); D=1.5;
- X Matrix A(15,15);
- X int i,j;
- X for (i=1;i<=15;i++) for (j=1;j<=15;j++) A(i,j)=i*i+j-150;
- X { A = A + D; }
- X ColumnVector B(15);
- X for (i=1;i<=15;i++) B(i)=i+i*i-150.0;
- X {
- X Tracer et1("Stage 1");
- X ColumnVector B1=B;
- X B=(A*2.0).i() * B1;
- X Matrix X = A*B-B1/2.0;
- X Clean(X, 0.000000001); Print(X);
- X A.ReDimension(3,5);
- X for (i=1; i<=3; i++) for (j=1; j<=5; j++) A(i,j) = i+100*j;
- X
- X B = A.AsColumn()+10000;
- X RowVector R = (A+10000).AsColumn().t();
- X Print( RowVector(R-B.t()) );
- X }
- X
- X {
- X Tracer et1("Stage 2");
- X B = A.AsColumn()+10000;
- X Matrix XR = (A+10000).AsMatrix(15,1).t();
- X Print( RowVector(XR-B.t()) );
- X }
- X
- X {
- X Tracer et1("Stage 3");
- X B = (A.AsMatrix(15,1)+A.AsColumn())/2.0+10000;
- X Matrix MR = (A+10000).AsColumn().t();
- X Print( RowVector(MR-B.t()) );
- X
- X B = (A.AsMatrix(15,1)+A.AsColumn())/2.0;
- X MR = A.AsColumn().t();
- X Print( RowVector(MR-B.t()) );
- X }
- X
- X {
- X Tracer et1("Stage 4");
- X B = (A.AsMatrix(15,1)+A.AsColumn())/2.0;
- X RowVector R = A.AsColumn().t();
- X Print( RowVector(R-B.t()) );
- X }
- X
- X {
- X Tracer et1("Stage 5");
- X RowVector R = (A.AsColumn()-5000).t();
- X B = ((R.t()+10000) - A.AsColumn())-5000;
- X Print( RowVector(B.t()) );
- X }
- X
- X {
- X Tracer et1("Stage 6");
- X B = A.AsColumn(); ColumnVector B1 = (A+10000).AsColumn() - 10000;
- X Print(ColumnVector(B1-B));
- X }
- X
- X {
- X Tracer et1("Stage 7");
- X Matrix X = B.AsMatrix(3,5); Print(Matrix(X-A));
- X for (i=1; i<=3; i++) for (j=1; j<=5; j++) B(5*(i-1)+j) -= i+100*j;
- X Print(B);
- X }
- X
- X {
- X Tracer et1("Stage 8");
- X A.ReDimension(7,7); D.ReDimension(7);
- X for (i=1; i<=7; i++) for (j=1; j<=7; j++) A(i,j) = i*j*j;
- X for (i=1; i<=7; i++) D(i,i) = i;
- X UpperTriangularMatrix U; U << A;
- X Matrix X = A; for (i=1; i<=7; i++) X(i,i) = i;
- X A.Inject(D); Print(Matrix(X-A));
- X X = U; U.Inject(D); A = U; for (i=1; i<=7; i++) X(i,i) = i;
- X Print(Matrix(X-A));
- X }
- X
- X {
- X Tracer et1("Stage 9");
- X A.ReDimension(7,5);
- X for (i=1; i<=7; i++) for (j=1; j<=5; j++) A(i,j) = i+100*j;
- X Matrix Y = A; Y = Y - ((const Matrix&)A); Print(Y);
- X Matrix X = A; // X.Release();
- X Y = A; Y = ((const Matrix&)X) - A; Print(Y); Y = 0.0;
- X Y = ((const Matrix&)X) - ((const Matrix&)A); Print(Y);
- X }
- X
- X// cout << "\nEnd of twelfth test\n";
- X}
- END_OF_FILE
- if test 2900 -ne `wc -c <'tmtc.cxx'`; then
- echo shar: \"'tmtc.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmtc.cxx'
- fi
- if test -f 'tmtd.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmtd.cxx'\"
- else
- echo shar: Extracting \"'tmtd.cxx'\" \(3508 characters\)
- sed "s/^X//" >'tmtd.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X
- X#include "include.h"
- X#include "newmatap.h"
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- Xvoid Clean(Matrix&, Real);
- X
- X
- X
- X
- Xvoid trymatd()
- X{
- X// cout << "\nThirteenth test of Matrix package\n";
- X Tracer et("Thirteenth test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X Matrix X(5,20);
- X int i,j;
- X for (j=1;j<=20;j++) X(1,j) = j+1;
- X for (i=2;i<=5;i++) for (j=1;j<=20; j++) X(i,j) = (long)X(i-1,j) * j % 1001;
- X SymmetricMatrix S; S << X * X.t();
- X LowerTriangularMatrix L = Cholesky(S);
- X Matrix Diff = L*L.t()-S; Clean(Diff, 0.000000001);
- X Print(Diff);
- X LowerTriangularMatrix L1(5);
- X HHDecompose(X,L1);
- X Diff = L - L1; Clean(Diff,0.000000001); Print(Diff);
- X
- X ColumnVector C1(5);
- X {
- X Tracer et1("Stage 1");
- X X.ReDimension(5,5);
- X for (j=1;j<=5;j++) X(1,j) = j+1;
- X for (i=2;i<=5;i++) for (j=1;j<=5; j++)
- X X(i,j) = (long)X(i-1,j) * j % 1001;
- X for (i=1;i<=5;i++) C1(i) = i*i;
- X CroutMatrix A = X;
- X ColumnVector C2 = A.i() * C1; C1 = X.i() * C1;
- X X = C1 - C2; Clean(X,0.000000001); Print(X);
- X }
- X
- X {
- X Tracer et1("Stage 2");
- X X.ReDimension(7,7);
- X for (j=1;j<=7;j++) X(1,j) = j+1;
- X for (i=2;i<=7;i++) for (j=1;j<=7; j++)
- X X(i,j) = (long)X(i-1,j) * j % 1001;
- X C1.ReDimension(7);
- X for (i=1;i<=7;i++) C1(i) = i*i;
- X RowVector R1 = C1.t();
- X Diff = R1 * X.i() - ( X.t().i() * R1.t() ).t(); Clean(Diff,0.000000001);
- X Print(Diff);
- X }
- X
- X {
- X Tracer et1("Stage 3");
- X X.ReDimension(5,5);
- X for (j=1;j<=5;j++) X(1,j) = j+1;
- X for (i=2;i<=5;i++) for (j=1;j<=5; j++)
- X X(i,j) = (long)X(i-1,j) * j % 1001;
- X C1.ReDimension(5);
- X for (i=1;i<=5;i++) C1(i) = i*i;
- X CroutMatrix A1 = X*X;
- X ColumnVector C2 = A1.i() * C1; C1 = X.i() * C1; C1 = X.i() * C1;
- X X = C1 - C2; Clean(X,0.000000001); Print(X);
- X }
- X
- X
- X {
- X Tracer et1("Stage 4");
- X int n = 40;
- X SymmetricBandMatrix B(n,2); B = 0.0;
- X for (i=1; i<=n; i++)
- X {
- X B(i,i) = 6;
- X if (i<=n-1) B(i,i+1) = -4;
- X if (i<=n-2) B(i,i+2) = 1;
- X }
- X B(1,1) = 5; B(n,n) = 5;
- X SymmetricMatrix A = B;
- X ColumnVector X(n); X = 0; X(1) = 1;
- X ColumnVector Y1 = A.i() * X;
- X LowerTriangularMatrix C1 = Cholesky(A);
- X ColumnVector Y2 = C1.t().i() * (C1.i() * X) - Y1;
- X Clean(Y2, 0.000000001); Print(Y2);
- X
- X LowerBandMatrix C2 = Cholesky(B);
- X Matrix M = C2 - C1; Clean(M, 0.000000001); Print(M);
- X ColumnVector Y3 = C2.t().i() * (C2.i() * X) - Y1;
- X Clean(Y3, 0.000000001); Print(Y3);
- X }
- X
- X {
- X Tracer et1("Stage 5");
- X SymmetricMatrix A(4), B(4);
- X#ifdef ATandT
- X Real a[10], b[10];
- X a[0] = 5; a[1] = 1; a[2] = 4; a[3] = 2; a[4] = 1;
- X a[5] = 6; a[6] = 1; a[7] = 0; a[8] = 1; a[9] = 7;
- X b[0] = 8; b[1] = 1; b[2] = 5; b[3] = 1; b[4] = 0;
- X b[5] = 9; b[6] = 2; b[7] = 1; b[8] = 0; b[9] = 6;
- X#else
- X Real a[] = { 5, 1, 4, 2, 1, 6, 1, 0, 1, 7 };
- X Real b[] = { 8, 1, 5, 1, 0, 9, 2, 1, 0, 6 };
- X#endif
- X A << a; B << b;
- X LowerTriangularMatrix AB = Cholesky(A) * Cholesky(B);
- X Matrix M = Cholesky(A) * B * Cholesky(A).t() - AB*AB.t();
- X Clean(M, 0.000000001); Print(M);
- X M = A * Cholesky(B); M = M * M.t() - A * B * A;
- X Clean(M, 0.000000001); Print(M);
- X }
- X
- X
- X// cout << "\nEnd of Thirteenth test\n";
- X}
- END_OF_FILE
- if test 3508 -ne `wc -c <'tmtd.cxx'`; then
- echo shar: \"'tmtd.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmtd.cxx'
- fi
- if test -f 'tmte.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmte.cxx'\"
- else
- echo shar: Extracting \"'tmte.cxx'\" \(5726 characters\)
- sed "s/^X//" >'tmte.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X#define WANT_MATH
- X
- X#include "include.h"
- X
- X#include "newmatap.h"
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- Xvoid Clean(Matrix&, Real);
- Xvoid Clean(DiagonalMatrix&, Real);
- X
- X
- X
- X
- Xvoid trymate()
- X{
- X// cout << "\nFourteenth test of Matrix package\n";
- X Tracer et("Fourteenth test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X
- X {
- X Tracer et1("Stage 1");
- X Matrix A(8,5);
- X#ifndef ATandT
- X Real a[] = { 22, 10, 2, 3, 7,
- X 14, 7, 10, 0, 8,
- X -1, 13, -1,-11, 3,
- X -3, -2, 13, -2, 4,
- X 9, 8, 1, -2, 4,
- X 9, 1, -7, 5, -1,
- X 2, -6, 6, 5, 1,
- X 4, 5, 0, -2, 2 };
- X#else
- X Real a[40];
- X a[ 0]=22; a[ 1]=10; a[ 2]= 2; a[ 3]= 3; a[ 4]= 7;
- X a[ 5]=14; a[ 6]= 7; a[ 7]=10; a[ 8]= 0; a[ 9]= 8;
- X a[10]=-1; a[11]=13; a[12]=-1; a[13]=-11;a[14]= 3;
- X a[15]=-3; a[16]=-2; a[17]=13; a[18]=-2; a[19]= 4;
- X a[20]= 9; a[21]= 8; a[22]= 1; a[23]=-2; a[24]= 4;
- X a[25]= 9; a[26]= 1; a[27]=-7; a[28]= 5; a[29]=-1;
- X a[30]= 2; a[31]=-6; a[32]= 6; a[33]= 5; a[34]= 1;
- X a[35]= 4; a[36]= 5; a[37]= 0; a[38]=-2; a[39]= 2;
- X#endif
- X A << a;
- X DiagonalMatrix D; Matrix U; Matrix V;
- X#ifdef ATandT
- X int anc = A.Ncols(); DiagonalMatrix I(anc); // AT&T 2.1 bug
- X#else
- X DiagonalMatrix I(A.Ncols());
- X#endif
- X I=1.0;
- X SymmetricMatrix S1; S1 << A.t() * A;
- X SymmetricMatrix S2; S2 << A * A.t();
- X Real zero = 0.0; SVD(A+zero,D,U,V);
- X DiagonalMatrix D1; SVD(A,D1); Print(DiagonalMatrix(D-D1));
- X Matrix SU = U.t() * U - I; Clean(SU,0.000000001); Print(SU);
- X Matrix SV = V.t() * V - I; Clean(SV,0.000000001); Print(SV);
- X Matrix B = U * D * V.t() - A; Clean(B,0.000000001);Print(B);
- X D1=0.0; SVD(A,D1,A); Print(Matrix(A-U));
- X SortDescending(D);
- X D(1) -= sqrt(1248); D(2) -= 20; D(3) -= sqrt(384);
- X Clean(D,0.000000001); Print(D);
- X Jacobi(S1, D, V);
- X V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
- X SortDescending(D); D(1)-=1248; D(2)-=400; D(3)-=384;
- X Clean(D,0.000000001); Print(D);
- X Jacobi(S2, D, V);
- X V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
- X SortDescending(D); D(1)-=1248; D(2)-=400; D(3)-=384;
- X Clean(D,0.000000001); Print(D);
- X
- X EigenValues(S1, D, V);
- X V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
- X D(5)-=1248; D(4)-=400; D(3)-=384;
- X Clean(D,0.000000001); Print(D);
- X EigenValues(S2, D, V);
- X V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
- X D(8)-=1248; D(7)-=400; D(6)-=384;
- X Clean(D,0.000000001); Print(D);
- X
- X EigenValues(S1, D);
- X D(5)-=1248; D(4)-=400; D(3)-=384;
- X Clean(D,0.000000001); Print(D);
- X EigenValues(S2, D);
- X D(8)-=1248; D(7)-=400; D(6)-=384;
- X Clean(D,0.000000001); Print(D);
- X }
- X {
- X Tracer et1("Stage 2");
- X Matrix A(20,21);
- X for (int i=1; i<=20; i++) for (int j=1; j<=21; j++)
- X { if (i>j) A(i,j) = 0; else if (i==j) A(i,j) = 21-i; else A(i,j) = -1; }
- X A = A.t();
- X SymmetricMatrix S1; S1 << A.t() * A;
- X SymmetricMatrix S2; S2 << A * A.t();
- X DiagonalMatrix D; Matrix U; Matrix V;
- X#ifdef ATandT
- X int anc = A.Ncols(); DiagonalMatrix I(anc); // AT&T 2.1 bug
- X#else
- X DiagonalMatrix I(A.Ncols());
- X#endif
- X I=1.0;
- X SVD(A,D,U,V);
- X Matrix SU = U.t() * U - I; Clean(SU,0.000000001); Print(SU);
- X Matrix SV = V.t() * V - I; Clean(SV,0.000000001); Print(SV);
- X Matrix B = U * D * V.t() - A; Clean(B,0.000000001); Print(B);
- X for (i=1; i<=20; i++) D(i) -= sqrt((22-i)*(21-i));
- X Clean(D,0.000000001); Print(D);
- X Jacobi(S1, D, V);
- X V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
- X SortDescending(D);
- X for (i=1; i<=20; i++) D(i) -= (22-i)*(21-i);
- X Clean(D,0.000000001); Print(D);
- X Jacobi(S2, D, V);
- X V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
- X SortDescending(D);
- X for (i=1; i<=20; i++) D(i) -= (22-i)*(21-i);
- X Clean(D,0.000000001); Print(D);
- X
- X EigenValues(S1, D, V);
- X V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
- X for (i=1; i<=20; i++) D(i) -= (i+1)*i;
- X Clean(D,0.000000001); Print(D);
- X EigenValues(S2, D, V);
- X V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
- X for (i=2; i<=21; i++) D(i) -= (i-1)*i;
- X Clean(D,0.000000001); Print(D);
- X
- X EigenValues(S1, D);
- X for (i=1; i<=20; i++) D(i) -= (i+1)*i;
- X Clean(D,0.000000001); Print(D);
- X EigenValues(S2, D);
- X for (i=2; i<=21; i++) D(i) -= (i-1)*i;
- X Clean(D,0.000000001); Print(D);
- X }
- X
- X {
- X Tracer et1("Stage 3");
- X Matrix A(30,30);
- X for (int i=1; i<=30; i++) for (int j=1; j<=30; j++)
- X { if (i>j) A(i,j) = 0; else if (i==j) A(i,j) = 1; else A(i,j) = -1; }
- X Real d1 = A.LogDeterminant().Value();
- X DiagonalMatrix D; Matrix U; Matrix V;
- X#ifdef ATandT
- X int anc = A.Ncols(); DiagonalMatrix I(anc); // AT&T 2.1 bug
- X#else
- X DiagonalMatrix I(A.Ncols());
- X#endif
- X I=1.0;
- X SVD(A,D,U,V);
- X Matrix SU = U.t() * U - I; Clean(SU,0.000000001); Print(SU);
- X Matrix SV = V.t() * V - I; Clean(SV,0.000000001); Print(SV);
- X Real d2 = D.LogDeterminant().Value();
- X Matrix B = U * D * V.t() - A; Clean(B,0.000000001); Print(B);
- X SortDescending(D); // Print(D);
- X Real d3 = D.LogDeterminant().Value();
- X ColumnVector Test(3);
- X Test(1) = d1 - 1; Test(2) = d2 - 1; Test(3) = d3 - 1;
- X Clean(Test,0.00000001); Print(Test); // only 8 decimal figures
- X }
- X
- X// cout << "\nEnd of Fourteenth test\n";
- X}
- END_OF_FILE
- if test 5726 -ne `wc -c <'tmte.cxx'`; then
- echo shar: \"'tmte.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmte.cxx'
- fi
- if test -f 'tmtf.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmtf.cxx'\"
- else
- echo shar: Extracting \"'tmtf.cxx'\" \(1681 characters\)
- sed "s/^X//" >'tmtf.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X#define WANT_MATH
- X
- X#include "include.h"
- X
- X#include "newmatap.h"
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- Xvoid Clean(Matrix&, Real);
- X
- X
- X
- X
- Xstatic void test(int n)
- X{
- X Tracer et("Test FFT");
- X
- X ColumnVector A(n), B(n), X, Y;
- X for (int i=0; i<n; i++)
- X {
- X Real f = 6.2831853071795864769*i/n;
- X A.element(i) = fabs(sin(7.0*f) + 0.5 * cos(19.0 * f));
- X B.element(i) = fabs(0.25 * cos(31.0 * f));
- X }
- X FFT(A, B, X, Y);
- X X=X/n; Y=-Y/n;
- X FFT(X, Y, X, Y);
- X Y = -Y;
- X X = X - A; Y = Y - B;
- X Clean(X,0.000000001); Clean(Y,0.000000001); Print(X); Print(Y);
- X}
- X
- Xvoid trymatf()
- X{
- X Tracer et("Fifteenth test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X// cout << "\nFifteenth test of Matrix package\n";
- X// cout << "\n";
- X
- X ColumnVector A(12), B(12);
- X for (int i = 1; i <=12; i++)
- X {
- X Real i1 = i - 1;
- X A(i) = .7
- X + .2 * cos(6.2831853071795864769 * 4.0 * i1 / 12)
- X + .3 * sin(6.2831853071795864769 * 3.0 * i1 / 12);
- X B(i) = .9
- X + .5 * sin(6.2831853071795864769 * 2.0 * i1 / 12)
- X + .4 * cos(6.2831853071795864769 * 1.0 * i1 / 12);
- X }
- X FFT(A, B, A, B);
- X A(1) -= 8.4; A(3) -= 3.0; A(5) -= 1.2; A(9) -= 1.2; A(11) += 3.0;
- X B(1) -= 10.8; B(2) -= 2.4; B(4) += 1.8; B(10) -= 1.8; B(12) -= 2.4;
- X Clean(A,0.000000001); Clean(B,0.000000001); Print(A); Print(B);
- X
- X test(2048);
- X test(2000);
- X test(27*81);
- X test(2310);
- X test(49*49);
- X test(13*13*13);
- X test(43*47);
- X
- X// cout << "\nEnd of Fifteenth test\n";
- X}
- END_OF_FILE
- if test 1681 -ne `wc -c <'tmtf.cxx'`; then
- echo shar: \"'tmtf.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmtf.cxx'
- fi
- if test -f 'tmtg.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmtg.cxx'\"
- else
- echo shar: Extracting \"'tmtg.cxx'\" \(4734 characters\)
- sed "s/^X//" >'tmtg.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X
- X#include "include.h"
- X
- X#include "newmatap.h"
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- Xvoid Clean(Matrix&, Real);
- X
- X
- X
- X
- Xvoid trymatg()
- X{
- X// cout << "\nSixteenth test of Matrix package\n";
- X// cout << "\n";
- X Tracer et("Sixteenth test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X
- X int i,j;
- X Matrix M(4,7);
- X for (i=1; i<=4; i++) for (j=1; j<=7; j++) M(i,j) = 100 * i + j;
- X ColumnVector CV = M.AsColumn();
- X {
- X Tracer et1("Stage 1");
- X RowVector test(7);
- X test(1) = SumSquare(M);
- X test(2) = SumSquare(CV);
- X test(3) = SumSquare(CV.t());
- X test(4) = SumSquare(CV.AsDiagonal());
- X test(5) = SumSquare(M.AsColumn());
- X test(6) = Matrix(CV.t()*CV)(1,1);
- X test(7) = (CV.t()*CV).AsScalar();
- X test = test - 2156560; Print(test);
- X }
- X
- X/*
- X MatrixType mt(MatrixType::US);
- X {
- X M.AsColumn().MatrixParameters(i,j,mt);
- X cout << (char*)mt <<" "<< i <<" "<< j <<" (ColV, 28, 1)\n";
- X }
- X {
- X M.AsRow().MatrixParameters(i,j,mt);
- X cout << (char*)mt <<" "<< i <<" "<< j <<" (RowV, 1, 28)\n";
- X }
- X {
- X M.AsDiagonal().MatrixParameters(i,j,mt);
- X cout << (char*)mt <<" "<< i <<" "<< j <<" (Diag, 28, 28)\n";
- X }
- X {
- X M.AsMatrix(2,14).MatrixParameters(i,j,mt);
- X cout << (char*)mt <<" "<< i <<" "<< j <<" (Rect, 2, 14)\n";
- X }
- X {
- X (M.t()*M*M.t()).MatrixParameters(i,j,mt);
- X cout << (char*)mt <<" "<< i <<" "<< j <<" (Rect, 7, 4)\n";
- X }
- X {
- X M.t().MatrixParameters(i,j,mt);
- X cout << (char*)mt <<" "<< i <<" "<< j <<" (Rect, 7, 4)\n";
- X }
- X {
- X (M/2.0).MatrixParameters(i,j,mt);
- X cout << (char*)mt <<" "<< i <<" "<< j <<" (Rect, 4, 7)\n";
- X }
- X*/
- X UpperTriangularMatrix U(6);
- X for (i=1; i<=6; i++) for (j=i; j<=6; j++) U(i,j) = i + (i-j) * (i-j);
- X M = U; DiagonalMatrix D; D << U;
- X LowerTriangularMatrix L = U.t(); SymmetricMatrix S; S << (L+U)/2.0;
- X {
- X Tracer et1("Stage 2");
- X RowVector test(6);
- X test(1) = U.Trace();
- X test(2) = L.Trace();
- X test(3) = D.Trace();
- X test(4) = S.Trace();
- X test(5) = M.Trace();
- X test(6) = ((L+U)/2.0).Trace();
- X test = test - 21; Print(test);
- X test.ReDimension(5);
- X test(1) = LogDeterminant(U).Value();
- X test(2) = LogDeterminant(L).Value();
- X test(3) = LogDeterminant(D).Value();
- X test(4) = LogDeterminant(D).Value();
- X test(5) = LogDeterminant((L+D)/2.0).Value(); // (!=Glockenspiel)
- X test = test - 720; Clean(test,0.000000001); Print(test);
- X }
- X
- X {
- X Tracer et1("Stage 3");
- X S << L*U; M = S;
- X RowVector test(4);
- X test(1) = LogDeterminant(S).Value();
- X test(2) = LogDeterminant(M).Value();
- X test(3) = LogDeterminant(L*U).Value(); // (!=Glockenspiel)
- X test(4) = LogDeterminant(Matrix(L*L)).Value(); // (!=Glockenspiel)
- X test = test - 720.0 * 720.0; Clean(test,0.000000001); Print(test);
- X }
- X
- X {
- X Tracer et1("Stage 4");
- X M = S * S;
- X Matrix SX = S;
- X RowVector test(3);
- X test(1) = SumSquare(S);
- X test(2) = SumSquare(SX);
- X test(3) = Trace(M);
- X test = test - 3925961; Print(test);
- X }
- X {
- X Tracer et1("Stage 5");
- X SymmetricMatrix SM(10), SN(10);
- X for (i=1; i<=10; i++) for (j=i; j<=10; j++)
- X {
- X SM(i,j) = 1.5 * i - j; SN(i,j) = SM(i,j) * SM(i,j);
- X if (SM(i,j) < 0.0) SN(i,j) = - SN(i,j);
- X }
- X Matrix M = SM, N = SN; RowVector test(4);
- X test(1) = SumAbsoluteValue(SN);
- X test(2) = SumAbsoluteValue(-SN);
- X test(3) = SumAbsoluteValue(N);
- X test(4) = SumAbsoluteValue(-N);
- X test = test - 1168.75; Print(test);
- X test(1) = MaximumAbsoluteValue(SM);
- X test(2) = MaximumAbsoluteValue(-SM);
- X test(3) = MaximumAbsoluteValue(M);
- X test(4) = MaximumAbsoluteValue(-M);
- X test = test - 8.5; Print(test);
- X }
- X {
- X Tracer et1("Stage 6");
- X Matrix M(15,20); Real value = 0.0;
- X for (i=1; i<=15; i++) { for (j=1; j<=20; j++) M(i,j) = 1.5 * i - j; }
- X for (i=1; i<=20; i++)
- X { Real v = SumAbsoluteValue(M.Column(i)); if (value<v) value = v; }
- X RowVector test(3);
- X test(1) = value;
- X test(2) = Norm1(M);
- X test(3) = NormInfinity(M.t());
- X test = test - 165; Print(test);
- X test.ReDimension(5);
- X ColumnVector CV = M.AsColumn(); M = CV.t();
- X test(1) = Norm1(CV.t());
- X test(2) = MaximumAbsoluteValue(M);
- X test(3) = NormInfinity(CV);
- X test(4) = Norm1(M);
- X test(5) = NormInfinity(M.t());
- X test = test - 21.5; Print(test);
- X }
- X
- X// cout << "\nEnd of Sixteenth test\n";
- X}
- END_OF_FILE
- if test 4734 -ne `wc -c <'tmtg.cxx'`; then
- echo shar: \"'tmtg.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmtg.cxx'
- fi
- if test -f 'tmth.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmth.cxx'\"
- else
- echo shar: Extracting \"'tmth.cxx'\" \(6079 characters\)
- sed "s/^X//" >'tmth.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X
- X#include "include.h"
- X
- X#include "newmatap.h"
- X//#include "newmatio.h"
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- Xvoid Clean(Matrix&, Real);
- X
- X
- X
- X
- Xvoid trymath()
- X{
- X// cout << "\nSeventeenth test of Matrix package\n";
- X// cout << "\n";
- X Tracer et("Seventeenth test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X
- X
- X {
- X Tracer et1("Stage 1");
- X int i, j;
- X BandMatrix B(8,3,1);
- X for (i=1; i<=8; i++) for (j=-3; j<=1; j++)
- X { int k = i+j; if (k>0 && k<=8) B(i,k) = i + k/64.0; }
- X
- X DiagonalMatrix I(8); I = 1; BandMatrix B1; B1 = I;
- X UpperTriangularMatrix UT = I; Print(Matrix(B1-UT));
- X Print(Matrix(B * B - B * 2 + I - (B - I) * (B - I)));
- X Matrix A = B; BandMatrix C; C = B;
- X Print(Matrix(B * A - C * 2 + I - (A - I) * (B - I)));
- X
- X C.ReDimension(8,2,2); C = 0.0; C.Inject(B);
- X Matrix X = A.t();
- X B1.ReDimension(8,2,2); B1.Inject(X); Print(Matrix(C.t()-B1));
- X
- X Matrix BI = B.i(); A = A.i()-BI; Clean(A,0.000000001); Print(A);
- X
- X
- X UpperBandMatrix U; U << B; LowerBandMatrix L; L << B;
- X DiagonalMatrix D; D << B;
- X Print( Matrix(L + (U - D - B)) );
- X
- X
- X
- X for (i=1; i<=8; i++) A.Column(i) << B.Column(i);
- X Print(Matrix(A-B));
- X }
- X {
- X Tracer et1("Stage 2");
- X BandMatrix A(7,2,2);
- X int i,j;
- X for (i=1; i<=7; i++) for (j=1; j<=7; j++)
- X {
- X int k=i-j; if (k<0) k = -k;
- X if (k==0) A(i,j)=6;
- X else if (k==1) A(i,j) = -4;
- X else if (k==2) A(i,j) = 1;
- X A(1,1) = A(7,7) = 5;
- X }
- X DiagonalMatrix D(7); D = 0.0; A = A - D;
- X BandLUMatrix B(A); Matrix M = A;
- X ColumnVector V(3);
- X V(1) = LogDeterminant(B).Value();
- X V(2) = LogDeterminant(A).Value();
- X V(3) = LogDeterminant(M).Value();
- X V = V / 64 - 1; Clean(V,0.000000001); Print(V);
- X ColumnVector X(7);
- X#ifdef ATandT
- X Real a[7]; a[0]=1; a[1]=2; a[2]=3; a[3]=4; a[4]=5; a[5]=6; a[6]=7;
- X#else
- X Real a[] = {1,2,3,4,5,6,7};
- X#endif
- X X << a;
- X M = (M.i()*X).t() - (B.i()*X).t() * 2.0 + (A.i()*X).t();
- X Clean(M,0.000000001); Print(M);
- X
- X
- X BandMatrix P(80,2,5); ColumnVector CX(80);
- X for (i=1; i<=80; i++) for (j=1; j<=80; j++)
- X { int d = i-j; if (d<=2 && d>=-5) P(i,j) = i + j/100.0; }
- X for (i=1; i<=80; i++) CX(i) = i*100.0;
- X Matrix MP = P;
- X ColumnVector V1 = P.i() * CX; ColumnVector V2 = MP.i() * CX;
- X V = V1 - V2; Clean(V,0.000000001); Print(V);
- X
- X V1 = P * V1; V2 = MP * V2; V = V1 - V2; Clean(V,0.000000001); Print(V);
- X RowVector XX(1);
- X XX = LogDeterminant(P).Value() / LogDeterminant(MP).Value() - 1.0;
- X Clean(XX,0.000000001); Print(XX);
- X
- X LowerBandMatrix LP(80,5);
- X for (i=1; i<=80; i++) for (j=1; j<=80; j++)
- X { int d = i-j; if (d<=5 && d>=0) LP(i,j) = i + j/100.0; }
- X MP = LP;
- X XX.ReDimension(4);
- X XX(1) = LogDeterminant(LP).Value();
- X XX(2) = LogDeterminant(MP).Value();
- X V1 = LP.i() * CX; V2 = MP.i() * CX;
- X V = V1 - V2; Clean(V,0.000000001); Print(V);
- X
- X UpperBandMatrix UP(80,4);
- X for (i=1; i<=80; i++) for (j=1; j<=80; j++)
- X { int d = i-j; if (d<=0 && d>=-4) UP(i,j) = i + j/100.0; }
- X MP = UP;
- X XX(3) = LogDeterminant(UP).Value();
- X XX(4) = LogDeterminant(MP).Value();
- X V1 = UP.i() * CX; V2 = MP.i() * CX;
- X V = V1 - V2; Clean(V,0.000000001); Print(V);
- X XX = XX / SumAbsoluteValue(XX) - .25; Clean(XX,0.000000001); Print(XX);
- X }
- X {
- X Tracer et1("Stage 3");
- X SymmetricBandMatrix SA(8,5);
- X int i,j;
- X for (i=1; i<=8; i++) for (j=1; j<=8; j++)
- X if (i-j<=5 && 0<=i-j) SA(i,j) =i + j/128.0;
- X DiagonalMatrix D(8); D = 10; SA = SA + D;
- X
- X Matrix MA1(8,8); Matrix MA2(8,8);
- X for (i=1; i<=8; i++)
- X { MA1.Column(i) << SA.Column(i); MA2.Row(i) << SA.Row(i); }
- X Print(Matrix(MA1-MA2));
- X
- X D = 10; SA = SA.t() + D; MA1 = MA1 + D;
- X Print(Matrix(MA1-SA));
- X
- X UpperBandMatrix UB(8,3); LowerBandMatrix LB(8,4);
- X D << SA; UB << SA; LB << SA;
- X SA = SA * 5.0; D = D * 5.0; LB = LB * 5.0; UB = UB * 5.0;
- X BandMatrix B = LB - D + UB - SA; Print(Matrix(B));
- X
- X SymmetricBandMatrix A(7,2); A = 100.0;
- X for (i=1; i<=7; i++) for (j=1; j<=7; j++)
- X {
- X int k=i-j;
- X if (k==0) A(i,j)=6;
- X else if (k==1) A(i,j) = -4;
- X else if (k==2) A(i,j) = 1;
- X A(1,1) = A(7,7) = 5;
- X }
- X BandLUMatrix C(A); Matrix M = A;
- X ColumnVector X(7);
- X X(1) = LogDeterminant(C).Value() - 64;
- X X(2) = LogDeterminant(A).Value() - 64;
- X X(3) = LogDeterminant(M).Value() - 64;
- X X(4) = SumSquare(M) - SumSquare(A);
- X X(5) = SumAbsoluteValue(M) - SumAbsoluteValue(A);
- X X(6) = MaximumAbsoluteValue(M) - MaximumAbsoluteValue(A);
- X X(7) = Trace(M) - Trace(A);
- X Clean(X,0.000000001); Print(X);
- X
- X#ifdef ATandT
- X Real a[7]; a[0]=1; a[1]=2; a[2]=3; a[3]=4; a[4]=5; a[5]=6; a[6]=7;
- X#else
- X Real a[] = {1,2,3,4,5,6,7};
- X#endif
- X X << a;
- X X = M.i()*X - C.i()*X * 2 + A.i()*X;
- X Clean(X,0.000000001); Print(X);
- X
- X
- X LB << A; UB << A; D << A;
- X BandMatrix XA = LB + (UB - D);
- X Print(Matrix(XA - A));
- X
- X for (i=1; i<=7; i++) for (j=1; j<=7; j++)
- X {
- X int k=i-j;
- X if (k==0) A(i,j)=6;
- X else if (k==1) A(i,j) = -4;
- X else if (k==2) A(i,j) = 1;
- X A(1,1) = A(7,7) = 5;
- X }
- X D = 1;
- X
- X M = LB.i() * LB - D; Clean(M,0.000000001); Print(M);
- X M = UB.i() * UB - D; Clean(M,0.000000001); Print(M);
- X M = XA.i() * XA - D; Clean(M,0.000000001); Print(M);
- X#ifdef GXX
- X Matrix MUB = UB; Matrix MLB = LB;
- X M = LB.i() * UB - LB.i() * MUB; Clean(M,0.000000001); Print(M);
- X M = UB.i() * LB - UB.i() * MLB; Clean(M,0.000000001); Print(M);
- X#else
- X M = LB.i() * UB - LB.i() * Matrix(UB); Clean(M,0.000000001); Print(M);
- X M = UB.i() * LB - UB.i() * Matrix(LB); Clean(M,0.000000001); Print(M);
- X#endif
- X }
- X
- X
- X// cout << "\nEnd of Seventeenth test\n";
- X}
- END_OF_FILE
- if test 6079 -ne `wc -c <'tmth.cxx'`; then
- echo shar: \"'tmth.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmth.cxx'
- fi
- if test -f 'tmti.cxx' -a "${1}" != "-c" ; then
- echo shar: Will not clobber existing file \"'tmti.cxx'\"
- else
- echo shar: Extracting \"'tmti.cxx'\" \(3525 characters\)
- sed "s/^X//" >'tmti.cxx' <<'END_OF_FILE'
- X
- X//#define WANT_STREAM
- X
- X#include "include.h"
- X
- X#include "newmatap.h"
- X//#include "newmatio.h"
- X
- Xvoid Print(const Matrix& X);
- Xvoid Print(const UpperTriangularMatrix& X);
- Xvoid Print(const DiagonalMatrix& X);
- Xvoid Print(const SymmetricMatrix& X);
- Xvoid Print(const LowerTriangularMatrix& X);
- X
- Xvoid Clean(Matrix&, Real);
- X
- Xvoid WillNotConverge()
- X{
- X Matrix A(10,10);
- X Throw(ConvergenceException(A));
- X}
- X
- Xvoid trymati()
- X{
- X Tracer et("Eighteenth test of Matrix package");
- X Exception::PrintTrace(TRUE);
- X ProgramException::SetAction(0); // turn of error messages
- X DataException::SetAction(0);
- X ConvergenceException::SetAction(0);
- X ColumnVector checks(9); checks = 1.0; checks(1) = 0.0;
- X
- X Try { WillNotConverge(); }
- X Catch(ConvergenceException) { checks(2) = 0; }
- X Catch(InternalException) { checks(1) = 1; }
- X Catch(ProgramException) { checks(1) = 1; }
- X Catch(DataException) { checks(1) = 1; }
- X Catch(SpaceException) { checks(1) = 1; }
- X CatchAndThrow;
- X
- X Try { Matrix M(10,10); SymmetricMatrix S = M; }
- X Catch(ConvergenceException) { checks(1) = 1; }
- X Catch(InternalException) { checks(1) = 1; }
- X Catch(ProgramException) { checks(3) = 0; }
- X Catch(DataException) { checks(1) = 1; }
- X Catch(SpaceException) { checks(1) = 1; }
- X CatchAndThrow;
- X
- X Try { Matrix M(10,10); M(10,11) = 2.0; }
- X Catch(ConvergenceException) { checks(1) = 1; }
- X Catch(InternalException) { checks(1) = 1; }
- X Catch(ProgramException) { checks(4) = 0; }
- X Catch(DataException) { checks(1) = 1; }
- X Catch(SpaceException) { checks(1) = 1; }
- X CatchAndThrow;
- X
- X Try { Matrix M(10,10); M = 0.0; M = M.i(); }
- X Catch(ConvergenceException) { checks(1) = 1; }
- X Catch(InternalException) { checks(1) = 1; }
- X Catch(ProgramException) { checks(1) = 1; }
- X Catch(DataException) { checks(5) = 0; }
- X Catch(SpaceException) { checks(1) = 1; }
- X CatchAndThrow;
- X
- X Try { ColumnVector A(30), B(50); A = 5; B = 3; FFT(A,B,A,B); }
- X Catch(ConvergenceException) { checks(1) = 1; }
- X Catch(InternalException) { checks(1) = 1; }
- X Catch(ProgramException) { checks(6) = 0; }
- X Catch(DataException) { checks(1) = 1; }
- X Catch(SpaceException) { checks(1) = 1; }
- X CatchAndThrow;
- X
- X Try
- X {
- X ColumnVector A(30); A = 5; Matrix At = A.t();
- X DiagonalMatrix D;
- X SVD(At,D);
- X }
- X Catch(ConvergenceException) { checks(1) = 1; }
- X Catch(InternalException) { checks(1) = 1; }
- X Catch(ProgramException) { checks(6) = 0; }
- X Catch(DataException) { checks(1) = 1; }
- X Catch(SpaceException) { checks(1) = 1; }
- X CatchAndThrow;
- X
- X Try { BandMatrix X(10,3,4); X(1,10) = 4.0; }
- X Catch(ConvergenceException) { checks(1) = 1; }
- X Catch(InternalException) { checks(1) = 1; }
- X Catch(ProgramException) { checks(7) = 0; }
- X Catch(DataException) { checks(1) = 1; }
- X Catch(SpaceException) { checks(1) = 1; }
- X CatchAndThrow;
- X
- X Try
- X {
- X SymmetricMatrix S(10); S = 5;
- X LowerTriangularMatrix L = Cholesky(S);
- X }
- X Catch(ConvergenceException) { checks(1) = 1; }
- X Catch(InternalException) { checks(1) = 1; }
- X Catch(ProgramException) { checks(1) = 1; }
- X Catch(DataException) { checks(8) = 0; }
- X Catch(SpaceException) { checks(1) = 1; }
- X CatchAndThrow;
- X
- X Try { BandMatrix M(10,3,5); M = 0.0; Matrix XM = M.i(); }
- X Catch(ConvergenceException) { checks(1) = 1; }
- X Catch(InternalException) { checks(1) = 1; }
- X Catch(ProgramException) { checks(1) = 1; }
- X Catch(DataException) { checks(9) = 0; }
- X Catch(SpaceException) { checks(1) = 1; }
- X CatchAndThrow;
- X
- X Print(checks);
- X
- X}
- END_OF_FILE
- if test 3525 -ne `wc -c <'tmti.cxx'`; then
- echo shar: \"'tmti.cxx'\" unpacked with wrong size!
- fi
- # end of 'tmti.cxx'
- fi
- echo shar: End of archive 7 \(of 7\).
- cp /dev/null ark7isdone
- MISSING=""
- for I in 1 2 3 4 5 6 7 ; do
- if test ! -f ark${I}isdone ; then
- MISSING="${MISSING} ${I}"
- fi
- done
- if test "${MISSING}" = "" ; then
- echo You have unpacked all 7 archives.
- rm -f ark[1-9]isdone
- else
- echo You still need to unpack the following archives:
- echo " " ${MISSING}
- fi
- ## End of shell archive.
- exit 0
-
- exit 0 # Just in case...
-