home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / lin_alg.lha / lin_alg / vmatrix.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1993-08-09  |  14.1 KB  |  460 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *            Verify Operations on Matrices
  6.  *
  7.  ************************************************************************
  8.  */
  9.  
  10. #include "LinAlg.h"
  11. #include <math.h>
  12. #include <builtin.h>
  13. #include <ostream.h>
  14.  
  15. /*
  16.  *------------------------------------------------------------------------
  17.  *            Service validation functions
  18.  */
  19.  
  20. static void verify_element_value(const Matrix& m,const REAL val)
  21. {
  22.   register imax = 0, jmax = 0;
  23.   register double max_dev = 0;
  24.   register int i,j;
  25.   for(i=m.q_row_lwb(); i<=m.q_row_upb(); i++)
  26.     for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
  27.     {
  28.       register double dev = abs(m(i,j)-val);
  29.       if( dev >= max_dev )
  30.     imax = i, jmax = j, max_dev = dev;
  31.     }
  32.  
  33.   if( max_dev == 0 )
  34.     return;
  35.   else if( max_dev < 1e-5 )
  36.     message("Element (%d,%d) with value %g differs the most from what\n"
  37.         "was expected, %g, though the deviation %g is small\n",
  38.         imax,jmax,m(imax,jmax),val,max_dev);
  39.   else
  40.     _error("A significant difference from the expected value %g\n"
  41.        "encountered for element (%d,%d) with value %g",
  42.        val,imax,jmax,m(imax,jmax));
  43. }
  44.  
  45. static void verify_matrix_identity(const Matrix& m1, const Matrix& m2)
  46. {
  47.   register imax = 0, jmax = 0;
  48.   register double max_dev = 0;
  49.   register int i,j;
  50.   are_compatible(m1,m2);
  51.   for(i=m1.q_row_lwb(); i<=m1.q_row_upb(); i++)
  52.     for(j=m1.q_col_lwb(); j<=m1.q_col_upb(); j++)
  53.     {
  54.       register double dev = abs(m1(i,j)-m2(i,j));
  55.       if( dev >= max_dev )
  56.     imax = i, jmax = j, max_dev = dev;
  57.     }
  58.  
  59.   if( max_dev == 0 )
  60.     return;
  61.   if( max_dev < 1e-5 )
  62.     message("Two (%d,%d) elements of matrices with values %g and %g\n"
  63.         "differ the most, though the deviation %g is small\n",
  64.         imax,jmax,m1(imax,jmax),m2(imax,jmax),max_dev);
  65.   else
  66.     _error("A significant difference between the matrices encountered\n"
  67.        "at (%d,%d) element, with values %g and %g",
  68.        imax,jmax,m1(imax,jmax),m2(imax,jmax));
  69. }
  70.  
  71. /*
  72.  *------------------------------------------------------------------------
  73.  *       Test allocation functions and compatibility check
  74.  */
  75.  
  76. static void test_allocation(void)
  77. {
  78.  
  79.   cout << "\n\n---> Test allocation and compatibility check\n";
  80.  
  81.   Matrix m1(4,20);    m1.set_name("Matrix m1");
  82.   Matrix m2(1,4,1,20);    m2.set_name("Matrix m2");
  83.   Matrix m3(1,4,0,19);    m3.set_name("Matrix m3");
  84.   Matrix m4(m1);    m4.set_name("Matrix m4");
  85.   cout << "The following matrices have been allocated\n";
  86.   m1.info(), m2.info(), m3.info(), m4.info();
  87.  
  88.   cout << "\nStatus information reported for matrix m3:\n";
  89.   cout << "  Row lower bound ... " << m3.q_row_lwb() << "\n";
  90.   cout << "  Row upper bound ... " << m3.q_row_upb() << "\n";
  91.   cout << "  Col lower bound ... " << m3.q_col_lwb() << "\n";
  92.   cout << "  Col upper bound ... " << m3.q_col_upb() << "\n";
  93.   cout << "  No. rows ..........." << m3.q_nrows() << "\n";
  94.   cout << "  No. cols ..........." << m3.q_ncols() << "\n";
  95.   cout << "  No. of elements ...." << m3.q_no_elems() << "\n";
  96.   cout << "  Name " << m3.q_name() << "\n";
  97.  
  98.   cout << "\nCheck matrices 1 & 2 for compatibility\n";
  99.   are_compatible(m1,m2);
  100.  
  101.   cout << "Check matrices 1 & 4 for compatibility\n";
  102.   are_compatible(m1,m4);
  103.  
  104. #if 0
  105.   cout << "Check matrices 1 & 3 for compatibility\n";
  106.   are_compatible(m1,m3);
  107. #endif
  108.  
  109.   cout << "m2 has to be compatible with m3 after resizing to m3\n";
  110.   m2.resize_to(m3);
  111.   are_compatible(m2,m3);
  112.  
  113.   Matrix m5(m1.q_nrows()+1,m1.q_ncols()+5); m5.set_name("Matrix m5");
  114.   cout << "m1 has to be compatible with m5 after resizing to m5\n";
  115.   m5.info();
  116.   m1.resize_to(m5.q_nrows(),m5.q_ncols());
  117.   are_compatible(m1,m5);
  118.  
  119. }
  120.  
  121. /*
  122.  *------------------------------------------------------------------------
  123.  *             Test uniform element operations
  124.  */
  125.  
  126. static void test_element_op(const int rsize, const int csize)
  127. {
  128.   const double pattern = 8.625;
  129.   register int i,j;
  130.  
  131.   cout << "\n\n---> Test operations that treat each element uniformly\n";
  132.  
  133.   Matrix m(-1,rsize-2,1,csize);
  134.   Matrix m1(m);
  135.  
  136.   cout << "Writing zeros to m...\n";
  137.   for(i=m.q_row_lwb(); i<=m.q_row_upb(); i++)
  138.     for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
  139.       m(i,j) = 0;
  140.   verify_element_value(m,0);
  141.  
  142.   cout << "Clearing m1 ...\n";
  143.   m1.clear();
  144.   verify_element_value(m1,0);
  145.  
  146.   cout << "Comparing m1 with 0 ...\n";
  147.   assure(m1 == 0, "m1 is not zero!");
  148.  
  149.   cout << "Writing a pattern " << pattern << " by assigning to m(i,j)...\n";
  150.   for(i=m.q_row_lwb(); i<=m.q_row_upb(); i++)
  151.     for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
  152.       m(i,j) = pattern;
  153.   verify_element_value(m,pattern);
  154.  
  155.   cout << "Writing the pattern by assigning to m1 as a whole ...\n";
  156.   m1 = pattern;
  157.   verify_element_value(m1,pattern);
  158.  
  159.   cout << "Comparing m and m1 ...\n";
  160.   assure(m == m1, "m and m1 are unexpectedly different!");
  161.   cout << "Comparing (m=0) and m1 ...\n";
  162.   assert(!(m.clear() == m1));
  163.  
  164.   cout << "\nClear m and add the pattern\n";
  165.   m.clear();
  166.   m += pattern;
  167.   verify_element_value(m,pattern);
  168.   cout << "   add the doubled pattern with the negative sign\n";
  169.   m += -2*pattern;
  170.   verify_element_value(m,-pattern);
  171.   cout << "   subtract the trippled pattern with the negative sign\n";
  172.   m -= -3*pattern;
  173.   verify_element_value(m,2*pattern);
  174.  
  175.   cout << "\nVerify comparison operations\n";
  176.   m = pattern;
  177.   cout << "   (m=pattern)>0\n";
  178.   assure( m > 0, "Unexpectedly, v is not greater than 0");
  179.   cout << "   (m=pattern)>-pattern\n";
  180.   assert( m > -pattern );
  181.   cout << "   (m=-pattern)<-pattern/2\n";
  182.   m -= 2*pattern;
  183.   assert( m  < -pattern/2 );
  184.  
  185.   cout << "\nAssign 2*pattern to m by repeating additions\n";
  186.   m = 0; m += pattern; m += pattern;
  187.   cout << "Assign 2*pattern to m1 by multiplying by two \n";
  188.   m1 = pattern; m1 *= 2;
  189.   verify_element_value(m1,2*pattern);
  190.   assert( m == m1 );
  191.   cout << "Multiply m1 by one half returning it to the 1*pattern\n";
  192.   m1 *= 1/2.;
  193.   verify_element_value(m1,pattern);
  194.  
  195.   cout << "\nAssign -pattern to m and m1\n";
  196.   m.clear(); m -= pattern; m1 = -pattern;
  197.   verify_element_value(m,-pattern);
  198.   assert( m == m1 );
  199.   cout << "m = sqrt(sqr(m)); m1 = abs(m1); Now m and m1 have to be the same\n";
  200.   m.sqr();
  201.   verify_element_value(m,pattern*pattern);
  202.   m.sqrt();
  203.   verify_element_value(m,pattern);
  204.   m1.abs();
  205.   verify_element_value(m1,pattern);
  206.   assert( m == m1 );
  207.  
  208.   cout << "\nCheck out to see that sin^2(x) + cos^2(x) = 1\n";
  209.   for(i=m.q_row_lwb(); i<=m.q_row_upb(); i++)
  210.     for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
  211.       m(i,j) = 4*PI/m.q_no_elems() * (i*m.q_ncols()+j);
  212.   m1 = m;
  213.   m.user_function(sin);
  214.   m1.user_function(cos);
  215.   m.sqr();
  216.   m1.sqr();
  217.   m += m1;
  218.   verify_element_value(m,1);
  219.   
  220.   cout << "\nDone\n\n";
  221. }
  222.  
  223. /*
  224.  *------------------------------------------------------------------------
  225.  *          Test binary matrix element-by-element operations
  226.  */
  227.  
  228. static void test_binary_ebe_op(const int rsize, const int csize)
  229. {
  230.   const double pattern = 4.25;
  231.   register int i,j;
  232.  
  233.   cout << "\n---> Test Binary Matrix element-by-element operations\n";
  234.  
  235.   Matrix m(2,rsize+1,0,csize-1);
  236.   Matrix m1(m);
  237.   Matrix mp(m);
  238.  
  239.   for(i=mp.q_row_lwb(); i<=mp.q_row_upb(); i++)
  240.     for(j=mp.q_col_lwb(); j<=mp.q_col_upb(); j++)
  241.       mp(i,j) = (i-m.q_nrows()/2.)*j*pattern;
  242.   
  243.  
  244.   cout << "\nVerify assignment of a matrix to the matrix\n";
  245.   m = pattern;
  246.   m1.clear();
  247.   m1 = m;
  248.   verify_element_value(m1,pattern);
  249.   assert( m1 == m );
  250.  
  251.   cout << "\nAdding the matrix to itself, uniform pattern " << pattern << "\n";
  252.   m.clear(); m = pattern;
  253.   m1 = m; m1 += m1;
  254.   verify_element_value(m1,2*pattern);
  255.   cout << "  subtracting two matrices ...\n";
  256.   m1 -= m;
  257.   verify_element_value(m1,pattern);
  258.   cout << "  subtracting the matrix from itself\n";
  259.   m1 -= m1;
  260.   verify_element_value(m1,0);
  261.   cout << "  adding two matrices together\n";
  262.   m1 += m;
  263.   verify_element_value(m1,pattern);
  264.  
  265.   cout << "\nArithmetic operations on matrices with not the same elements\n";
  266.   cout << "   adding mp to the zero matrix...\n";
  267.   m.clear(); m += mp;
  268.   verify_matrix_identity(m,mp);
  269.   m1 = m;
  270.   cout << "   making m = 3*mp and m1 = 3*mp, via add() and succesive mult\n";
  271.   add(m,2,mp);
  272.   m1 += m1; m1 += mp;
  273.   verify_matrix_identity(m,m1);
  274.   cout << "   clear both m and m1, by subtracting from itself and via add()\n";
  275.   m1 -= m1;
  276.   add(m,-3,mp);
  277.   verify_matrix_identity(m,m1);
  278.  
  279.   cout << "\nTesting element-by-element multiplications and divisions\n";
  280.   cout << "   squaring each element with sqr() and via multiplication\n";
  281.   m = mp; m1 = mp;
  282.   m.sqr();
  283.   element_mult(m1,m1);
  284.   verify_matrix_identity(m,m1);
  285.   cout << "   compare (m = pattern^2)/pattern with pattern\n";
  286.   m = pattern; m1 = pattern;
  287.   m.sqr();
  288.   element_div(m,m1);
  289.   verify_element_value(m,pattern);
  290.   compare(m1,m,"Original vector and vector after squaring and dividing");
  291.  
  292.   cout << "\nDone\n";
  293. }
  294.  
  295. /*
  296.  *------------------------------------------------------------------------
  297.  *            Verify the determinant evaluation
  298.  */
  299.  
  300. static void test_determinant(const int msize)
  301. {
  302.   cout << "\n---> Verify determinant evaluation\n"
  303.           "for a square matrix of size " << msize << "\n";
  304.  
  305.   Matrix m(msize,msize);
  306.  
  307.   cout << "\nCheck to see that the determinant of the unit matrix is one";
  308.   m.unit_matrix();
  309.   cout << "\n    determinant is " << m.determinant();
  310.   assert( m.determinant() == 1 );
  311.  
  312.   const double pattern = 2.5;
  313.   cout << "\nCheck the determinant for the matrix with " << pattern <<
  314.           "\n    at the diagonal";
  315.   register int i,j;
  316.   for(i=m.q_row_lwb(); i<=m.q_row_upb(); i++)
  317.     for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
  318.       m(i,j) = ( i==j ? pattern : 0 );
  319.   cout << "\n    determinant is " << m.determinant();
  320.   assert( m.determinant() == pow(pattern,m.q_nrows()) );
  321.  
  322.   cout << "\nCheck the determinant for the matrix with " << pattern <<
  323.           "\n    at the anti-diagonal";
  324.   for(i=m.q_row_lwb(); i<=m.q_row_upb(); i++)
  325.     for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
  326.       m(i,j) = ( i==(m.q_col_upb()+m.q_col_lwb()-j) ? pattern : 0 );
  327.   assert( m.determinant() == pow(pattern,m.q_nrows()) * 
  328.      ( m.q_nrows()*(m.q_nrows()-1)/2 & 1 ? -1 : 1 ) );
  329.  
  330.   cout << "\nCheck the determinant for the singular matrix"
  331.           "\n    defined as above with zero first row";
  332.   m.clear();
  333.   for(i=m.q_row_lwb()+1; i<=m.q_row_upb(); i++)
  334.     for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
  335.       m(i,j) = ( i==(m.q_col_upb()+m.q_col_lwb()-j) ? pattern : 0 );
  336.   cout << "\n    determinant is " << m.determinant();
  337.   assert( m.determinant() == 0 );
  338.  
  339.   cout << "\nCheck out the determinant of the Hilbert matrix";
  340.   Matrix H(3,3);
  341.   H.hilbert_matrix();
  342.   cout << "\n    3x3 Hilbert matrix: exact determinant 1/2160 ";
  343.   cout << "\n                              computed    1/"<< 1/H.determinant();
  344.  
  345.   H.resize_to(4,4);
  346.   H.hilbert_matrix();
  347.   cout << "\n    4x4 Hilbert matrix: exact determinant 1/6048000 ";
  348.   cout << "\n                              computed    1/"<< 1/H.determinant();
  349.  
  350.   H.resize_to(5,5);
  351.   H.hilbert_matrix();
  352.   cout << "\n    5x5 Hilbert matrix: exact determinant 3.749295e-12";
  353.   cout << "\n                              computed    "<< H.determinant();
  354.  
  355.   H.resize_to(7,7);
  356.   H.hilbert_matrix();
  357.   cout << "\n    7x7 Hilbert matrix: exact determinant 4.8358e-25";
  358.   cout << "\n                              computed    "<< H.determinant();
  359.  
  360.   H.resize_to(9,9);
  361.   H.hilbert_matrix();
  362.   cout << "\n    9x9 Hilbert matrix: exact determinant 9.72023e-43";
  363.   cout << "\n                              computed    "<< H.determinant();
  364.  
  365.   H.resize_to(10,10);
  366.   H.hilbert_matrix();
  367.   cout << "\n    10x10 Hilbert matrix: exact determinant 2.16418e-53";
  368.   cout << "\n                              computed    "<< H.determinant();
  369.  
  370.   cout << "\nDone\n";
  371. }
  372.  
  373. /*
  374.  *------------------------------------------------------------------------
  375.  *            Verify the norm calculation
  376.  */
  377.  
  378. static void test_norms(const int rsize, const int csize)
  379. {
  380.   cout << "\n---> Verify norm calculations\n";
  381.  
  382.   const double pattern = 10.25;
  383.  
  384.   if( rsize % 2 == 1 || csize %2 == 1 )
  385.     _error("Sorry, size of the matrix to test must be even for this test\n");
  386.  
  387.   Matrix m(rsize,csize);
  388.   Matrix m1(m);
  389.  
  390.   cout << "\nAssign " << pattern << " to all the elements and check norms\n";
  391.   m = pattern;
  392.   cout << "  1. (col) norm should be pattern*nrows\n";
  393.   assert( m.norm_1() == pattern*m.q_nrows() );
  394.   assert( m.norm_1() == m.col_norm() );
  395.   cout << "  Inf (row) norm should be pattern*ncols\n";
  396.   assert( m.norm_inf() == pattern*m.q_ncols() );
  397.   assert( m.norm_inf() == m.row_norm() );
  398.   cout << "  Square of the Eucl norm has got to be pattern^2 * no_elems\n";
  399.   assert( m.e2_norm() == sqr(pattern)*m.q_no_elems() );
  400.   m1.clear();
  401.   assert( m.e2_norm() == e2_norm(m,m1) );
  402.  
  403. #if 0
  404.   double ap_step = 1;
  405.   double ap_a0   = -pattern;
  406.   int n = m.q_no_elems();
  407.   cout << "\nAssign the arithm progression with 1. term " << ap_a0 <<
  408.           "\nand the difference " << ap_step << "\n";
  409.   register int i;
  410.   for(i=v.q_lwb(); i<=v.q_upb(); i++)
  411.     v(i) = (i-v.q_lwb())*ap_step + ap_a0;
  412.   int l = min(max((int)ceil(-ap_a0/ap_step),0),n);
  413.   double norm = (2*ap_a0 + (l+n-1)*ap_step)/2*(n-l) +
  414.     (-2*ap_a0-(l-1)*ap_step)/2*l;
  415.   cout << "  1. norm should be " << norm << "\n";
  416.   assert( v.norm_1() == norm );
  417.   norm = n*( sqr(ap_a0) + ap_a0*ap_step*(n-1) + sqr(ap_step)*(n-1)*(2*n-1)/6);
  418.   cout << "  Square of the 2. norm has got to be "
  419.           "n*[ a0^2 + a0*q*(n-1) + q^2/6*(n-1)*(2n-1) ], or " << norm << "\n";
  420.   assert( v.norm_2_sqr() == norm );
  421.   norm = max(abs(v(v.q_lwb())),abs(v(v.q_upb())));
  422.   cout << "  Inf norm should be max(abs(a0),abs(a0+(n-1)*q)), ie " << norm <<
  423.           "\n";
  424.   assert( v.norm_inf() == norm );
  425.  
  426.   m1.clear();
  427.   compare(m,m1,"Compare the matrix m with a zero matrix");
  428.  
  429.   cout << "\nConstruct v1 to be orthogonal to v as v(n), -v(n-1), v(n-2)...\n";
  430.   for(i=0; i<v1.q_no_elems(); i++)
  431.     v1(i+v1.q_lwb()) = v(v.q_upb()-i) * ( i % 2 == 1 ? -1 : 1 );
  432.   cout << "||v1|| has got to be equal ||v|| regardless of the norm def\n";
  433.   assert( v1.norm_1() == v.norm_1() );
  434.   assert( v1.norm_2_sqr() == v.norm_2_sqr() );
  435.   assert( v1.norm_inf() == v.norm_inf() );
  436. #endif
  437.  
  438.   cout << "\nDone\n";
  439. }
  440.  
  441. /*
  442.  *------------------------------------------------------------------------
  443.  *                Root module
  444.  */
  445.  
  446. main()
  447. {
  448.   cout << "\n\n" << _Minuses << 
  449.           "\n\t\tVerify Operations on Matrices\n";
  450.  
  451. #if 1
  452.   test_allocation();
  453.   test_element_op(20,10);
  454.   test_binary_ebe_op(10,20);
  455.   test_norms(40,20);
  456. #endif
  457.   test_determinant(20);
  458. }
  459.  
  460.