MOAB: Mesh Oriented datABase  (version 5.2.1)
MsqMatrixTest.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2006 Sandia National Laboratories.  Developed at the
00005     University of Wisconsin--Madison under SNL contract number
00006     624796.  The U.S. Government and the University of Wisconsin
00007     retian certain rights to this software.
00008 
00009     This library is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU Lesser General Public
00011     License as published by the Free Software Foundation; either
00012     version 2.1 of the License, or (at your option) any later version.
00013 
00014     This library is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017     Lesser General Public License for more details.
00018 
00019     You should have received a copy of the GNU Lesser General Public License
00020     (lgpl.txt) along with this library; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 
00023     (2006) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file MsqMatrixTest.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "MsqMatrix.hpp"
00034 #include <sstream>
00035 #include "cppunit/extensions/HelperMacros.h"
00036 #include "UnitUtil.hpp"
00037 
00038 using namespace MBMesquite;
00039 using namespace std;
00040 
00041 class MsqMatrixTest : public CppUnit::TestFixture
00042 {
00043   private:
00044     CPPUNIT_TEST_SUITE( MsqMatrixTest );
00045 
00046     CPPUNIT_TEST( test_initialize );
00047     CPPUNIT_TEST( test_assign );
00048     CPPUNIT_TEST( test_matrix_multiply );
00049     CPPUNIT_TEST( test_scalar_multiply );
00050     CPPUNIT_TEST( test_matrix_add );
00051     CPPUNIT_TEST( test_scalar_add );
00052     CPPUNIT_TEST( test_matrix_subtract );
00053     CPPUNIT_TEST( test_scalar_subtract );
00054     CPPUNIT_TEST( test_get_set_row );
00055     CPPUNIT_TEST( test_get_set_column );
00056     CPPUNIT_TEST( test_inverse );
00057     CPPUNIT_TEST( test_qr );
00058     CPPUNIT_TEST( test_frobenius );
00059     CPPUNIT_TEST( test_vec_length );
00060     CPPUNIT_TEST( test_determinant );
00061     CPPUNIT_TEST( test_vec_outer_product );
00062 
00063     CPPUNIT_TEST_SUITE_END();
00064 
00065   public:
00066     void test_initialize();
00067     void test_assign();
00068     void test_matrix_multiply();
00069     void test_scalar_multiply();
00070     void test_matrix_add();
00071     void test_scalar_add();
00072     void test_matrix_subtract();
00073     void test_scalar_subtract();
00074     void test_get_set_row();
00075     void test_get_set_column();
00076     void test_inverse();
00077     void test_qr();
00078     void test_frobenius();
00079     void test_vec_length();
00080     void test_determinant();
00081     void test_vec_outer_product();
00082 };
00083 
00084 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( MsqMatrixTest, "MsqMatrixTest" );
00085 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( MsqMatrixTest, "Unit" );
00086 
00087 void MsqMatrixTest::test_initialize()
00088 {
00089     MsqMatrix< 4, 4 > I44( 1.0 );
00090     ASSERT_IDENTITY_MATRIX( I44 );
00091 
00092     double values[] = { 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 };
00093     MsqMatrix< 3, 3 > I33( values );
00094     ASSERT_IDENTITY_MATRIX( I33 );
00095 
00096     double c0v[] = { 1.0, 0.0 }, c1v[] = { 0.0, 1.0 };
00097     MsqMatrix< 2, 1 > cols[] = { c0v, c1v };
00098     MsqMatrix< 2, 2 > I22( cols );
00099     ASSERT_IDENTITY_MATRIX( I22 );
00100 
00101     double r0v[] = { 1.0, 2.0 }, r1v[] = { 3.0, 4.0 };
00102     MsqMatrix< 1, 2 > rows[] = { r0v, r1v };
00103     MsqMatrix< 2, 2 > I( rows );
00104     CPPUNIT_ASSERT_DOUBLES_EQUAL( r0v[0], I( 0, 0 ), 1e-6 );
00105     CPPUNIT_ASSERT_DOUBLES_EQUAL( r0v[1], I( 0, 1 ), 1e-6 );
00106     CPPUNIT_ASSERT_DOUBLES_EQUAL( r1v[0], I( 1, 0 ), 1e-6 );
00107     CPPUNIT_ASSERT_DOUBLES_EQUAL( r1v[1], I( 1, 1 ), 1e-6 );
00108 
00109     std::string str = "-1 -2 1 2 3 4";
00110     MsqMatrix< 2, 3 > S1( str );
00111     CPPUNIT_ASSERT_DOUBLES_EQUAL( -1, S1( 0, 0 ), 1e-6 );
00112     CPPUNIT_ASSERT_DOUBLES_EQUAL( -2, S1( 0, 1 ), 1e-6 );
00113     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1, S1( 0, 2 ), 1e-6 );
00114     CPPUNIT_ASSERT_DOUBLES_EQUAL( 2, S1( 1, 0 ), 1e-6 );
00115     CPPUNIT_ASSERT_DOUBLES_EQUAL( 3, S1( 1, 1 ), 1e-6 );
00116     CPPUNIT_ASSERT_DOUBLES_EQUAL( 4, S1( 1, 2 ), 1e-6 );
00117 
00118     MsqMatrix< 2, 3 > S2( str.c_str() );
00119     CPPUNIT_ASSERT_DOUBLES_EQUAL( -1, S2( 0, 0 ), 1e-6 );
00120     CPPUNIT_ASSERT_DOUBLES_EQUAL( -2, S2( 0, 1 ), 1e-6 );
00121     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1, S2( 0, 2 ), 1e-6 );
00122     CPPUNIT_ASSERT_DOUBLES_EQUAL( 2, S2( 1, 0 ), 1e-6 );
00123     CPPUNIT_ASSERT_DOUBLES_EQUAL( 3, S2( 1, 1 ), 1e-6 );
00124     CPPUNIT_ASSERT_DOUBLES_EQUAL( 4, S2( 1, 2 ), 1e-6 );
00125 
00126     MsqMatrix< 2, 2 > II( I33, 1, 1 );
00127     ASSERT_IDENTITY_MATRIX( II );
00128 }
00129 
00130 void MsqMatrixTest::test_assign()
00131 {
00132     MsqMatrix< 3, 2 > m;
00133     for( int i = 0; i < 6; ++i )
00134     {
00135         m( i / 2, i % 2 ) = (double)i;
00136         CPPUNIT_ASSERT_DOUBLES_EQUAL( (double)i, m( i / 2, i % 2 ), DBL_EPSILON );
00137     }
00138 }
00139 
00140 void MsqMatrixTest::test_matrix_multiply()
00141 {
00142     double v1[] = { 1, 2, 3, -3, -2, -1 };
00143     double v2[] = { 4, 3, 5, 10, 2, 0.5 };
00144 
00145     MsqMatrix< 2, 3 > A( v1 );
00146     MsqMatrix< 3, 2 > B( v2 );
00147     MsqMatrix< 2, 2 > R = A * B;
00148 
00149     CPPUNIT_ASSERT_DOUBLES_EQUAL( 20.0, R( 0, 0 ), 1e-6 );
00150     CPPUNIT_ASSERT_DOUBLES_EQUAL( 24.5, R( 0, 1 ), 1e-6 );
00151     CPPUNIT_ASSERT_DOUBLES_EQUAL( -24.0, R( 1, 0 ), 1e-6 );
00152     CPPUNIT_ASSERT_DOUBLES_EQUAL( -29.5, R( 1, 1 ), 1e-6 );
00153 }
00154 
00155 void MsqMatrixTest::test_scalar_multiply()
00156 {
00157     MsqMatrix< 4, 4 > A( 33.0 );
00158     MsqMatrix< 4, 4 > B( 1.0 / 33 * A );
00159     A *= 1.0 / 33.0;
00160     ASSERT_IDENTITY_MATRIX( A );
00161     ASSERT_IDENTITY_MATRIX( B );
00162 }
00163 
00164 void MsqMatrixTest::test_matrix_add()
00165 {
00166     double v1[] = { 3, 4, 5, 2, 1, 0, 7, 8, 9 };
00167     double v2[] = { -2, -4, -5, -2, 0, 0, -7, -8, -8 };
00168     MsqMatrix< 3, 3 > A( v1 ), B( v2 );
00169     MsqMatrix< 3, 3 > R = A + B;
00170     ASSERT_IDENTITY_MATRIX( R );
00171     A += B;
00172     ASSERT_IDENTITY_MATRIX( A );
00173 }
00174 
00175 void MsqMatrixTest::test_scalar_add()
00176 {
00177     double v[] = { -6, -7, -7, -6 };
00178     MsqMatrix< 2, 2 > A( v ), B( A + 7 );
00179     A += 7;
00180     ASSERT_IDENTITY_MATRIX( A );
00181     ASSERT_IDENTITY_MATRIX( B );
00182 }
00183 
00184 void MsqMatrixTest::test_matrix_subtract()
00185 {
00186     double v1[] = { 3, 4, 5, 2, 1, 0, 7, 8, 9 };
00187     double v2[] = { 2, 4, 5, 2, 0, 0, 7, 8, 8 };
00188     MsqMatrix< 3, 3 > A( v1 ), B( v2 );
00189     MsqMatrix< 3, 3 > R = A - B;
00190     ASSERT_IDENTITY_MATRIX( R );
00191     A -= B;
00192     ASSERT_IDENTITY_MATRIX( A );
00193 }
00194 
00195 void MsqMatrixTest::test_scalar_subtract()
00196 {
00197     double v[] = { -6, -7, -7, -6 };
00198     MsqMatrix< 2, 2 > A( v ), B( A - -7.0 );
00199     A -= -7.0;
00200     ASSERT_IDENTITY_MATRIX( A );
00201     ASSERT_IDENTITY_MATRIX( B );
00202 }
00203 
00204 void MsqMatrixTest::test_get_set_row()
00205 {
00206     double values[] = { -1, 1, -2, 2, -3, 3, -4, 4, -5 };
00207     MsqMatrix< 3, 3 > A( values );
00208     MsqMatrix< 1, 3 > r0 = A.row( 0 );
00209     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[0], r0( 0, 0 ), DBL_EPSILON );
00210     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[1], r0( 0, 1 ), DBL_EPSILON );
00211     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[2], r0( 0, 2 ), DBL_EPSILON );
00212     MsqMatrix< 1, 3 > r1 = A.row( 1 );
00213     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[3], r1( 0, 0 ), DBL_EPSILON );
00214     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[4], r1( 0, 1 ), DBL_EPSILON );
00215     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[5], r1( 0, 2 ), DBL_EPSILON );
00216     MsqMatrix< 1, 3 > r2 = A.row( 2 );
00217     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[6], r2( 0, 0 ), DBL_EPSILON );
00218     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[7], r2( 0, 1 ), DBL_EPSILON );
00219     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[8], r2( 0, 2 ), DBL_EPSILON );
00220 
00221     r0( 0, 0 ) = 20;
00222     r0( 0, 1 ) = 30;
00223     r0( 0, 2 ) = 10;
00224     A.set_row( 0, r0 );
00225     CPPUNIT_ASSERT_DOUBLES_EQUAL( A( 0, 0 ), 20, DBL_EPSILON );
00226     CPPUNIT_ASSERT_DOUBLES_EQUAL( A( 0, 1 ), 30, DBL_EPSILON );
00227     CPPUNIT_ASSERT_DOUBLES_EQUAL( A( 0, 2 ), 10, DBL_EPSILON );
00228 }
00229 
00230 void MsqMatrixTest::test_get_set_column()
00231 {
00232     double values[] = { -1, 1, -2, 2, -3, 3, -4, 4, -5 };
00233     MsqMatrix< 3, 3 > A( values );
00234     MsqMatrix< 3, 1 > c0 = A.column( 0 );
00235     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[0], c0( 0, 0 ), DBL_EPSILON );
00236     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[3], c0( 1, 0 ), DBL_EPSILON );
00237     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[6], c0( 2, 0 ), DBL_EPSILON );
00238     MsqMatrix< 3, 1 > c1 = A.column( 1 );
00239     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[1], c1( 0, 0 ), DBL_EPSILON );
00240     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[4], c1( 1, 0 ), DBL_EPSILON );
00241     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[7], c1( 2, 0 ), DBL_EPSILON );
00242     MsqMatrix< 3, 1 > c2 = A.column( 2 );
00243     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[2], c2( 0, 0 ), DBL_EPSILON );
00244     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[5], c2( 1, 0 ), DBL_EPSILON );
00245     CPPUNIT_ASSERT_DOUBLES_EQUAL( values[8], c2( 2, 0 ), DBL_EPSILON );
00246 
00247     c0( 0, 0 ) = 20;
00248     c0( 1, 0 ) = 30;
00249     c0( 2, 0 ) = 10;
00250     A.set_column( 0, c0 );
00251     CPPUNIT_ASSERT_DOUBLES_EQUAL( A( 0, 0 ), 20, DBL_EPSILON );
00252     CPPUNIT_ASSERT_DOUBLES_EQUAL( A( 1, 0 ), 30, DBL_EPSILON );
00253     CPPUNIT_ASSERT_DOUBLES_EQUAL( A( 2, 0 ), 10, DBL_EPSILON );
00254 }
00255 
00256 void MsqMatrixTest::test_inverse()
00257 {
00258     double v2[] = { 4, 5, 6, 7 };
00259     MsqMatrix< 2, 2 > M2( v2 );
00260     MsqMatrix< 2, 2 > M2I = inverse( M2 );
00261     MsqMatrix< 2, 2 > M2P = M2 * M2I;
00262     ASSERT_IDENTITY_MATRIX( M2P );
00263 
00264     double v3[] = { 0.1, 0.5, 0.3, 1, 2, 3, 0.5, 0.1, 10 };
00265     MsqMatrix< 3, 3 > M3( v3 );
00266     MsqMatrix< 3, 3 > M3I = inverse( M3 );
00267     MsqMatrix< 3, 3 > M3P = M3 * M3I;
00268     ASSERT_IDENTITY_MATRIX( M3P );
00269 }
00270 
00271 void MsqMatrixTest::test_qr()
00272 {
00273     double v2[] = { 4, 5, 6, 7 };
00274     MsqMatrix< 2, 2 > M2( v2 ), Q2, R2;
00275     QR( M2, Q2, R2 );
00276     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, R2( 1, 0 ), 1e-6 );
00277     MsqMatrix< 2, 2 > P2 = Q2 * R2;
00278     ASSERT_MATRICES_EQUAL( M2, P2, 1e-6 );
00279     MsqMatrix< 2, 2 > I2 = Q2 * transpose( Q2 );
00280     ASSERT_IDENTITY_MATRIX( I2 );
00281 
00282     double v3[] = { 0.1, 0.5, 0.3, 1, 2, 3, 0.5, 0.1, 10 };
00283     MsqMatrix< 3, 3 > M3( v3 ), Q3, R3;
00284     QR( M3, Q3, R3 );
00285     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, R3( 1, 0 ), 1e-6 );
00286     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, R3( 2, 0 ), 1e-6 );
00287     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, R3( 2, 1 ), 1e-6 );
00288     MsqMatrix< 3, 3 > P3 = Q3 * R3;
00289     ASSERT_MATRICES_EQUAL( M3, P3, 1e-6 );
00290     MsqMatrix< 3, 3 > I3 = Q3 * transpose( Q3 );
00291     ASSERT_IDENTITY_MATRIX( I3 );
00292 }
00293 
00294 void MsqMatrixTest::test_frobenius()
00295 {
00296     double v3[] = { 0.1, 0.5, 0.3, 1, 2, 3, 0.5, 0.1, 10 };
00297     MsqMatrix< 3, 3 > M3( v3 );
00298     double exp = 0.0;
00299     for( unsigned i = 0; i < 9; ++i )
00300         exp += v3[i] * v3[i];
00301     CPPUNIT_ASSERT_DOUBLES_EQUAL( exp, sqr_Frobenius( M3 ), 1e-10 );
00302 }
00303 
00304 void MsqMatrixTest::test_vec_length()
00305 {
00306     double values[] = { 2, 4, 4 };
00307     MsqMatrix< 3, 1 > c( values );
00308     MsqMatrix< 1, 3 > r( values );
00309     CPPUNIT_ASSERT_DOUBLES_EQUAL( 36.0, sqr_length( c ), 1e-12 );
00310     CPPUNIT_ASSERT_DOUBLES_EQUAL( 36.0, sqr_length( r ), 1e-12 );
00311     CPPUNIT_ASSERT_DOUBLES_EQUAL( 6.0, length( c ), 1e-12 );
00312     CPPUNIT_ASSERT_DOUBLES_EQUAL( 6.0, length( r ), 1e-12 );
00313 }
00314 
00315 void MsqMatrixTest::test_determinant()
00316 {
00317     MsqMatrix< 5, 5 > m5( 1.0 );
00318     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( m5 ), 1e-8 );
00319 
00320     double m4vals[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 };
00321     MsqMatrix< 4, 4 > m4( m4vals );
00322     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, det( m4 ), 1e-8 );
00323 
00324     double m4vals2[] = { 2, 5, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1 };
00325     MsqMatrix< 4, 4 > m4b( m4vals2 );
00326     CPPUNIT_ASSERT_DOUBLES_EQUAL( -3.0, det( m4b ), 1e-8 );
00327 
00328     double m3vals[] = { 1, 2, -1, -1, 0, 1, 1, 1, 3 };
00329     MsqMatrix< 3, 3 > m3( m3vals );
00330     CPPUNIT_ASSERT_DOUBLES_EQUAL( 8.0, det( m3 ), 1e-8 );
00331 
00332     double m2vals[] = { sqrt( 2.0 ), 1, 1, sqrt( 2.0 ) };
00333     MsqMatrix< 2, 2 > m2( m2vals );
00334     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( m2 ), 1e-8 );
00335 }
00336 
00337 void MsqMatrixTest::test_vec_outer_product()
00338 {
00339     double v1[] = { 1, 3, -1, 12 };
00340     double v2[] = { -10, 4, 0.5, -9, 2 };
00341     MsqMatrix< 4, 1 > V1( v1 );
00342     MsqMatrix< 5, 1 > V2( v2 );
00343 
00344     MsqMatrix< 4, 5 > M1 = V1 * transpose( V2 );
00345     MsqMatrix< 4, 5 > M2 = outer( V1, V2 );
00346     ASSERT_MATRICES_EQUAL( M1, M2, 1e-8 );
00347 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines