00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
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
00023
00024 namespace linalg
00025 {
00026 using namespace linalg;
00027
00028
00029
00030
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)
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
00060 void Matrix::set_name(const char * new_name)
00061 {
00062 if( name != 0 && name[0] != '\0' )
00063 delete const_cast<char*>(name);
00064
00065 if( new_name == 0 || new_name[0] == '\0' )
00066 name = "";
00067 else
00068 name = new char[strlen(new_name)+1],
00069 strcpy(const_cast<char *>(name), new_name);
00070 }
00071
00072
00073
00074
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
00095
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
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
00128
00129 REAL * const * MatrixDABase::build_index(const Matrix& m)
00130 {
00131 m.is_valid();
00132 if( m.ncols == 1 )
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
00152
00153
00154
00155
00156
00157
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
00171
00172
00173
00174
00175
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
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
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++)
00209 *cp++ = norm_factor;
00210
00211
00212
00213
00214
00215
00216
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
00245
00246
00247
00248
00249
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
00268
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
00295
00296
00297
00298 void ElementWise::abs(void)
00299 {
00300 for(register REAL * ep=start_ptr; ep < end_ptr; ep++)
00301 *ep = fabs(*ep);
00302 }
00303
00304
00305 void ElementWise::sqr(void)
00306 {
00307 for(register REAL * ep=start_ptr; ep < end_ptr; ep++)
00308 *ep = *ep * *ep;
00309 }
00310
00311
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
00323
00324
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
00353
00354
00355
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
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
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
00408
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
00425
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
00481
00482
00483
00484
00485
00486
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)
00494 {
00495 register int j;
00496 register double sum = 0;
00497 for(j=0; j<ncols; j++,ep+=nrows)
00498 sum += fabs(*ep);
00499 ep -= nelems - 1;
00500 norm = max(norm,sum);
00501 }
00502 assert( ep == elements + nrows );
00503
00504 return norm;
00505 }
00506
00507
00508
00509
00510
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)
00518 {
00519 register int i;
00520 register double sum = 0;
00521 for(i=0; i<nrows; i++)
00522 sum += fabs(*ep++);
00523 norm = max(norm,sum);
00524 }
00525 assert( ep == elements + nelems );
00526
00527 return norm;
00528 }
00529
00530
00531
00532
00533 double Matrix::e2_norm(void) const
00534 {
00535 return of_every(*this).sum_squares();
00536 }
00537
00538
00539
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
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
00605 {
00606 is_valid();
00607 cerr << "\nMatrix " << static_cast<const DimSpec&>(*this) << ' '
00608 << name << ' ' << ref_counter << endl;
00609 }
00610
00611
00612
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
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
00632
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(
00664 const Matrix& matrix1,
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;
00675 double ndiff = 0;
00676 REAL difmax = -1;
00677
00678 AREALMark max_mark;
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
00715
00716
00717
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;
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 }