home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / LinAlg 3.1 / LinAlg / matrix1.cc < prev    next >
Encoding:
Text File  |  1995-12-21  |  17.1 KB  |  691 lines  |  [TEXT/ALFA]
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *            Linear Algebra Package
  6.  *
  7.  *        Basic linear algebra operations, level 1
  8.  *              Element-wise operations
  9.  *
  10.  * $Id: matrix1.cc,v 3.0 1995/01/07 19:06:12 oleg Exp $
  11.  *
  12.  ************************************************************************
  13.  */
  14.  
  15. #pragma implementation "LinAlg.h"
  16.  
  17. #include "LinAlg.h"
  18. #include <math.h>
  19.  
  20.  
  21. /*
  22.  *------------------------------------------------------------------------
  23.  *            Constructors and destructors
  24.  */
  25.  
  26. void Matrix::allocate(
  27.     const int no_rows,        // No. of rows
  28.     const int no_cols,        // No. of cols
  29.     const int row_lwb,        // Row index lower bound
  30.     const int col_lwb        // Col index lower bound
  31. )
  32. {
  33.   valid_code = MATRIX_val_code;
  34.  
  35.   assure((nrows=no_rows) > 0, "No. of matrix cols has got to be positive");
  36.   assure((ncols=no_cols) > 0, "No. of matrix rows has got to be positive");
  37.  
  38.   Matrix::row_lwb = row_lwb;
  39.   Matrix::col_lwb = col_lwb;
  40.  
  41.   name = "";
  42.  
  43.   nelems = nrows * ncols;
  44.  
  45.   assert( (elements = (REAL *)calloc(nelems,sizeof(REAL))) != 0 );
  46.  
  47.   if( ncols == 1 )        // Only one col - index is dummy actually
  48.   {
  49.     index = &elements;
  50.     return;
  51.   }
  52.   
  53.   assert( (index  = (REAL **)calloc(ncols,sizeof(REAL *))) != 0 );
  54.   register int i;
  55.   register REAL * col_p;
  56.   for(i=0,col_p=&elements[0]; i<ncols; i++,col_p += nrows)
  57.     index[i] = col_p;
  58. }
  59.  
  60.  
  61. Matrix::~Matrix(void)        // Dispose the Matrix struct
  62. {
  63.   is_valid();
  64.   if( ncols != 1 )
  65.     free(index);
  66.   free(elements);
  67.   if( name[0] != '\0' )
  68.     delete name;
  69.   valid_code = 0;
  70. }
  71.  
  72.                 // Set a new Matrix name
  73. void Matrix::set_name(const char * new_name)
  74. {
  75.   if( name != 0 && name[0] != '\0' )    // Dispose of the previous matrix name
  76.     delete name;
  77.  
  78.   if( new_name == 0 || new_name[0] == '\0' )
  79.     name = "";                // Matrix is anonymous now
  80.   else
  81.     name = new char[strlen(new_name)+1], strcpy(name,new_name);
  82. }
  83.  
  84.                 // Erase the old matrix and create a
  85.                 // new one according to new boundaries
  86.                 // with indexation starting at 1
  87. void Matrix::resize_to(const int nrows, const int ncols)
  88. {
  89.   is_valid();
  90.   if( nrows == Matrix::nrows && ncols == Matrix::ncols )
  91.     return;
  92.  
  93.   if( ncols != 1 )
  94.     free(index);
  95.   free(elements);
  96.  
  97.   char * old_name = name;
  98.   allocate(nrows,ncols);
  99.   name = old_name;
  100. }
  101.                 // Erase the old matrix and create a
  102.                 // new one according to new boundaries
  103. void Matrix::resize_to(const int row_lwb, const int row_upb,
  104.                const int col_lwb, const int col_upb)
  105. {
  106.   is_valid();
  107.   const int new_nrows = row_upb-row_lwb+1;
  108.   const int new_ncols = col_upb-col_lwb+1;
  109.   Matrix::row_lwb = row_lwb;
  110.   Matrix::col_lwb = col_lwb;
  111.   if( new_nrows == Matrix::nrows && new_ncols == Matrix::ncols )
  112.     return;
  113.  
  114.   if( ncols != 1 )
  115.     free(index);
  116.   free(elements);
  117.  
  118.   char * old_name = name;
  119.   allocate(new_nrows,new_ncols,row_lwb,col_lwb);
  120.   name = old_name;
  121. }
  122.  
  123.                     // Routing constructor module
  124. Matrix::Matrix(const MATRIX_CREATORS_1op op, const Matrix& prototype)
  125. {
  126.   prototype.is_valid();
  127.   switch(op)
  128.   {
  129.     case Zero:
  130.          allocate(prototype.nrows,prototype.ncols,
  131.           prototype.row_lwb,prototype.col_lwb);
  132.      break;
  133.  
  134.     case Unit:
  135.          allocate(prototype.nrows,prototype.ncols,
  136.           prototype.row_lwb,prototype.col_lwb);
  137.      unit_matrix();
  138.      break;
  139.  
  140.     case Transposed:
  141.      _transpose(prototype);
  142.      break;
  143.  
  144.     case Inverted:
  145.      _invert(prototype);
  146.      break;
  147.  
  148.     default:
  149.      _error("Operation %d is not yet implemented",op);
  150.   }
  151. }
  152.  
  153.                     // Routing constructor module
  154. Matrix::Matrix(const Matrix& A, const MATRIX_CREATORS_2op op, const Matrix& B)
  155. {
  156.   A.is_valid();
  157.   B.is_valid();
  158.   switch(op)
  159.   {
  160.     case Mult:
  161.          _AmultB(A,B);
  162.      break;
  163.  
  164.     case TransposeMult:
  165.      _AtmultB(A,B);
  166.      break;
  167.  
  168.     default:
  169.      _error("Operation %d is not yet implemented",op);
  170.   }
  171. }
  172.  
  173. /*
  174.  *------------------------------------------------------------------------
  175.  *             Making a matrix of a special kind    
  176.  */
  177.  
  178.                 // Make a unit matrix
  179.                 // (Matrix needn't be a square one)
  180.                 // The matrix is traversed in the
  181.                 // natural (that is, col by col) order
  182. Matrix& Matrix::unit_matrix(void)
  183. {
  184.   is_valid();
  185.   register REAL *ep = elements;
  186.   register int i,j;
  187.  
  188.   for(j=0; j < ncols; j++)
  189.     for(i=0; i < nrows; i++)
  190.         *ep++ = ( i==j ? 1.0 : 0.0 );
  191.  
  192.   return *this;
  193. }
  194.  
  195.                 // Make a Hilbert matrix
  196.                 // Hilb[i,j] = 1/(i+j-1), i,j=1...max, OR
  197.                 // Hilb[i,j] = 1/(i+j+1), i,j=0...max-1
  198.                 // (Matrix needn't be a square one)
  199.                 // The matrix is traversed in the
  200.                 // natural (that is, col by col) order
  201. Matrix& Matrix::hilbert_matrix(void)
  202. {
  203.   is_valid();
  204.   register REAL *ep = elements;
  205.   register int i,j;
  206.  
  207.   for(j=0; j < ncols; j++)
  208.     for(i=0; i < nrows; i++)
  209.         *ep++ = 1./(i+j+1);
  210.  
  211.   return *this;
  212. }
  213.  
  214.             // Create an orthonormal (2^n)*(no_cols) Haar
  215.             // (sub)matrix, whose columns are Haar functions
  216.             // If no_cols is 0, create the complete matrix
  217.             // with 2^n columns
  218.             // E.g., the complete Haar matrix of the second order
  219.             // is
  220.             // column 1: [ 1  1  1  1]/2
  221.             // column 2: [ 1  1 -1 -1]/2
  222.             // column 3: [ 1 -1  0  0]/sqrt(2)
  223.             // column 4: [ 0  0  1 -1]/sqrt(2)
  224.             // Matrix m is assumed to be zero originally
  225. void _make_haar_matrix(Matrix& m)
  226. {
  227.   m.is_valid();
  228.   assert( m.ncols <= m.nrows && m.ncols > 0 );
  229.   register REAL * cp = m.elements;
  230.      const REAL * m_end = m.elements + m.nelems;
  231.   register int i;
  232.  
  233.   double norm_factor = 1/sqrt((double)m.nrows);
  234.  
  235.   for(i=0; i<m.nrows; i++)    // First column is always 1 (up to normali- 
  236.     *cp++ = norm_factor;    // zation)
  237.  
  238.                 // The other functions are kind of steps:
  239.                 // stretch of 1 followed by the equally
  240.                 // long stretch of -1
  241.                 // The functions can be grouped in families
  242.                 // according to their order (step size),
  243.                 // differing only in the location of the step
  244.   int step_length = m.nrows/2;
  245.   while( cp < m_end && step_length > 0 )
  246.   {
  247.     for(register int step_position=0; cp < m_end && step_position < m.nrows;
  248.     step_position += 2*step_length, cp += m.nrows)
  249.     {
  250.       register REAL * ccp = cp + step_position;
  251.       for(i=0; i<step_length; i++)
  252.     *ccp++ = norm_factor;
  253.       for(i=0; i<step_length; i++)
  254.     *ccp++ = -norm_factor;
  255.     }
  256.     step_length /= 2;
  257.     norm_factor *= sqrt(2.0);
  258.   }
  259.   assert( step_length != 0 || cp == m_end );
  260.   assert( m.nrows != m.ncols || step_length == 0 );
  261. }
  262.  
  263. haar_matrix::haar_matrix(const int order, const int no_cols)
  264.     : LazyMatrix(1<<order,no_cols == 0 ? 1<<order : no_cols)
  265. {
  266.   assert(order > 0 && no_cols >= 0);
  267. }
  268.  
  269. /*
  270.  *------------------------------------------------------------------------
  271.  *             Matrix-scalar arithmetics
  272.  *          Modify every element according to the operation
  273.  */
  274.  
  275.                 // For every element, do `elem OP value`
  276. #define COMPUTED_VAL_ASSIGNMENT(OP,VALTYPE)                \
  277.                                     \
  278. Matrix& Matrix::operator OP (const VALTYPE val)                \
  279. {                                    \
  280.   is_valid();                                \
  281.   register REAL * ep = elements;                    \
  282.   while( ep < elements+nelems )                        \
  283.     *ep++ OP val;                            \
  284.                                     \
  285.   return *this;                                \
  286. }                                    \
  287.  
  288. COMPUTED_VAL_ASSIGNMENT(=,REAL)
  289. COMPUTED_VAL_ASSIGNMENT(+=,double)
  290. COMPUTED_VAL_ASSIGNMENT(-=,double)
  291. COMPUTED_VAL_ASSIGNMENT(*=,double)
  292.  
  293. #undef COMPUTED_VAL_ASSIGNMENT
  294.  
  295.  
  296.                 // is "element OP val" true for all
  297.                 // matrix elements?
  298.  
  299. #define COMPARISON_WITH_SCALAR(OP)                    \
  300.                                     \
  301. bool Matrix::operator OP (const REAL val) const                \
  302. {                                    \
  303.   is_valid();                                \
  304.   register REAL * ep = elements;                    \
  305.   while( ep < elements + nelems )                    \
  306.     if( !(*ep++ OP val) )                        \
  307.       return false;                            \
  308.                                     \
  309.   return true;                                \
  310. }                                    \
  311.  
  312.  
  313. COMPARISON_WITH_SCALAR(==)
  314. COMPARISON_WITH_SCALAR(!=)
  315. COMPARISON_WITH_SCALAR(<)
  316. COMPARISON_WITH_SCALAR(<=)
  317. COMPARISON_WITH_SCALAR(>)
  318. COMPARISON_WITH_SCALAR(>=)
  319.  
  320. #undef COMPARISON_WITH_SCALAR
  321.  
  322. /*
  323.  *------------------------------------------------------------------------
  324.  *        Apply algebraic functions to all the matrix elements
  325.  */
  326.  
  327.                 // Take an absolute value of a matrix
  328. Matrix& Matrix::abs(void)
  329. {
  330.   is_valid();
  331.   register REAL * ep;
  332.   for(ep=elements; ep < elements+nelems; ep++)
  333.     *ep = ::abs(*ep);
  334.  
  335.   return *this;
  336. }
  337.  
  338.                 // Square each element
  339. Matrix& Matrix::sqr(void)
  340. {
  341.   is_valid();
  342.   register REAL * ep;
  343.   for(ep=elements; ep < elements+nelems; ep++)
  344.     *ep = *ep * *ep;
  345.  
  346.   return *this;
  347. }
  348.  
  349.                 // Take a square root of all the elements
  350. Matrix& Matrix::sqrt(void)
  351. {
  352.   is_valid();
  353.   register REAL * ep;
  354.   for(ep=elements; ep < elements+nelems; ep++)
  355.     if( *ep >= 0 )
  356.       *ep = ::sqrt(*ep);
  357.     else
  358.       info(),
  359.       _error("(%d,%d)-th element, %g, is negative. Can't take the square root",
  360.          (ep-elements) % nrows + row_lwb,
  361.          (ep-elements) / nrows + col_lwb, *ep );
  362.  
  363.   return *this;
  364. }
  365.  
  366.                 // Apply a user-defined action to each matrix
  367.                 // element. The matrix is traversed in the
  368.                 // natural (that is, col by col) order
  369. Matrix& Matrix::apply(ElementAction& action)
  370. {
  371.   is_valid();
  372.   register REAL * ep = elements;
  373.   for(action.j=col_lwb; action.j<col_lwb+ncols; action.j++)
  374.     for(action.i=row_lwb; action.i<row_lwb+nrows; action.i++)
  375.       action.operation(*ep++);
  376.   assert( ep == elements+nelems );
  377.  
  378.   return *this;
  379. }
  380.  
  381. /*
  382.  *------------------------------------------------------------------------
  383.  *            Matrix-Matrix element-wise operations
  384.  */
  385.  
  386.                 // Check to see if two matrices are identical
  387. bool operator == (const Matrix& im1, const Matrix& im2)
  388. {
  389.   are_compatible(im1,im2);
  390.   return (memcmp(im1.elements,im2.elements,im1.nelems*sizeof(REAL)) == 0);
  391. }
  392.  
  393.                 // Add the source to the target
  394. Matrix& operator += (Matrix& target, const Matrix& source)
  395. {
  396.   are_compatible(target,source);
  397.  
  398.   register REAL * sp = source.elements;
  399.   register REAL * tp = target.elements;
  400.   for(; tp < target.elements+target.nelems;)
  401.     *tp++ += *sp++;
  402.   
  403.   return target;
  404. }
  405.   
  406.                 // Subtract the source from the target
  407. Matrix& operator -= (Matrix& target, const Matrix& source)
  408. {
  409.   are_compatible(target,source);
  410.   register REAL * sp = source.elements;
  411.   register REAL * tp = target.elements;
  412.   for(; tp < target.elements+target.nelems;)
  413.     *tp++ -= *sp++;
  414.   
  415.   return target;
  416. }
  417.  
  418.                 // Modified addition
  419.                 //    Target += scalar*Source
  420. Matrix& add(Matrix& target, const double scalar,const Matrix& source)
  421. {
  422.   are_compatible(target,source);
  423.  
  424.   register REAL * sp = source.elements;
  425.   register REAL * tp = target.elements;
  426.   for(; tp < target.elements+target.nelems;)
  427.     *tp++ += scalar * *sp++;
  428.   
  429.   return target;
  430. }
  431.  
  432.                 // Multiply target by the source
  433.                 // element-by-element
  434. Matrix& element_mult(Matrix& target, const Matrix& source)
  435. {
  436.   are_compatible(target,source);
  437.   register REAL * sp = source.elements;
  438.   register REAL * tp = target.elements;
  439.   for(; tp < target.elements+target.nelems;)
  440.     *tp++ *= *sp++;
  441.   
  442.   return target;
  443. }
  444.  
  445.                 // Divide target by the source
  446.                 // element-by-element
  447. Matrix& element_div(Matrix& target, const Matrix& source)
  448. {
  449.   are_compatible(target,source);
  450.   register REAL * sp = source.elements;
  451.   register REAL * tp = target.elements;
  452.   for(; tp < target.elements+target.nelems;)
  453.     *tp++ /= *sp++;
  454.   
  455.   return target;
  456. }
  457.  
  458. /*
  459.  *------------------------------------------------------------------------
  460.  *            Compute matrix norms
  461.  */
  462.  
  463.                 // Row matrix norm
  464.                 // MAX{ SUM{ |M(i,j)|, over j}, over i}
  465.                 // The norm is induced by the infinity
  466.                 // vector norm
  467. double Matrix::row_norm(void) const
  468. {
  469.   is_valid();
  470.   register REAL * ep = elements;
  471.   register double norm = 0;
  472.  
  473.   while(ep < elements+nrows)        // Scan the matrix row-after-row
  474.   {
  475.     register int j;
  476.     register double sum = 0;
  477.     for(j=0; j<ncols; j++,ep+=nrows)    // Scan a row to compute the sum
  478.       sum += ::abs(*ep);
  479.     ep -= nelems - 1;            // Point ep to the beg of the next row
  480.     norm = ::max(norm,sum);
  481.   }
  482.   assert( ep == elements + nrows );
  483.  
  484.   return norm;
  485. }
  486.  
  487.                 // Column matrix norm
  488.                 // MAX{ SUM{ |M(i,j)|, over i}, over j}
  489.                 // The norm is induced by the 1.
  490.                 // vector norm
  491. double Matrix::col_norm(void) const
  492. {
  493.   is_valid();
  494.   register REAL * ep = elements;
  495.   register double norm = 0;
  496.  
  497.   while(ep < elements+nelems)        // Scan the matrix col-after-col
  498.   {                    // (i.e. in the natural order of elems)
  499.     register int i;
  500.     register double sum = 0;
  501.     for(i=0; i<nrows; i++)        // Scan a col to compute the sum
  502.       sum += ::abs(*ep++);
  503.     norm = ::max(norm,sum);
  504.   }
  505.   assert( ep == elements + nelems );
  506.  
  507.   return norm;
  508. }
  509.  
  510.  
  511.                 // Square of the Euclidian norm
  512.                 // SUM{ m(i,j)^2 }
  513. double Matrix::e2_norm(void) const
  514. {
  515.   is_valid();
  516.   register REAL * ep;
  517.   register double sum = 0;
  518.   for(ep=elements; ep < elements+nelems;)
  519.     sum += ::sqr(*ep++);
  520.  
  521.   return sum;
  522. }
  523.  
  524.                 // Square of the Euclidian norm of the
  525.                 // difference between two matrices
  526. double e2_norm(const Matrix& m1, const Matrix& m2)
  527. {
  528.   are_compatible(m1,m2);
  529.   register REAL * mp1 = m1.elements;
  530.   register REAL * mp2 = m2.elements;
  531.   register double sum = 0;
  532.   for(; mp1 < m1.elements+m1.nelems;)
  533.     sum += sqr(*mp1++ - *mp2++);
  534.  
  535.   return sum;
  536. }
  537.  
  538. /*
  539.  *------------------------------------------------------------------------
  540.  *             Some service operations
  541.  */
  542.  
  543. void Matrix::info(void) const    // Print some information about the matrix
  544. {
  545.   is_valid();
  546.   message("\nMatrix %d:%dx%d:%d '%s'",row_lwb,nrows+row_lwb-1,
  547.       col_lwb,ncols+col_lwb-1,name);
  548. }
  549.  
  550.                 // Print the Matrix as a table of elements
  551.                 // (zeros are printed as dots)
  552. void Matrix::print(const char * title) const
  553. {
  554.   is_valid();
  555.   message("\nMatrix %dx%d '%s' is as follows",nrows,ncols,title);
  556.  
  557.   const int cols_per_sheet = 6;
  558.   register int sheet_counter;
  559.   
  560.   for(sheet_counter=1; sheet_counter<=ncols; sheet_counter +=cols_per_sheet)
  561.   {
  562.     message("\n\n     |");
  563.     register int i,j;
  564.     for(j=sheet_counter; j<sheet_counter+cols_per_sheet && j<=ncols; j++)
  565.       message("   %6d  |",j+col_lwb-1);
  566.     message("\n%s\n",_Minuses);
  567.     for(i=1; i<=nrows; i++)
  568.     {
  569.       message("%4d |",i+row_lwb-1);
  570.       for(j=sheet_counter; j<sheet_counter+cols_per_sheet && j<=ncols; j++)
  571.     message("%11.4g  ",(*this)(i+row_lwb-1,j+col_lwb-1));
  572.       message("\n");
  573.     }
  574.   }
  575.   message("Done\n");
  576. }
  577.  
  578. //#include <builtin.h>
  579. void compare(            // Compare the two Matrixs
  580.     const Matrix& matrix1,    // and print out the result of comparison
  581.     const Matrix& matrix2,
  582.     const char * title )
  583. {
  584.   register int i,j;
  585.  
  586.   are_compatible(matrix1,matrix2);
  587.  
  588.   message("\n\nComparison of two Matrices:\n\t%s",title);
  589.   matrix1.info();
  590.   matrix2.info();
  591.  
  592.   double norm1 = 0, norm2 = 0;        // Norm of the Matrixs
  593.   double ndiff = 0;            // Norm of the difference
  594.   int imax=0,jmax=0;            // For the elements that differ most
  595.   REAL difmax = -1;
  596.   register REAL *mp1 = matrix1.elements; // Matrix element pointers
  597.   register REAL *mp2 = matrix2.elements;
  598.  
  599.   for(j=0; j < matrix1.ncols; j++)    // Due to the column-wise arrangement,
  600.     for(i=0; i < matrix1.nrows; i++)    // the row index changes first
  601.     {
  602.       REAL mv1 = *mp1++;
  603.       REAL mv2 = *mp2++;
  604.       REAL diff = abs(mv1-mv2);
  605.  
  606.       if( diff > difmax )
  607.       {
  608.     difmax = diff;
  609.     imax = i;
  610.     jmax = j;
  611.       }
  612.       norm1 += abs(mv1);
  613.       norm2 += abs(mv2);
  614.       ndiff += abs(diff);
  615.     }
  616.  
  617.   imax += matrix1.row_lwb, jmax += matrix1.col_lwb;
  618.   message("\nMaximal discrepancy    \t\t%g",difmax);
  619.   message("\n   occured at the point\t\t(%d,%d)",imax,jmax);
  620.   const REAL mv1 = matrix1(imax,jmax);
  621.   const REAL mv2 = matrix2(imax,jmax);
  622.   message("\n Matrix 1 element is    \t\t%g",mv1);
  623.   message("\n Matrix 2 element is    \t\t%g",mv2);
  624.   message("\n Absolute error v2[i]-v1[i]\t\t%g",mv2-mv1);
  625.   message("\n Relative error\t\t\t\t%g\n",
  626.      (mv2-mv1)/max(abs(mv2+mv1)/2,1e-7) );
  627.  
  628.   message("\n||Matrix 1||   \t\t\t%g",norm1);
  629.   message("\n||Matrix 2||   \t\t\t%g",norm2);
  630.   message("\n||Matrix1-Matrix2||\t\t\t\t%g",ndiff);
  631.   message("\n||Matrix1-Matrix2||/sqrt(||Matrix1|| ||Matrix2||)\t%g\n\n",
  632.       ndiff/max( sqrt(norm1*norm2), 1e-7 )         );
  633.  
  634. }
  635.  
  636. /*
  637.  *------------------------------------------------------------------------
  638.  *            Service validation functions
  639.  */
  640.  
  641. void verify_element_value(const Matrix& m,const REAL val)
  642. {
  643.   register imax = 0, jmax = 0;
  644.   register double max_dev = 0;
  645.   register int i,j;
  646.   for(i=m.q_row_lwb(); i<=m.q_row_upb(); i++)
  647.     for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
  648.     {
  649.       register double dev = abs(m(i,j)-val);
  650.       if( dev >= max_dev )
  651.     imax = i, jmax = j, max_dev = dev;
  652.     }
  653.  
  654.   if( max_dev == 0 )
  655.     return;
  656.   else if( max_dev < 1e-5 )
  657.     message("Element (%d,%d) with value %g differs the most from what\n"
  658.         "was expected, %g, though the deviation %g is small\n",
  659.         imax,jmax,m(imax,jmax),val,max_dev);
  660.   else
  661.     _error("A significant difference from the expected value %g\n"
  662.        "encountered for element (%d,%d) with value %g",
  663.        val,imax,jmax,m(imax,jmax));
  664. }
  665.  
  666. void verify_matrix_identity(const Matrix& m1, const Matrix& m2)
  667. {
  668.   register imax = 0, jmax = 0;
  669.   register double max_dev = 0;
  670.   register int i,j;
  671.   are_compatible(m1,m2);
  672.   for(i=m1.q_row_lwb(); i<=m1.q_row_upb(); i++)
  673.     for(j=m1.q_col_lwb(); j<=m1.q_col_upb(); j++)
  674.     {
  675.       register double dev = abs(m1(i,j)-m2(i,j));
  676.       if( dev >= max_dev )
  677.     imax = i, jmax = j, max_dev = dev;
  678.     }
  679.  
  680.   if( max_dev == 0 )
  681.     return;
  682.   if( max_dev < 1e-5 )
  683.     message("Two (%d,%d) elements of matrices with values %g and %g\n"
  684.         "differ the most, though the deviation %g is small\n",
  685.         imax,jmax,m1(imax,jmax),m2(imax,jmax),max_dev);
  686.   else
  687.     _error("A significant difference between the matrices encountered\n"
  688.        "at (%d,%d) element, with values %g and %g",
  689.        imax,jmax,m1(imax,jmax),m2(imax,jmax));
  690. }
  691.