Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

matrix1.cpp

Go to the documentation of this file.
00001 // This may look like C code, but it is really -*- C++ -*-
00002 /*
00003  ************************************************************************
00004  *
00005  *                      Linear Algebra Package
00006  *
00007  *              Basic linear algebra operations, level 1
00008  *                    Element-wise operations
00009  *
00010  * $Id: matrix1.cpp,v 1.1 2004/05/21 21:02:52 maxx Exp $
00011  *
00012  ************************************************************************
00013  */
00014 
00015 #ifdef __GNUC__
00016 #pragma implementation "LinAlg.h"
00017 #endif
00018 
00019 #include "LAStreams.h"
00020 #include <math.h>
00021 #include "iostream.h"
00022 //#include <builtin.h>
00023 
00024 namespace linalg 
00025 {
00026     using namespace linalg;
00027     
00028 /*
00029  *------------------------------------------------------------------------
00030  *                      Constructors and destructors
00031  */
00032 
00033 void Matrix::allocate(void)
00034 {
00035   valid_code = MATRIX_val_code;
00036 
00037   assure( (nelems = nrows * ncols) > 0,
00038         "The number of matrix cols and rows has got to be positive");
00039 
00040   assure( !ref_counter.q_engaged(),
00041         "An attempt to allocate Matrix data which are in use" );
00042   
00043   name = "";
00044 
00045   assert( (elements = (REAL *)calloc(nelems,sizeof(REAL))) != 0 );
00046 
00047 }
00048 
00049 
00050 Matrix::~Matrix(void)           // Dispose of the Matrix struct
00051 {
00052   is_valid();
00053   free(elements);
00054   if( name[0] != '\0' )
00055     delete const_cast<char*>(name);
00056   valid_code = 0;
00057 }
00058 
00059                                 // Set a new Matrix name
00060 void Matrix::set_name(const char * new_name)
00061 {
00062   if( name != 0 && name[0] != '\0' )    // Dispose of the previous matrix name
00063     delete const_cast<char*>(name);
00064 
00065   if( new_name == 0 || new_name[0] == '\0' )
00066     name = "";                          // Matrix is anonymous now
00067   else
00068     name = new char[strlen(new_name)+1],
00069       strcpy(const_cast<char *>(name), new_name);
00070 }
00071 
00072                                 // Erase the old matrix and create a
00073                                 // new one according to new boundaries
00074                                 // with indexation starting at 1
00075 void Matrix::resize_to(const int nrows, const int ncols)
00076 {
00077   is_valid();
00078   if( nrows == Matrix::nrows && ncols == Matrix::ncols )
00079     return;
00080 
00081 #if 0
00082   if( ncols != 1 )
00083     free(index);
00084 #endif
00085   assert( !ref_counter.q_engaged() );
00086   free(elements);
00087 
00088   const char * const old_name = name;
00089   Matrix::nrows = nrows;
00090   Matrix::ncols = ncols;
00091   allocate();
00092   name = old_name;
00093 }
00094                                 // Erase the old matrix and create a
00095                                 // new one according to new boundaries
00096 void Matrix::resize_to(const DimSpec& dimension_specs)
00097 {
00098   is_valid();
00099   Matrix::row_lwb = dimension_specs.q_row_lwb();
00100   Matrix::col_lwb = dimension_specs.q_col_lwb();
00101   resize_to(dimension_specs.q_nrows(),dimension_specs.q_ncols());
00102 }
00103 
00104 #if 0
00105                                         // Routing constructor module
00106 Matrix::Matrix(const Matrix& A, const MATRIX_CREATORS_2op op, const Matrix& B)
00107 {
00108   A.is_valid();
00109   B.is_valid();
00110   switch(op)
00111   {
00112     case Mult:
00113          _AmultB(A,B);
00114          break;
00115 
00116     case TransposeMult:
00117          _AtmultB(A,B);
00118          break;
00119 
00120     default:
00121          _error("Operation %d is not yet implemented",op);
00122   }
00123 }
00124 #endif
00125 
00126 
00127                         // Build a column index to facilitate direct
00128                         // access to matrix elements
00129 REAL * const * MatrixDABase::build_index(const Matrix& m)
00130 {
00131   m.is_valid();
00132   if( m.ncols == 1 )            // Only one col - index is dummy actually
00133     return const_cast<REAL * const *>(&m.elements);
00134 
00135   REAL ** const index = (REAL **)calloc(m.ncols,sizeof(REAL *));
00136   assert( index != 0 );
00137   register REAL * col_p = &m.elements[0];
00138   for(register REAL** ip = index; ip<index+m.ncols; col_p += m.nrows)
00139     *ip++ = col_p;
00140   return index;
00141 }
00142 
00143 MatrixDABase::~MatrixDABase(void)
00144 {
00145   if( ncols != 1 )
00146     free((void *)index);
00147 }
00148 
00149 /*
00150  *------------------------------------------------------------------------
00151  *                  Making a matrix of a special kind   
00152  */
00153 
00154                                 // Make a unit matrix
00155                                 // (Matrix needn't be a square one)
00156                                 // The matrix is traversed in the
00157                                 // natural (that is, col by col) order
00158 Matrix& Matrix::unit_matrix(void)
00159 {
00160   is_valid();
00161   register REAL *ep = elements;
00162 
00163   for(register int j=0; j < ncols; j++)
00164     for(register int i=0; i < nrows; i++)
00165         *ep++ = ( i==j ? 1.0 : 0.0 );
00166 
00167   return *this;
00168 }
00169 
00170                                 // Make a Hilbert matrix
00171                                 // Hilb[i,j] = 1/(i+j-1), i,j=1...max, OR
00172                                 // Hilb[i,j] = 1/(i+j+1), i,j=0...max-1
00173                                 // (Matrix needn't be a square one)
00174                                 // The matrix is traversed in the
00175                                 // natural (that is, col by col) order
00176 Matrix& Matrix::hilbert_matrix(void)
00177 {
00178   is_valid();
00179   register REAL *ep = elements;
00180 
00181   for(register int j=0; j < ncols; j++)
00182     for(register int i=0; i < nrows; i++)
00183         *ep++ = 1./(i+j+1);
00184 
00185   return *this;
00186 }
00187 
00188                         // Create an orthonormal (2^n)*(no_cols) Haar
00189                         // (sub)matrix, whose columns are Haar functions
00190                         // If no_cols is 0, create the complete matrix
00191                         // with 2^n columns
00192                         // E.g., the complete Haar matrix of the second order
00193                         // is
00194                         // column 1: [ 1  1  1  1]/2
00195                         // column 2: [ 1  1 -1 -1]/2
00196                         // column 3: [ 1 -1  0  0]/sqrt(2)
00197                         // column 4: [ 0  0  1 -1]/sqrt(2)
00198                         // Matrix m is assumed to be zero originally
00199 void haar_matrix::fill_in(Matrix& m) const
00200 {
00201   m.is_valid();
00202   assert( m.ncols <= m.nrows && m.ncols > 0 );
00203   register REAL * cp = m.elements;
00204      const REAL * const m_end = m.elements + m.nelems;
00205 
00206   double norm_factor = 1/sqrt((double)m.nrows);
00207 
00208   for(register int i=0; i<m.nrows; i++) // First column is always 1
00209     *cp++ = norm_factor;        // (up to a normalization)
00210 
00211                                 // The other functions are kind of steps:
00212                                 // stretch of 1 followed by the equally
00213                                 // long stretch of -1
00214                                 // The functions can be grouped in families
00215                                 // according to their order (step size),
00216                                 // differing only in the location of the step
00217   int step_length = m.nrows/2;
00218   while( cp < m_end && step_length > 0 )
00219   {
00220     for(register int step_position=0; cp < m_end && step_position < m.nrows;
00221         step_position += 2*step_length, cp += m.nrows)
00222     {
00223       register REAL * ccp = cp + step_position;
00224       for(register int i=0; i<step_length; i++)
00225         *ccp++ = norm_factor;
00226       for(register int j=0; j<step_length; j++)
00227         *ccp++ = -norm_factor;
00228     }
00229     step_length /= 2;
00230     norm_factor *= sqrt(2.0);
00231   }
00232   assert( step_length != 0 || cp == m_end );
00233   assert( m.nrows != m.ncols || step_length == 0 );
00234 }
00235 
00236 haar_matrix::haar_matrix(const int order, const int no_cols)
00237     : LazyMatrix(1<<order,no_cols == 0 ? 1<<order : no_cols)
00238 {
00239   assert(order > 0 && no_cols >= 0);
00240 }
00241 
00242 /*
00243  *------------------------------------------------------------------------
00244  *                      Collection-scalar arithmetics
00245  *      walk through the collection and do something with each
00246  *                      element in turn, and the scalar
00247  */
00248 
00249                                 // For every element, do `elem OP value`
00250 #define COMPUTED_VAL_ASSIGNMENT(OP,VALTYPE)                             \
00251                                                                         \
00252 void ElementWise::operator OP (const VALTYPE val)                       \
00253 {                                                                       \
00254   register REAL * ep = start_ptr;                                       \
00255   while( ep < end_ptr )                                                 \
00256     *ep++ OP val;                                                       \
00257 }                                                                       \
00258 
00259 COMPUTED_VAL_ASSIGNMENT(=,REAL)
00260 COMPUTED_VAL_ASSIGNMENT(+=,double)
00261 COMPUTED_VAL_ASSIGNMENT(-=,double)
00262 COMPUTED_VAL_ASSIGNMENT(*=,double)
00263 
00264 #undef COMPUTED_VAL_ASSIGNMENT
00265 
00266 
00267                                 // is "element OP val" true for all
00268                                 // elements in the collection?
00269 
00270 #define COMPARISON_WITH_SCALAR(OP)                                      \
00271                                                                         \
00272 bool ElementWiseConst::operator OP (const REAL val) const               \
00273 {                                                                       \
00274   register const REAL * ep = start_ptr;                                 \
00275   while( ep < end_ptr )                                                 \
00276     if( !(*ep++ OP val) )                                               \
00277       return false;                                                     \
00278                                                                         \
00279   return true;                                                          \
00280 }                                                                       \
00281 
00282 
00283 COMPARISON_WITH_SCALAR(==)
00284 COMPARISON_WITH_SCALAR(!=)
00285 COMPARISON_WITH_SCALAR(<)
00286 COMPARISON_WITH_SCALAR(<=)
00287 COMPARISON_WITH_SCALAR(>)
00288 COMPARISON_WITH_SCALAR(>=)
00289 
00290 #undef COMPARISON_WITH_SCALAR
00291 
00292 /*
00293  *------------------------------------------------------------------------
00294  *      Apply algebraic functions to all elements of a collection
00295  */
00296 
00297                                 // Take an absolute value of a matrix
00298 void ElementWise::abs(void)
00299 {
00300   for(register REAL * ep=start_ptr; ep < end_ptr; ep++)
00301     *ep = fabs(*ep);
00302 }
00303 
00304                                 // Square each element
00305 void ElementWise::sqr(void)
00306 {
00307   for(register REAL * ep=start_ptr; ep < end_ptr; ep++)
00308     *ep = *ep * *ep;
00309 }
00310 
00311                                 // Take the square root of all the elements
00312 void ElementWise::sqrt(void)
00313 {
00314   for(register REAL * ep=start_ptr; ep < end_ptr; ep++)
00315     if( *ep >= 0 )
00316       *ep = ::sqrt(*ep);
00317     else
00318       _error("%d-th element, %g, is negative. Can't take the square root",
00319              (ep-start_ptr), *ep );
00320 }
00321 
00322                                 // Apply a user-defined action to each matrix
00323                                 // element. The matrix is traversed in the
00324                                 // natural (that is, col by col) order
00325 Matrix& Matrix::apply(ElementAction& action)
00326 {
00327   is_valid();
00328   register REAL * ep = elements;
00329   for(register int j=col_lwb; j<col_lwb+ncols; j++)
00330     for(register int i=row_lwb; i<row_lwb+nrows; i++)
00331       action.operation(*ep++,i,j);
00332   assert( ep == elements+nelems );
00333 
00334   return *this;
00335 }
00336 
00337 const Matrix& Matrix::apply(ConstElementAction& action) const
00338 {
00339   is_valid();
00340   register const REAL * ep = elements;
00341   for(register int j=col_lwb; j<col_lwb+ncols; j++)
00342     for(register int i=row_lwb; i<row_lwb+nrows; i++)
00343       action.operation(*ep++,i,j);
00344   assert( ep == elements+nelems );
00345 
00346   return *this;
00347 }
00348 
00349 
00350 /*
00351  *------------------------------------------------------------------------
00352  *              Element-wise operations on two groups of elements
00353  */
00354 
00355                                 // Check to see if two matrices are identical
00356 bool Matrix::operator == (const Matrix& m2) const
00357 {
00358   are_compatible(*this,m2);
00359   return (memcmp(elements,m2.elements,nelems*sizeof(REAL)) == 0);
00360 }
00361 
00362                                 // For every element, do `elem OP another.elem`
00363 #define TWO_GROUP_COMP(OP)                                      \
00364                                                                         \
00365 bool ElementWiseConst::operator OP (const ElementWiseConst& another) const \
00366 {                                                                       \
00367   sure_compatible_with(another);                                        \
00368   register const REAL * sp = another.start_ptr;                         \
00369   register const REAL * tp = start_ptr;                                 \
00370   while( tp < end_ptr )                                                 \
00371     if( !(*tp++ OP *sp++) )                                             \
00372       return false;                                                     \
00373                                                                         \
00374   return true;                                                          \
00375 }                                                                       \
00376 
00377 TWO_GROUP_COMP(==)
00378 TWO_GROUP_COMP(!=)
00379 TWO_GROUP_COMP(<)
00380 TWO_GROUP_COMP(<=)
00381 TWO_GROUP_COMP(>)
00382 TWO_GROUP_COMP(>=)
00383 
00384 #undef TWO_GROUP_COMP
00385 
00386                                 // For every element, do `elem OP another.elem`
00387 #define TWO_GROUP_OP(OP)                                        \
00388                                                                         \
00389 void ElementWise::operator OP (const ElementWiseConst& another)         \
00390 {                                                                       \
00391   sure_compatible_with(another);                                        \
00392   register const REAL * sp = another.start_ptr;                         \
00393   register REAL * tp = start_ptr;                                       \
00394   while( tp < end_ptr )                                                 \
00395     *tp++ OP *sp++;                                                     \
00396 }                                                                       \
00397 
00398 TWO_GROUP_OP(=)
00399 TWO_GROUP_OP(+=)
00400 TWO_GROUP_OP(-=)
00401 TWO_GROUP_OP(*=)
00402 TWO_GROUP_OP(/=)
00403 
00404 #undef TWO_GROUP_OP
00405 
00406 #if 0
00407                                 // Modified addition
00408                                 //      Target += scalar*Source
00409 Matrix& add(Matrix& target, const double scalar,const Matrix& source)
00410 {
00411   are_compatible(target,source);
00412 
00413   register REAL * sp = source.elements;
00414   register REAL * tp = target.elements;
00415   for(; tp < target.elements+target.nelems;)
00416     *tp++ += scalar * *sp++;
00417   
00418   return target;
00419 }
00420 #endif
00421 
00422 /*
00423  *------------------------------------------------------------------------
00424  *      Reduce a collection or a difference between two collections
00425  *              to a single number: a "norm"
00426  */
00427 
00428 #define REDUCE_SUM(X,VAL) X += (VAL)
00429 #define REDUCE_SUMSQ(X,VAL) X += sqr(VAL)
00430 #define REDUCE_SUMABS(X,VAL) X += fabs(VAL)
00431 #define REDUCE_MAXABS(X,VAL) X = max((REAL)X,fabs(VAL))
00432 
00433 #define REDUCE_ONE(NAME,OP)                                                     \
00434                                                                         \
00435 double ElementWiseConst::NAME (void) const                              \
00436 {                                                                       \
00437   register double norm = 0;                                             \
00438   register const REAL * ep = start_ptr;                                 \
00439   while( ep < end_ptr )                                                 \
00440     OP(norm,*ep++);                                                     \
00441   return norm;                                                          \
00442 }                                                                       \
00443 
00444 
00445 
00446 REDUCE_ONE(sum,REDUCE_SUM)
00447 REDUCE_ONE(sum_squares,REDUCE_SUMSQ)
00448 REDUCE_ONE(sum_abs,REDUCE_SUMABS)
00449 REDUCE_ONE(max_abs,REDUCE_MAXABS)
00450 
00451 #undef REDUCE_ONE
00452 
00453 #define REDUCE_DIFF_OF_TWO(NAME,OP)                                     \
00454                                                                         \
00455 double ElementWiseConst::NAME (const ElementWiseConst& another) const   \
00456 {                                                                       \
00457   sure_compatible_with(another);                                        \
00458   register double norm = 0;                                             \
00459   register const REAL * sp = another.start_ptr;                         \
00460   register const REAL * tp = start_ptr;                                 \
00461   while( tp < end_ptr )                                                 \
00462     OP(norm,*tp++ - *sp++);                                             \
00463   return norm;                                                          \
00464 }                                                                       \
00465 
00466 REDUCE_DIFF_OF_TWO(sum_squares,REDUCE_SUMSQ)
00467 REDUCE_DIFF_OF_TWO(sum_abs,REDUCE_SUMABS)
00468 REDUCE_DIFF_OF_TWO(max_abs,REDUCE_MAXABS)
00469 
00470 #undef REDUCE_DIFF_OF_TWO
00471 
00472 #undef REDUCE_SUM
00473 #undef REDUCE_SUMSQ
00474 #undef REDUCE_SUMABS
00475 #undef REDUCE_MAXABS
00476 
00477 
00478 /*
00479  *------------------------------------------------------------------------
00480  *                      Compute matrix norms
00481  */
00482 
00483                                 // Row matrix norm
00484                                 // MAX{ SUM{ |M(i,j)|, over j}, over i}
00485                                 // The norm is induced by the infinity
00486                                 // vector norm
00487 double Matrix::row_norm(void) const
00488 {
00489   is_valid();
00490   register const REAL * ep = elements;
00491   register double norm = 0;
00492 
00493   while(ep < elements+nrows)            // Scan the matrix row-after-row
00494   {
00495     register int j;
00496     register double sum = 0;
00497     for(j=0; j<ncols; j++,ep+=nrows)    // Scan a row to compute the sum
00498       sum += fabs(*ep);
00499     ep -= nelems - 1;                   // Point ep to the beg of the next row
00500     norm = max(norm,sum);
00501   }
00502   assert( ep == elements + nrows );
00503 
00504   return norm;
00505 }
00506 
00507                                 // Column matrix norm
00508                                 // MAX{ SUM{ |M(i,j)|, over i}, over j}
00509                                 // The norm is induced by the 1.
00510                                 // vector norm
00511 double Matrix::col_norm(void) const
00512 {
00513   is_valid();
00514   register const REAL * ep = elements;
00515   register double norm = 0;
00516 
00517   while(ep < elements+nelems)           // Scan the matrix col-after-col
00518   {                                     // (i.e. in the natural order of elems)
00519     register int i;
00520     register double sum = 0;
00521     for(i=0; i<nrows; i++)              // Scan a col to compute the sum
00522       sum += fabs(*ep++);
00523     norm = max(norm,sum);
00524   }
00525   assert( ep == elements + nelems );
00526 
00527   return norm;
00528 }
00529 
00530 
00531                                 // Square of the Euclidian norm
00532                                 // SUM{ m(i,j)^2 }
00533 double Matrix::e2_norm(void) const
00534 {
00535   return of_every(*this).sum_squares();
00536 }
00537 
00538                                 // Square of the Euclidian norm of the
00539                                 // difference between two matrices
00540 double e2_norm(const Matrix& m1, const Matrix& m2)
00541 {
00542   return of_every(m1).sum_squares(of_every(m2));
00543 }
00544 
00545 /*
00546  *------------------------------------------------------------------------
00547  *                      Some service operations
00548  */
00549 
00550 ostream& operator << (ostream& os, const DimSpec& dimspec)
00551 {
00552   return os << dimspec.row_lwb << ':' << dimspec.nrows+dimspec.row_lwb-1 << 'x'
00553             << dimspec.col_lwb << ':' << dimspec.ncols+dimspec.col_lwb-1;
00554 }
00555 
00556 ostream& operator << (ostream& os, const AREALMark& mark)
00557 {
00558   return (bool)mark ?
00559         os << "stream mark at " << mark.offset
00560       : os << "invalid stream mark";
00561 }
00562 
00563 ostream& operator << (ostream& os, const rowcol& rc)
00564 {
00565   return os << '(' << rc.row << ',' << rc.col << ')';
00566 }
00567 
00568 ostream& operator << (ostream& os, const IRange range)
00569 {
00570   os << '[';
00571   if( range.lwb == -IRange::INF )
00572     os << "-INF";
00573   else if( range.lwb == IRange::INF )
00574     os << "INF";
00575   else
00576     os << range.lwb;
00577   os << ',';
00578   if( range.upb == -IRange::INF )
00579     os << "-INF";
00580   else if( range.upb == IRange::INF )
00581     os << "INF";
00582   else
00583     os << range.upb;
00584   return os << ']';
00585 }
00586 
00587 ostream& operator << (ostream& os, const RWWatchDog& wd)
00588 {
00589   if( wd.q_exclusive() )
00590     return os << "is being exclusively held";
00591   else if ( !wd.q_engaged() )
00592     return os << "is not engaged";
00593   else
00594     return os << "is shared among " << wd.ref_count << " views";
00595 }
00596 
00597 volatile void RWWatchDog::access_violation(const char reason [])
00598 {
00599   cerr << "RWWatchDog growls at " << reason
00600        << "\nwhile the guarded object " << *this << endl;
00601   _error("The program is terminated due to a reason above");
00602 }
00603 
00604 void Matrix::info(void) const   // Print some information about the matrix
00605 {
00606   is_valid();
00607   cerr << "\nMatrix " << static_cast<const DimSpec&>(*this) << ' '
00608        << name << ' ' << ref_counter << endl;
00609 }
00610 
00611                         // The matrix in question is not valid
00612                         // ==> crash the program
00613 volatile void Matrix::invalid_matrix(void) const
00614 {
00615   _error("An operation attempted on an invalid matrix reference: 0x%0x\n",
00616         (long)this);
00617 }
00618 
00619                         // Dump the current status of the stream
00620 ostream& AREALBlockStreamIn::dump(ostream& os) const
00621 {
00622   return
00623   os << "\nAREALBlockStreamIn: first_el_p " << hex << ((long int)first_el_p) << dec
00624       << ", last_el_p +" << (last_el_p - first_el_p)
00625       << "\ncurrent pos " << (curr_el_p - first_el_p)
00626       << ", last_col_el_p " << (last_col_el_p - first_el_p)
00627       << ", col_size " << col_size << ", eoc_jump " << eoc_jump << endl;
00628 }
00629 
00630 
00631                                 // Print the Matrix as a table of elements
00632                                 // (zeros are printed as dots)
00633 void Matrix::print(const char title []) const
00634 {
00635   is_valid();
00636   cerr << title << " follows " << endl;
00637   info();
00638 
00639   const int cols_per_sheet = 6;
00640   ConstMatrixDA accessor(*this);
00641   
00642   for(register int sheet_counter=1;
00643       sheet_counter<=q_ncols(); sheet_counter += cols_per_sheet)
00644   {
00645     message("\n\n     |");
00646     for(register int j=q_col_lwb() + sheet_counter - 1;
00647         j<q_col_lwb() - 1 + sheet_counter+cols_per_sheet && j<=q_col_upb(); j++)
00648       message("   %6d  |",j);
00649     message("\n%s\n",_Minuses);
00650     for(register int i=q_row_lwb(); i<=q_row_upb(); i++)
00651     {
00652       message("%4d |",i);
00653       for(register int j=q_col_lwb() + sheet_counter - 1;
00654           j<q_col_lwb() - 1 + sheet_counter+cols_per_sheet && 
00655           j<=q_col_upb(); j++)
00656         message("%11.4g  ",accessor(i,j));
00657       message("\n");
00658     }
00659   }
00660   message("Done\n");
00661 }
00662 
00663 void compare(                   // Compare the two Matrices
00664         const Matrix& matrix1,  // and print out the result of comparison
00665         const Matrix& matrix2,
00666         const char title [] )
00667 {
00668   are_compatible(matrix1,matrix2);
00669 
00670   cerr << "\n\nComparison of two Matrices:\n\t" << title;
00671   matrix1.info();
00672   matrix2.info();
00673 
00674   double norm1 = 0, norm2 = 0;          // Norm of the Matrices
00675   double ndiff = 0;                     // Norm of the difference
00676   REAL difmax = -1;
00677 
00678   AREALMark max_mark;           // For the elements that differ most
00679   LAStreamIn m1_stream(matrix1), m2_stream(matrix2);
00680   while( !m1_stream.eof() )
00681   {
00682     const REAL mv1 = m1_stream.get();
00683     const REAL mv2 = m2_stream.get();
00684     const REAL diff = fabs(mv1-mv2);
00685 
00686     if( diff > difmax )
00687       difmax = diff, max_mark = m1_stream.tell_prev();
00688     norm1 += fabs(mv1);
00689     norm2 += fabs(mv2);
00690     ndiff += diff;
00691   }
00692   assert( m2_stream.eof() );
00693   
00694   cerr << "\nMaximal discrepancy    \t\t" << difmax;
00695   cerr <<"\n   occured at the point\t\t" << m1_stream.get_pos(max_mark);
00696   const REAL mv1 = m1_stream.seek(max_mark).peek();
00697   const REAL mv2 = m2_stream.seek(max_mark).peek();
00698   cerr << "\n Matrix 1 element is    \t\t" << mv1;
00699   cerr << "\n Matrix 2 element is    \t\t" << mv2;
00700   cerr << "\n Absolute error v2[i]-v1[i]\t\t" << (mv2-mv1);
00701   cerr << "\n Relative error\t\t\t\t" << (mv2-mv1)/max((double)fabs(mv2+mv1)/2,1e-7)
00702        << endl;
00703 
00704   cerr << "\n||Matrix 1||   \t\t\t" << norm1;
00705   cerr << "\n||Matrix 2||   \t\t\t" << norm2;
00706   cerr << "\n||Matrix1-Matrix2||\t\t\t\t" << ndiff;
00707   cerr << "\n||Matrix1-Matrix2||/sqrt(||Matrix1|| ||Matrix2||)\t" 
00708        << ndiff/max( sqrt(norm1*norm2), 1e-7 ) << endl;
00709 
00710 }
00711 
00712 /*
00713  *------------------------------------------------------------------------
00714  *                      Service validation functions
00715  */
00716 
00717 //SB: had to take this out of verify_element_value
00718 struct MaxDev : public ConstElementAction
00719   {
00720     int imax, jmax;
00721     double max_dev;
00722     const REAL expected_val;
00723     REAL found_val;
00724     MaxDev(const REAL val) : imax(-1), jmax(-1), max_dev(0), 
00725         expected_val(val), found_val(0) {}
00726     void operation(const REAL element, const int i, const int j)
00727         { const double dev = fabs(element - expected_val);
00728           if( dev >= max_dev )
00729             max_dev = dev, found_val = element, imax =i, jmax = j; }
00730   };
00731 
00732 void verify_element_value(const Matrix& m,const REAL val)
00733 {
00734 
00735   MaxDev maxdev(val);
00736   m.apply(maxdev);
00737   
00738   if( maxdev.max_dev == 0 )
00739     return;
00740   else if( maxdev.max_dev < 1e-5 )
00741     message("Element (%d,%d) with value %g differs the most from what\n"
00742             "was expected, %g, though the deviation %g is small\n",
00743             maxdev.imax,maxdev.jmax,maxdev.found_val,val,maxdev.max_dev);
00744   else
00745     _error("A significant difference from the expected value %g\n"
00746            "encountered for element (%d,%d) with value %g",
00747            val,maxdev.imax,maxdev.jmax,maxdev.found_val);
00748 
00749 }
00750 
00751 
00752 void verify_matrix_identity(const Matrix& m1, const Matrix& m2)
00753 {
00754   are_compatible(m1,m2);
00755   
00756   AREALMark max_mark;           // For the elements that differ most
00757   register double max_dev = 0;
00758   LAStreamIn m1_stream(m1), m2_stream(m2);
00759   while( !m1_stream.eof() )
00760   {
00761     register const double dev = fabs(m1_stream.get() - m2_stream.get());
00762     if( dev >= max_dev )
00763       max_dev = dev, max_mark = m1_stream.tell_prev();
00764   }
00765 
00766   if( max_dev == 0 )
00767     return;
00768   cerr << "Two " << m1_stream.get_pos(max_mark)
00769        << " elements of matrices with values "
00770        << m1_stream.seek(max_mark).peek() << " and "
00771        << m2_stream.seek(max_mark).peek() << endl;
00772   if( max_dev < 1e-5 )
00773     cerr << "differ the most, although the deviation "
00774          << max_dev << " is small" << endl;
00775   else
00776     _error("with a significant difference %g\n",max_dev);
00777 }
00778 }

Generated on Wed Dec 15 21:20:29 2004 for vuVolume by  doxygen 1.3.9.1