//----------------------------------------------------------- // TOOLS // C++ Templates & Tools, 2nd Edition //----------------------------------------------------------- // // TestMatrix.cpp // // Procedures for testing matrix classes // //----------------------------------------------------------- // Copyright 1996 by Scott Robert Ladd. All rights reserved. //----------------------------------------------------------- #include "strstream.h" #include "iomanip.h" #include "limits.h" #include "matrix.h" //using namespace Coyote; // utility routine to display a matrix template void ShowMatrix ( const Matrix & m ) { const size_t MaxR = m.GetRows(); const size_t MaxC = m.GetCols(); size_t r, c; for (r = 0; r < MaxR; ++r) { cout << "[ "; for (c = 0; c < MaxC; ++c) cout << setw(3) << m.Get(r,c) << " "; cout << "]\r\n"; } } int Times2(const int & x) { return x * 2; } //-------------- // Test matrices //-------------- void main ( void * args ) { // display a headline cout << "Matrix test\r\n===========\r\n"; Matrix A(3,3), B(3,3), C(3,3), D(3,3); A(0,0) = 1; A(0,1) = 3; A(0,2) = -4; A(1,0) = 1; A(1,1) = 1; A(1,2) = -2; A(2,0) = -1; A(2,1) = -2; A(2,2) = 5; B(0,0) = 8; B(0,1) = 3; B(0,2) = 0; B(1,0) = 3; B(1,1) = 10; B(1,2) = 2; B(2,0) = 0; B(2,1) = 2; B(2,2) = 6; D(0,0) = 1; D(0,1) = 2; D(0,2) = -1; D(1,0) = 2; D(1,1) = -1; D(1,2) = -3; D(2,0) = 0; D(2,1) = -2; D(2,2) = 4; cout << "\r\nMatrix A = \r\n"; ShowMatrix(A); cout << "\r\nMatrix B = \r\n"; ShowMatrix(B); C = A % B; cout << "\r\nMatrix C (A % B) = \r\n"; ShowMatrix(C); C = A + B; cout << "\r\nMatrix C (A + B) = \r\n"; ShowMatrix(C); C = A; C += B; cout << "\r\nMatrix C (= A, += B) =\r\n"; ShowMatrix(C); C = A + 1; cout << "\r\nMatrix C (= A + 1) =\r\n"; ShowMatrix(C); C += 1; cout << "\r\nMatrix C (+= 1) =\r\n"; ShowMatrix(C); C = A - B; cout << "\r\nMatrix C (A - B) = \r\n"; ShowMatrix(C); C = A; C -= B; cout << "\r\nMatrix C (= A, -= B) =\r\n"; ShowMatrix(C); C = A - 1; cout << "\r\nMatrix C (= A - 1) =\r\n"; ShowMatrix(C); C -= 1; cout << "\r\nMatrix C (-= 1) =\r\n"; ShowMatrix(C); C = A * B; cout << "\r\nMatrix C (A * B) = \r\n"; ShowMatrix(C); C = A; C *= B; cout << "\r\nMatrix C (= A, *= B) =\r\n"; ShowMatrix(C); C = A * 2; cout << "\r\nMatrix C (= A * 2) =\r\n"; ShowMatrix(C); C *= 2; cout << "\r\nMatrix C (*= 2) =\r\n"; ShowMatrix(C); C = B / A; cout << "\r\nMatrix C (B / A) = \r\n"; ShowMatrix(C); C = B; C /= A; cout << "\r\nMatrix C (= B, /= A) =\r\n"; ShowMatrix(C); C = A / 2; cout << "\r\nMatrix C (= A / 2) =\r\n"; ShowMatrix(C); C /= 2; cout << "\r\nMatrix C (/= 2) =\r\n"; ShowMatrix(C); C = -A; cout << "\r\nMatrix C (-A) = \r\n"; ShowMatrix(C); // test comparisons cout << "\r\nMatrix A = \r\n"; ShowMatrix(A); cout << "\r\nMatrix D = \r\n"; ShowMatrix(D); if (A.Equals(D)) cout << "\r\nERROR: A should not equal D"; else cout << "\r\nOKAY: A not equal D"; C = A; if (A.Equals(C)) cout << "\r\nOKAY: A equals C\r\n"; else cout << "\r\nERROR: A should equal C\r\n"; Matrix I(3,3); I = (A == D); cout << "\r\nMatrix I = (A == D)\r\n"; ShowMatrix(I); I = (A != D); cout << "\r\nMatrix I = (A != D)\r\n"; ShowMatrix(I); I = (A < D); cout << "\r\nMatrix I = (A < D)\r\n"; ShowMatrix(I); I = (A <= D); cout << "\r\nMatrix I = (A <= D)\r\n"; ShowMatrix(I); I = (A > D); cout << "\r\nMatrix I = (A > D)\r\n"; ShowMatrix(I); I = (A >= D); cout << "\r\nMatrix I = (A >= D)\r\n"; ShowMatrix(I); // check fill function C.Fill(9); cout << "\r\nC filled with 9 =\r\n"; ShowMatrix(C); // check Apply functions C = Apply(A, Times2); cout << "\r\nC = A.Apply(Times2)\r\n"; ShowMatrix(C); C.Apply(Times2); cout << "\r\nApply(C,Times2)\r\n"; ShowMatrix(C); // check row and column vector functions Matrix S(1,1); Matrix r1A(3,1); Matrix c0B(1,3); r1A = A.VectorRow(1); c0B = B.VectorCol(0); cout << "\r\nMatrix S = \r\n"; ShowMatrix(S); cout << "\r\nMatrix R1A = \r\n"; ShowMatrix(r1A); cout << "\r\nMatrix C0B = \r\n"; ShowMatrix(c0B); if (r1A.IsRowVector()) cout << "\r\nOKAY: R1A is row vector"; else cout << "\r\nERROR: R1A should be a row vector"; if (!r1A.IsColVector()) cout << "\r\nOKAY: R1A is not a column vector"; else cout << "\r\nERROR: R1A should not be a column vector"; if (!c0B.IsRowVector()) cout << "\r\nOKAY: C0B is not a row vector"; else cout << "\r\nERROR: C0B should not be a row vector"; if (c0B.IsColVector()) cout << "\r\nOKAY: C0B is column vector"; else cout << "\r\nERROR: C0B should be a column vector"; if (c0B.IsVector()) cout << "\r\nOKAY: C0B is a vector"; else cout << "\r\nERROR: C0B should be a vector"; if (!A.IsVector()) cout << "\r\nOKAY: A is not a vector"; else cout << "\r\nERROR: A should not be a vector"; if (!c0B.IsSquare()) cout << "\r\nOKAY: C0B is not square"; else cout << "\r\nERROR: C0B should not be square"; if (A.IsSquare()) cout << "\r\nOKAY: A is square"; else cout << "\r\nERROR: A should be square"; B.Fill(0.0); if (B.IsZero()) cout << "\r\nOKAY: B is zero"; else cout << "\r\nERROR: B should be zero"; if (!A.IsZero()) cout << "\r\nOKAY: A is not zero"; else cout << "\r\nERROR: A should not be zero"; // test inner product int ip = r1A.InnerProduct(c0B); cout << "\r\n\r\ninner product of R1A and C0B = " << ip << "\r\n"; // make some bigger matrices Matrix M1(5,5), M2(5,5,3), M3(5,5), M4(5,5); const int junk[] = { 1, 5, 3, 0, 1, 0, 2, 0, 4, 5, 1, 0, 0, 2, 3, 7, 1, 3, 0, 0, 2, 1, 0, 4, 6 }; const int ident[] = { 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1 }; const int tridi[] = { 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1 }; const int utri[] = { 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1 }; const int ltri[] = { 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1 }; const int perm[] = { 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0 }; const int det[] = { 3, 5, 3, 8, 1, 2, 6, 3, 4, 5, 1, 4, 5, 2, 3, 7, 1, 3, 6, 8, 2, 4, 1, 4, 9 }; M1 = ident; M3 = M1 * 2; M4 = junk; cout << "\r\nmatrix M1 = \r\n"; ShowMatrix(M1); cout << "\r\nmatrix M2 = \r\n"; ShowMatrix(M2); cout << "\r\nmatrix M3 = \r\n"; ShowMatrix(M3); cout << "\r\nmatrix M4 = \r\n"; ShowMatrix(M4); if (M1.IsDiagonal()) cout << "\r\nOKAY: M1 is diagonal"; else cout << "\r\nERROR: M1 should be diagonal"; if (M1.IsIdentity()) cout << "\r\nOKAY: M1 is an identity matrix"; else cout << "\r\nERROR: M1 should be an identity matrix"; if (!M2.IsDiagonal()) cout << "\r\nOKAY: M2 is not diagonal"; else cout << "\r\nERROR: M2 should not be diagonal"; if (!M2.IsIdentity()) cout << "\r\nOKAY: M2 is not an identity matrix"; else cout << "\r\nERROR: M2 should not be an identity matrix"; if (M3.IsDiagonal()) cout << "\r\nOKAY: M3 is diagonal"; else cout << "\r\nERROR: M3 should be diagonal"; if (!M3.IsIdentity()) cout << "\r\nOKAY: M3 is not an identity matrix"; else cout << "\r\nERROR: M3 should not be an identity matrix"; if (!M4.IsDiagonal()) cout << "\r\nOKAY: M4 is not diagonal"; else cout << "\r\nERROR: M4 should not be diagonal"; if (!M4.IsIdentity()) cout << "\r\nOKAY: M4 is not an identity matrix"; else cout << "\r\nERROR: M4 should not be an identity matrix"; // tridiagonal tests M1 = tridi; cout << "\r\n\r\nmatrix M1 = \r\n"; ShowMatrix(M1); if (M1.IsTridiagonal()) cout << "\r\nOKAY: M1 is tridiagonal"; else cout << "\r\nERROR: M1 should be tridiagonal"; if (!M4.IsTridiagonal()) cout << "\r\nOKAY: M4 is not tridiagonal"; else cout << "\r\nERROR: M1 should not be tridiagonal"; // upper triangular tests M1 = utri; cout << "\r\n\r\nmatrix M1 = \r\n"; ShowMatrix(M1); if (M1.IsUpperTriangular()) cout << "\r\nOKAY: M1 is upper-triangular"; else cout << "\r\nERROR: M1 should be upper-triangular"; if (!M4.IsUpperTriangular()) cout << "\r\nOKAY: M4 is not upper-triangular"; else cout << "\r\nERROR: M4 should not be upper-triangular"; // lower triangular tests M1 = ltri; cout << "\r\n\r\nmatrix M1 = \r\n"; ShowMatrix(M1); if (M1.IsLowerTriangular()) cout << "\r\nOKAY: M1 is lower-triangular"; else cout << "\r\nERROR: M1 should be lower-triangular"; if (!M4.IsLowerTriangular()) cout << "\r\nOKAY: M4 is not lower-triangular"; else cout << "\r\nERROR: M4 should not be lower-triangular"; // permutation tests M1 = perm; cout << "\r\n\r\nmatrix M1 = \r\n"; ShowMatrix(M1); M2 = ident; cout << "\r\n\r\nmatrix M2 = \r\n"; ShowMatrix(M2); if (M1.IsPermutation()) cout << "\r\nOKAY: M1 is permutation matrix"; else cout << "\r\nERROR: M1 should be permutation"; if (M2.IsPermutation()) cout << "\r\nOKAY: M2 is permutation matrix"; else cout << "\r\nERROR: M2 should be permutation"; if (!M4.IsPermutation()) cout << "\r\nOKAY: M4 is not permutation"; else cout << "\r\nERROR: M4 should not be permutation"; // check singularity function M1(0,1) = 0; cout << "\r\n\r\nmatrix M1 = \r\n"; ShowMatrix(M1); if (M1.IsSingular()) cout << "\r\nOKAY: M1 is singular"; else cout << "\r\nERROR: M1 should be singular"; if (!M2.IsSingular()) cout << "\r\nOKAY: M2 is not singular"; else cout << "\r\nERROR: M2 should not be singular"; if (!M4.IsSingular()) cout << "\r\nOKAY: M4 is not singular"; else cout << "\r\nERROR: M4 should not be singular"; // test minors and determinants cout << "\r\n\r\nmatrix M4 = \r\n"; ShowMatrix(M4); cout << "\r\nminor M4(1,1) = \r\n"; ShowMatrix(M4.Minor(1,1)); cout << "\r\nminor M4(0,4) = \r\n"; ShowMatrix(M4.Minor(0,4)); Matrix M5(2,2), M6(3,3); M5(0,0) = 1; M5(0,1) = 2; M5(1,0) = 3; M5(1,1) = 4; M6(0,0) = 1; M6(0,1) = 3; M6(0,2) = 2; M6(1,0) = 5; M6(1,1) = 4; M6(1,2) = 7; M6(2,0) = 6; M6(2,1) = 9; M6(2,2) = 8; M4 = det; Matrix T4(5,5), T5(2,2), T6(3,3); T4 = M4.Transpose(); T5 = M5.Transpose(); T6 = M6.Transpose(); cout << "\r\nmatrix M5 = \r\n"; ShowMatrix(M5); cout << "\r\ndeterminant of M5 = " << M5.Determinant() << "\r\n"; cout << "\r\nmatrix T5 = \r\n"; ShowMatrix(T5); cout << "\r\ndeterminant of T5 = " << T5.Determinant() << "\r\n"; cout << "\r\nmatrix M6 = \r\n"; ShowMatrix(M6); cout << "\r\ndeterminant of M6 = " << M6.Determinant() << "\r\n"; cout << "\r\nmatrix T6 = \r\n"; ShowMatrix(T6); cout << "\r\ndeterminant of T6 = " << T6.Determinant() << "\r\n"; cout << "\r\nmatrix M4 = \r\n"; ShowMatrix(M4); cout << "\r\ndeterminant of M4 = " << M4.Determinant() << "\r\n"; cout << "\r\nmatrix T4 = \r\n"; ShowMatrix(T4); cout << "\r\ndeterminant of T4 = " << T4.Determinant() << "\r\n"; Matrix R; cout << "\r\nMatrix R (def. constr.) = \r\n"; ShowMatrix(R); R.Resize(10,10); cout << "\r\nMatrix R (now 10x10) = \r\n"; ShowMatrix(R); // check Matrix cout << "\r\nFLOATING POINT!"; Matrix X(3,4), Y(4,3), Z(3,3); X(0,0) = 1.0; X(1,0) = 5.0; X(2,0) = 2.0; X(0,1) = 2.0; X(1,1) = 2.0; X(2,1) = 4.0; X(0,2) = 0.0; X(1,2) = 3.0; X(2,2) = 3.0; X(0,3) = 1.0; X(1,3) = 2.0; X(2,3) = 1.0; Y(0,0) = 0.0; Y(2,0) = 1.0; Y(0,1) = 1.0; Y(2,1) = 0.0; Y(0,2) = 2.0; Y(2,2) = 5.0; Y(1,0) = 1.0; Y(3,0) = 3.0; Y(1,1) = 3.0; Y(3,1) = 1.0; Y(1,2) = 2.0; Y(3,2) = 2.0; cout << "\r\nMatrix X = \r\n"; ShowMatrix(X); cout << "\r\nMatrix Y = \r\n"; ShowMatrix(Y); Z = X % Y; cout << "\r\nMatrix Z (X % Y) = \r\n"; ShowMatrix(Z); // check transposition Matrix tX; tX = X.Transpose(); cout << "\r\nOriginal X =\r\n"; ShowMatrix(X); cout << "\r\nTranspose X =\r\n"; ShowMatrix(tX); X(0,0) = 1; X(0,1) = 3; X(0,2) = -4; X(0,3) = 8; X(1,0) = 1; X(1,1) = 1; X(1,2) = -2; X(1,3) = 2; X(2,0) = -1; X(2,1) = -2; X(2,2) = 5; X(2,3) = -1; cout << "\r\nOriginal X =\r\n"; ShowMatrix(X); Matrix lX(X.LinSolve()); cout << "\r\nX after elimination =\r\n"; ShowMatrix(X); cout << "\r\nlinear equation solution =\r\n"; ShowMatrix(lX); X(0,0) = 1.0; X(1,0) = 3.0; X(2,0) = 5.0; X(0,1) = 2.0; X(1,1) = 5.0; X(2,1) = 6.0; X(0,2) = 0.0; X(1,2) = 4.0; X(2,2) = 3.0; X(0,3) = 0.1; X(1,3) = 12.5; X(2,3) = 10.3; cout << "\r\nOriginal X =\r\n"; ShowMatrix(X); lX = X.LinSolve(); cout << "\r\nX after elimination =\r\n"; ShowMatrix(X); cout << "\r\nlinear equation solution =\r\n"; ShowMatrix(lX); Matrix Adbl(3,3), Bdbl(3,1); Adbl(0,0) = 1.0; Adbl(0,1) = 2.0; Adbl(0,2) = 0.0; Adbl(1,0) = 3.0; Adbl(1,1) = 5.0; Adbl(1,2) = 4.0; Adbl(2,0) = 5.0; Adbl(2,1) = 6.0; Adbl(2,2) = 3.0; Bdbl(0,0) = 0.1; Bdbl(1,0) = 12.5; Bdbl(2,0) = 10.3; cout << "\r\n\r\nmatrix Adbl = \r\n"; ShowMatrix(Adbl); cout << "\r\nmatrix Bdbl = \r\n"; ShowMatrix(Bdbl); Matrix alup(Adbl); // copy Adbl before LUP decomp cout << "\r\nLU decomp of Adbl (before) = \r\n"; ShowMatrix(alup); Matrix aperm = alup.LUPDecompose(); cout << "\r\nLU decomp of Adbl (after) = \r\n"; ShowMatrix(alup); cout << "\r\nPermutation of Adbl = \r\n"; ShowMatrix(aperm); Matrix asol = alup.LUPSolve(aperm,Bdbl); cout << "\r\nlinear solution of Adbl and Bdbl = \r\n"; ShowMatrix(asol); Matrix ainv = alup.LUPInvert(aperm); cout << "\r\ninverse of Adbl and Bdbl = \r\n"; ShowMatrix(ainv); Matrix aid = Adbl % ainv; cout << "\r\ninverse dot Adbl = \r\n"; ShowMatrix(aid); Grid iperm = ainv.LUPDecompose(); Matrix invinv = ainv.LUPInvert(iperm); cout << "\r\ninverse of inverse =\r\n"; ShowMatrix(invinv); }