00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 #ifndef __GNUC__
00060 #pragma once
00061 #endif
00062 #ifndef _LinAlg_h
00063 #define _LinAlg_h 1
00064
00065 #if defined(__GNUC__)
00066 #pragma interface
00067 #define __unused(X) X __attribute__((unused))
00068
00069 #include <iostream.h>
00070 using namespace std;
00071
00072 #else
00073 #include <iostream.h>
00074 using namespace std;
00075 #define __unused(X) X
00076 #endif
00077
00078 #include "myenv.h"
00079 #include "std.h"
00080 #if defined(__GNUC__)
00081 #include <values.h>
00082 #ifndef MAXINT
00083 #include <limits.h>
00084 #define MAXINT INT_MAX
00085 #endif
00086 #else
00087 #include <limits.h>
00088 #define MAXINT INT_MAX
00089 #endif
00090 #include <math.h>
00091 #include "builtin.h"
00092 #include "minmax.h"
00093 #if defined(index)
00094 #undef index
00095 #endif
00096
00097
00098
00099
00100
00101
00102 typedef float REAL;
00103
00104 namespace linalg
00105 {
00106 using namespace linalg;
00107
00108
00109
00110 class DimSpec
00111 {
00112 protected:
00113 int nrows;
00114 int ncols;
00115 int row_lwb;
00116 int col_lwb;
00117
00118 public:
00119
00120 DimSpec(const int _nrows, const int _ncols)
00121 : nrows(_nrows), ncols(_ncols), row_lwb(1), col_lwb(1)
00122 { assert(nrows>0 && ncols>0); }
00123 DimSpec(const int _row_lwb, const int _row_upb,
00124 const int _col_lwb, const int _col_upb)
00125 : nrows(_row_upb-_row_lwb+1),
00126 ncols(_col_upb-_col_lwb+1),
00127 row_lwb(_row_lwb), col_lwb(_col_lwb)
00128 { assert(nrows>0 && ncols>0); }
00129
00130
00131 int q_row_lwb(void) const { return row_lwb; }
00132 int q_row_upb(void) const { return nrows+row_lwb-1; }
00133 int q_nrows(void) const { return nrows; }
00134 int q_col_lwb(void) const { return col_lwb; }
00135 int q_col_upb(void) const { return ncols+col_lwb-1; }
00136 int q_ncols(void) const { return ncols; }
00137
00138 bool operator == (const DimSpec& another_dim_spec) const
00139 { return nrows == another_dim_spec.nrows &&
00140 ncols == another_dim_spec.ncols &&
00141 row_lwb == another_dim_spec.row_lwb &&
00142 col_lwb == another_dim_spec.col_lwb; }
00143
00144 friend ostream& operator << (ostream& os, const DimSpec& dimspec);
00145 };
00146
00147
00148
00149 struct rowcol
00150 {
00151 int row, col;
00152 rowcol(const int _row, const int _col) : row(_row), col(_col) {}
00153 bool operator == (const rowcol another) const
00154 { return row == another.row && col == another.col; }
00155 friend ostream& operator << (ostream& os, const rowcol& rc);
00156 };
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166 struct IRange
00167 {
00168 const int lwb, upb;
00169 enum { INF=MAXINT };
00170 IRange(const int _lwb, const int _upb) : lwb(_lwb), upb(_upb) {}
00171 static IRange from(const int lwb) { return IRange(lwb,INF); }
00172 static IRange through(const int upb) { return IRange(-INF,upb); }
00173 bool q_empty(void) const { return lwb > upb; }
00174 bool q_proper(void) const
00175
00176 { return fabs((double)lwb) != INF && fabs((double)upb) != INF && lwb <= upb; }
00177 friend ostream& operator << (ostream& os, const IRange range);
00178 int no_less_than(const int strict_lwb) const
00179 { if( lwb == -IRange::INF ) return strict_lwb;
00180 assert(lwb >= strict_lwb); return lwb; }
00181 int no_more_than(const int strict_upb) const
00182 { if( upb == IRange::INF ) return strict_upb;
00183 assert(upb <= strict_upb); return upb; }
00184 };
00185
00186
00187 struct IRangeStride
00188 {
00189 const int lwb, upb, stride;
00190 IRangeStride(const int _lwb, const int _upb, const int _stride)
00191 : lwb(_lwb), upb(_upb), stride(_stride) {}
00192 IRangeStride(const IRange range)
00193 : lwb(range.lwb), upb(range.upb), stride(1) {}
00194 bool q_empty(void) const { return lwb > upb; }
00195 bool q_proper(void) const
00196
00197 { return fabs((double)lwb) != IRange::INF && fabs((double)upb) != IRange::INF
00198 && lwb <= upb && stride > 0; }
00199 };
00200
00201
00202
00203 class DimSpecSubranged : public DimSpec
00204 {
00205 protected:
00206 int min_offset;
00207 int max_offset;
00208 const int original_nrows;
00209
00210 public:
00211 DimSpecSubranged(const DimSpec& orig, const IRange row_range,
00212 const IRange col_range)
00213 : DimSpec(row_range.no_less_than(orig.q_row_lwb()),
00214 row_range.no_more_than(orig.q_row_upb()),
00215 col_range.no_less_than(orig.q_col_lwb()),
00216 col_range.no_more_than(orig.q_col_upb())),
00217 original_nrows(orig.q_nrows())
00218 { min_offset = (col_lwb - orig.q_col_lwb())*original_nrows +
00219 (row_lwb - orig.q_row_lwb());
00220 max_offset = min_offset + (ncols-1)*original_nrows + nrows - 1;
00221 }
00222
00223
00224
00225 int q_min_offset(void) const { return min_offset; }
00226 int q_max_offset(void) const { return max_offset; }
00227 int q_row_diff(void) const { return original_nrows - nrows; }
00228 };
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240 class ConstElementAction
00241 {
00242 public:
00243 virtual void operation(const REAL element, const int i, const int j) = 0;
00244 };
00245
00246
00247
00248
00249 class ElementAction
00250 {
00251 public:
00252 virtual void operation(REAL& element, const int i, const int j) = 0;
00253
00254
00255
00256
00257
00258 };
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272 class ConstElementAction2
00273 {
00274 public:
00275 virtual void operation(const REAL element1, const REAL element2,
00276 const int i, const int j) = 0;
00277 };
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287 class LazyMatrix : public DimSpec
00288 {
00289 friend class Matrix;
00290 virtual void fill_in(Matrix& m) const = 0;
00291
00292
00293
00294 LazyMatrix& operator = (const LazyMatrix&);
00295 public:
00296 LazyMatrix(const int nrows=1, const int ncols=1)
00297 : DimSpec(nrows,ncols) {}
00298 LazyMatrix(const int row_lwb, const int row_upb,
00299 const int col_lwb, const int col_upb)
00300 : DimSpec(row_lwb,row_upb,col_lwb,col_upb) {}
00301 LazyMatrix(const DimSpec& dims) : DimSpec(dims) {}
00302 };
00303
00304
00305
00306 class MinMax
00307 {
00308 REAL min_val, max_val;
00309 public:
00310 MinMax(REAL _min_val, REAL _max_val)
00311 : min_val(_min_val), max_val(_max_val) {}
00312 explicit MinMax(REAL val)
00313 : min_val(val), max_val(val) {}
00314 MinMax& operator << (REAL val)
00315 { if( val > max_val ) max_val = val;
00316 else if(val < min_val) min_val = val; return *this; }
00317 REAL min(void) const { return min_val; }
00318 REAL max(void) const { return max_val; }
00319 double ratio(void) const { return max_val/min_val; }
00320 friend ostream& operator << (ostream& os, const MinMax& mx);
00321 };
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347 class RWWatchDog
00348 {
00349 int ref_count;
00350
00351 RWWatchDog(const RWWatchDog&);
00352 void operator = (const RWWatchDog&);
00353
00354 volatile void access_violation(const char reason []);
00355
00356 public:
00357 bool q_engaged(void) const { return ref_count != 0; }
00358 bool q_shared(void) const { return ref_count > 0; }
00359 bool q_exclusive(void) const { return ref_count == -1; }
00360
00361 void get_exclusive(void)
00362 { if( q_engaged() ) access_violation("get_exclusive"); ref_count = -1; }
00363 void release_exclusive(void)
00364 { if(!q_exclusive()) access_violation("release_exclusive");
00365 ref_count = 0; }
00366
00367 void get_shared(void)
00368 { if( q_exclusive() ) access_violation("get_shared"); ref_count++; }
00369 void release_shared(void)
00370 { if( !q_shared() ) access_violation("release_shared"); ref_count--; }
00371
00372 RWWatchDog(void) : ref_count(0) {}
00373 ~RWWatchDog(void)
00374 { if( ref_count!=0 ) access_violation("destroying");
00375 ref_count = -1; }
00376
00377
00378 friend ostream& operator << (ostream& os, const RWWatchDog& wd);
00379 };
00380
00381 #if 0
00382 class MinMaxLoc : public MinMax
00383 {
00384 };
00385 #endif
00386
00387
00388
00389
00390
00391
00392
00393 class Matrix;
00394 class MatrixColumn;
00395 class ConstMatrixColumn;
00396 class MatrixRow;
00397 class MatrixDiag;
00398
00399
00400
00401
00402
00403 class haar_matrix : public LazyMatrix
00404 {
00405 public:
00406 void fill_in(Matrix& m) const;
00407 haar_matrix(const int n, const int no_cols=0);
00408 };
00409
00410
00411 class LazyTransposedMatrix : public LazyMatrix
00412 {
00413 const Matrix& proto;
00414 void fill_in(Matrix& m) const;
00415 inline LazyTransposedMatrix(const Matrix & _proto);
00416 public:
00417 friend inline const LazyTransposedMatrix transposed(const Matrix & proto);
00418 };
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432 class Matrix : public DimSpec
00433 {
00434 public:
00435 friend class ElementWiseConst;
00436 friend class ElementWiseStrideConst;
00437 friend class LAStreamIn;
00438 friend class LAStreamOut;
00439 friend class LAStrideStreamIn;
00440 friend class LAStrideStreamOut;
00441 friend class LABlockStreamIn;
00442 friend class LABlockStreamOut;
00443 friend class MatrixDABase;
00444 class ConstReference;
00445 class Reference;
00446 friend class ConstReference;
00447
00448 friend class Vector;
00449 friend class ConstMatrixColumn;
00450 friend class ConstMatrixRow;
00451 friend class ConstMatrixDiag;
00452
00453 private:
00454 int valid_code;
00455 enum { MATRIX_val_code = 5757 };
00456 RWWatchDog ref_counter;
00457
00458 protected:
00459 const char * name;
00460 int nelems;
00461 REAL * elements;
00462
00463 void allocate(void);
00464
00465 friend void haar_matrix::fill_in(Matrix& m) const;
00466 friend void LazyTransposedMatrix::fill_in(Matrix& m) const;
00467
00468 public:
00469
00470
00471
00472 Matrix(const int nrows, const int ncols);
00473 Matrix(const int row_lwb, const int row_upb,
00474 const int col_lwb, const int col_upb);
00475 Matrix(const DimSpec& dimension_specs);
00476 inline Matrix(const Matrix& another);
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490 Matrix(const LazyMatrix& lazy_constructor);
00491 explicit Matrix(const char * file_name);
00492
00493
00494 ~Matrix(void);
00495
00496 void set_name(const char * name);
00497
00498
00499
00500 void resize_to(const int nrows, const int ncols);
00501 void resize_to(const int row_lwb, const int row_upb,
00502 const int col_lwb, const int col_upb);
00503 void resize_to(const DimSpec& dimension_specs);
00504
00505
00506 void is_valid(void) const
00507 { if( valid_code != MATRIX_val_code ) invalid_matrix(); }
00508
00509
00510 int q_no_elems(void) const { return nelems; }
00511 const char * q_name(void) const { return name; }
00512
00513 class ConstReference
00514 {
00515 friend class Reference;
00516 Matrix& m;
00517 ConstReference(const ConstReference&);
00518 void operator = (const ConstReference&);
00519 public:
00520 ConstReference(const Matrix& _m)
00521 : m(const_cast<Matrix&>(_m))
00522 { m.is_valid(); m.ref_counter.get_shared(); }
00523 ~ConstReference(void) { m.is_valid(); m.ref_counter.release_shared(); }
00524 operator const Matrix& (void) const { return m; }
00525 const Matrix& ref(void) const { return m; }
00526 };
00527
00528 class Reference : public ConstReference
00529 {
00530 Reference(const Reference&);
00531 void operator = (const Reference&);
00532 public:
00533 Reference(Matrix& _m) : ConstReference(_m) {}
00534 operator Matrix& (void) const { return m; }
00535 Matrix& ref(void) const { return m; }
00536 };
00537
00538
00539 Matrix& apply(ElementAction& action);
00540 const Matrix& apply(ConstElementAction& action) const;
00541
00542
00543
00544
00545
00546 inline Matrix& operator = (const REAL val);
00547
00548
00549
00550
00551
00552
00553
00554 static double dummy_determinant_ref;
00555 Matrix& invert(double &determ_ref = dummy_determinant_ref);
00556
00557
00558 inline Matrix& operator = (const Matrix& source);
00559 inline Matrix& operator = (const ConstReference& ref)
00560 { return *this = ref.ref(); }
00561 Matrix& operator = (const LazyMatrix& source);
00562 Matrix& clear(void);
00563
00564
00565
00566 bool operator == (const Matrix& im2) const;
00567 friend inline void are_compatible(const Matrix& im1, const Matrix& im2);
00568
00569
00570
00571
00572 Matrix& operator *= (const Matrix& source);
00573
00574 inline Matrix& operator *= (const ConstReference& ref)
00575 { return *this *= ref.ref(); }
00576
00577
00578 void mult(const Matrix& A, const Matrix& B);
00579
00580
00581 double row_norm(void) const;
00582 double norm_inf(void) const
00583 { return row_norm(); }
00584 double col_norm(void) const;
00585 double norm_1(void) const
00586 { return col_norm(); }
00587 double e2_norm(void) const;
00588
00589
00590 friend double e2_norm(const Matrix& m1, const Matrix& m2);
00591
00592
00593 double determinant(void) const;
00594
00595 double asymmetry_index(void) const;
00596
00597
00598
00599
00600 Matrix& unit_matrix(void);
00601 Matrix& hilbert_matrix(void);
00602
00603
00604
00605
00606
00607 void write(const char * file_name,const char * title = "") const;
00608 void info(void) const;
00609 void print(const char title []) const;
00610 volatile void invalid_matrix(void) const;
00611
00612 };
00613
00614 void compare(const Matrix& im1, const Matrix& im2, const char title[]);
00615
00616
00617
00618
00619
00620
00621 class Vector : public Matrix
00622 {
00623 public:
00624 Vector(const int n);
00625
00626
00627 Vector(const int lwb, const int upb);
00628
00629
00630
00631 Vector(const int lwb, const int upb,
00632 double iv1, ...
00633 );
00634
00635
00636 Vector(const LazyMatrix& lazy_constructor);
00637
00638
00639
00640 void resize_to(const int n);
00641 void resize_to(const int lwb, const int upb);
00642 void resize_to(const Vector& v);
00643
00644
00645
00646
00647
00648 inline REAL& operator () (const int index);
00649 inline REAL operator () (const int index) const
00650 { return (const_cast<Vector&>(*this))(index); }
00651
00652
00653
00654
00655
00656 int q_lwb(void) const { return q_row_lwb(); }
00657 int q_upb(void) const { return q_row_upb(); }
00658
00659
00660 friend double operator * (const Vector& v1, const Vector& v2);
00661
00662
00663
00664
00665
00666 Vector& operator *= (const Matrix& A);
00667
00668
00669 inline double norm_1(void) const;
00670 inline double norm_2_sqr(void) const;
00671 inline double norm_inf(void) const;
00672
00673 Vector& operator = (const Vector& v)
00674 { Matrix::operator =(v); return *this; }
00675 Vector& operator = (const LazyMatrix& source);
00676
00677
00678
00679
00680
00681 inline Vector& operator = (const ElementWiseConst& stream);
00682 inline Vector& operator = (const ElementWiseStrideConst& stream);
00683 inline Vector& operator = (const REAL val)
00684 { Matrix::operator =(val); return *this; }
00685 };
00686
00687
00688
00689
00690 void verify_element_value(const Matrix& m,const REAL val);
00691 void verify_matrix_identity(const Matrix& m1, const Matrix& m2);
00692 inline void verify_element_value(const Matrix::ConstReference& m,const REAL val)
00693 { verify_element_value(m.ref(),val); }
00694 inline void verify_matrix_identity(const Matrix::ConstReference& m1,
00695 const Matrix& m2)
00696 { verify_matrix_identity(m1.ref(),m2); }
00697 inline void verify_matrix_identity(const Matrix& m1,
00698 const Matrix::ConstReference& m2)
00699 { verify_matrix_identity(m1,m2.ref()); }
00700 inline void verify_matrix_identity(const Matrix::ConstReference& m1,
00701 const Matrix::ConstReference& m2)
00702 { verify_matrix_identity(m1.ref(),m2.ref()); }
00703
00704
00705
00706 Matrix& operator *= (Matrix& m, const ConstMatrixDiag& diag);
00707
00708
00709 class LazyMatrixProduct : public linalg::LazyMatrix
00710 {
00711 private:
00712 const Matrix& A;
00713 const Matrix& B;
00714 void fill_in(Matrix& m) const { m.mult(A,B); };
00715
00716 LazyMatrixProduct(const Matrix & __A,
00717 const Matrix & __B)
00718 : LazyMatrix(__A.q_row_lwb(),__A.q_row_upb(),
00719 __B.q_col_lwb(),__B.q_col_upb()),
00720 A(__A), B(__B) {};
00721
00722
00723 public:
00724 friend const LazyMatrixProduct operator * (const Matrix& A, const Matrix& B)
00725 { return LazyMatrixProduct(A,B); };
00726
00727 };
00728
00729
00730 class LazyZeroMatrix : public LazyMatrix
00731 {
00732 void fill_in(Matrix& m) const { __unused(Matrix& m1) = m; }
00733 LazyZeroMatrix(const Matrix& proto)
00734 : LazyMatrix(static_cast<const DimSpec&>(proto)) {}
00735 public:
00736 friend const LazyZeroMatrix zero(const Matrix& proto)
00737 { proto.is_valid(); return proto; }
00738 };
00739
00740
00741 class LazyUnitMatrix : public LazyMatrix
00742 {
00743 void fill_in(Matrix& m) const { m.unit_matrix(); }
00744 LazyUnitMatrix(const Matrix& proto)
00745 : LazyMatrix(static_cast<const DimSpec&>(proto)) {}
00746 public:
00747 friend const LazyUnitMatrix unit(const Matrix& proto)
00748 { proto.is_valid(); return proto; }
00749 };
00750
00751
00752 class LazyInverseMatrix : public LazyMatrix
00753 {
00754 const Matrix& proto;
00755 void fill_in(Matrix& m) const { m = proto; m.invert(); }
00756 LazyInverseMatrix(const Matrix& _proto)
00757 : LazyMatrix(static_cast<const DimSpec&>(_proto)),
00758 proto(_proto) {}
00759 public:
00760 friend const LazyInverseMatrix inverse(const Matrix& proto)
00761 { proto.is_valid(); return proto; }
00762 };
00763
00764
00765
00766 inline LazyTransposedMatrix::LazyTransposedMatrix(const Matrix & _proto)
00767 : LazyMatrix(_proto.q_col_lwb(),_proto.q_col_upb(),
00768 _proto.q_row_lwb(),_proto.q_row_upb()),
00769 proto(_proto) {}
00770 inline const LazyTransposedMatrix transposed(const Matrix & proto)
00771 { proto.is_valid(); return proto; }
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811 class ElementWiseConstAction
00812 {
00813
00814
00815 ElementWiseConstAction(const ElementWiseConstAction&);
00816 ElementWiseConstAction& operator= (const ElementWiseConstAction&);
00817 public:
00818 ElementWiseConstAction(void) {}
00819 virtual void operator () (const REAL element) = 0;
00820 };
00821
00822
00823
00824
00825 class ElementWiseAction
00826 {
00827
00828
00829 ElementWiseAction(const ElementWiseAction&);
00830 ElementWiseAction& operator= (const ElementWiseAction&);
00831 public:
00832 ElementWiseAction(void) {}
00833 virtual void operator () (REAL& element) = 0;
00834 };
00835
00836 class ElementWiseConst
00837 {
00838 public:
00839 friend class ElementWise;
00840 friend class ElementWiseStrideConst;
00841 friend class ElementWiseStride;
00842
00843 REAL * const start_ptr;
00844
00845 REAL * const end_ptr;
00846
00847 void operator=(const ElementWiseConst&);
00848
00849 ElementWiseConst(const ElementWiseConst&);
00850
00851
00852
00853
00854
00855
00856
00857 ElementWiseConst(const Matrix& m)
00858 : start_ptr(const_cast<REAL*>(m.elements)),
00859 end_ptr(const_cast<REAL*>(m.elements+m.nelems))
00860 { m.is_valid(); }
00861
00862 inline ElementWiseConst(const ConstMatrixColumn& mc);
00863
00864
00865
00866
00867 struct ElementWiseConstRef
00868 {
00869 const ElementWiseConst& ref;
00870 ElementWiseConstRef(const ElementWiseConstRef&);
00871 void operator = (const ElementWiseConstRef&);
00872 ElementWiseConstRef(const ElementWiseConst& _ref) : ref(_ref) {}
00873 };
00874
00875 public:
00876
00877
00878
00879
00880
00881 friend inline ElementWiseConst of_every(const Matrix& m) { return m; }
00882 friend inline ElementWiseConst of_every(const Matrix::ConstReference& m)
00883 { return m.ref(); }
00884
00885 friend inline ElementWiseConst of_every(const ConstMatrixColumn& mc);
00886
00887
00888
00889 inline operator ElementWiseStrideConst (void);
00890
00891
00892
00893
00894
00895 bool operator == (const REAL val) const;
00896 bool operator != (const REAL val) const;
00897 bool operator < (const REAL val) const;
00898 bool operator <= (const REAL val) const;
00899 bool operator > (const REAL val) const;
00900 bool operator >= (const REAL val) const;
00901
00902
00903
00904
00905 bool operator == (const ElementWiseConst& another) const;
00906 bool operator != (const ElementWiseConst& another) const;
00907 bool operator < (const ElementWiseConst& another) const;
00908 bool operator <= (const ElementWiseConst& another) const;
00909 bool operator > (const ElementWiseConst& another) const;
00910 bool operator >= (const ElementWiseConst& another) const;
00911
00912 void sure_compatible_with(const ElementWiseConst& another) const
00913 { assure(end_ptr-start_ptr == another.end_ptr-another.start_ptr,
00914 "Two groups of elements to operate on have different sizes"); }
00915
00916
00917
00918 double sum(void) const;
00919 double sum_squares(void) const;
00920 double sum_abs(void) const;
00921 double max_abs(void) const;
00922
00923
00924
00925
00926
00927 double sum_squares(const ElementWiseConst& another) const;
00928 double sum_abs(const ElementWiseConst& another) const;
00929 double max_abs(const ElementWiseConst& another) const;
00930
00931
00932
00933
00934
00935 ElementWiseConstAction& apply(ElementWiseConstAction& functor) const;
00936 };
00937
00938
00939 class ElementWise : public ElementWiseConst
00940 {
00941 void operator=(const ElementWise&);
00942
00943 ElementWise(const ElementWise&);
00944
00945
00946
00947
00948
00949
00950
00951 ElementWise(Matrix& m) : ElementWiseConst(m) {}
00952 inline ElementWise(const MatrixColumn& mc);
00953 public:
00954
00955
00956
00957
00958
00959 friend inline ElementWise to_every(Matrix& m) { return m; }
00960 friend inline ElementWise to_every(const Matrix::Reference& m)
00961 { return m.ref(); }
00962
00963 friend inline ElementWise to_every(const MatrixColumn& mc);
00964
00965
00966
00967 inline operator ElementWiseStride (void);
00968
00969
00970
00971 inline ElementWiseStride with_stride (void);
00972
00973
00974
00975
00976 void operator = (const REAL val);
00977 void operator += (const double val);
00978 void operator -= (const double val);
00979 void operator *= (const double val);
00980
00981
00982 void abs(void);
00983 void sqr(void);
00984 void sqrt(void);
00985
00986 void operator = (const ElementWiseConst& another);
00987 void operator += (const ElementWiseConst& another);
00988 void operator -= (const ElementWiseConst& another);
00989 void operator *= (const ElementWiseConst& another);
00990 void operator /= (const ElementWiseConst& another);
00991
00992 ElementWiseAction& apply(ElementWiseAction& functor);
00993 };
00994
00995 inline Vector& Vector::operator = (const ElementWiseConst& stream)
00996 { to_every(*this) = stream; return *this; }
00997
00998
00999 inline double Vector::norm_1(void) const
01000 { return of_every(*this).sum_abs(); }
01001 inline double Vector::norm_2_sqr(void) const
01002 { return of_every(*this).sum_squares(); }
01003 inline double Vector::norm_inf(void) const
01004 { return of_every(*this).max_abs(); }
01005
01006
01007
01008
01009 class ElementWiseStrideConst
01010 {
01011 friend class ElementWiseStride;
01012 friend class ElementWiseConst;
01013
01014 REAL * const start_ptr;
01015
01016 REAL * const end_ptr;
01017 const int stride;
01018
01019 void operator=(const ElementWiseStrideConst&);
01020
01021 ElementWiseStrideConst(const ElementWiseStrideConst&);
01022
01023
01024
01025
01026
01027
01028
01029 ElementWiseStrideConst(const ElementWiseConst::ElementWiseConstRef& ewc)
01030 : start_ptr(ewc.ref.start_ptr),
01031 end_ptr(ewc.ref.end_ptr),
01032 stride(1)
01033 {}
01034
01035 inline ElementWiseStrideConst(const ConstMatrixRow& mr);
01036 inline ElementWiseStrideConst(const ConstMatrixDiag& md);
01037
01038 public:
01039
01040
01041
01042 friend inline ElementWiseStrideConst of_every(const ConstMatrixRow& mr);
01043 friend inline ElementWiseStrideConst of_every(const ConstMatrixDiag& md);
01044
01045
01046
01047
01048
01049 bool operator == (const REAL val) const;
01050 bool operator != (const REAL val) const;
01051 bool operator < (const REAL val) const;
01052 bool operator <= (const REAL val) const;
01053 bool operator > (const REAL val) const;
01054 bool operator >= (const REAL val) const;
01055
01056
01057
01058
01059 bool operator == (const ElementWiseStrideConst& another) const;
01060 bool operator != (const ElementWiseStrideConst& another) const;
01061 bool operator < (const ElementWiseStrideConst& another) const;
01062 bool operator <= (const ElementWiseStrideConst& another) const;
01063 bool operator > (const ElementWiseStrideConst& another) const;
01064 bool operator >= (const ElementWiseStrideConst& another) const;
01065
01066
01067
01068
01069
01070
01071
01072 double sum(void) const;
01073 double sum_squares(void) const;
01074 double sum_abs(void) const;
01075 double max_abs(void) const;
01076
01077
01078
01079
01080
01081 double sum_squares(const ElementWiseStrideConst& another) const;
01082 double sum_abs(const ElementWiseStrideConst& another) const;
01083 double max_abs(const ElementWiseStrideConst& another) const;
01084
01085
01086
01087
01088
01089 ElementWiseConstAction& apply(ElementWiseConstAction& functor) const;
01090 };
01091
01092 class ElementWiseStride : public ElementWiseStrideConst
01093 {
01094 friend class ElementWise;
01095 void operator=(const ElementWiseStride&);
01096
01097 ElementWiseStride(const ElementWiseStride&);
01098
01099
01100
01101
01102
01103
01104
01105 inline ElementWiseStride(const MatrixRow& mr);
01106 inline ElementWiseStride(const MatrixDiag& md);
01107 ElementWiseStride(const ElementWiseConst::ElementWiseConstRef& ewc)
01108 : ElementWiseStrideConst(ewc) {}
01109 public:
01110
01111
01112
01113
01114
01115 friend inline ElementWiseStride to_every(const MatrixRow& mr);
01116 friend inline ElementWiseStride to_every(const MatrixDiag& md);
01117
01118
01119
01120
01121 void operator = (const REAL val);
01122 void operator += (const double val);
01123 void operator -= (const double val);
01124 void operator *= (const double val);
01125
01126
01127 void abs(void);
01128 void sqr(void);
01129 void sqrt(void);
01130
01131 void operator = (const ElementWiseStrideConst& another);
01132 void operator += (const ElementWiseStrideConst& another);
01133 void operator -= (const ElementWiseStrideConst& another);
01134 void operator *= (const ElementWiseStrideConst& another);
01135 void operator /= (const ElementWiseStrideConst& another);
01136
01137 ElementWiseAction& apply(ElementWiseAction& functor);
01138 };
01139
01140 inline ElementWiseConst::operator ElementWiseStrideConst (void)
01141 { return ElementWiseConstRef(*this); }
01142
01143 inline ElementWise::operator ElementWiseStride (void)
01144 { return ElementWiseConstRef(*this); }
01145
01146 inline ElementWiseStride ElementWise::with_stride (void)
01147 { return ElementWiseConstRef(*this); }
01148
01149 inline Vector& Vector::operator = (const ElementWiseStrideConst& stream)
01150 { to_every(*this).with_stride() = stream; return *this; }
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161 class MatrixDABase : public DimSpec
01162 {
01163 REAL * const * const index;
01164 MatrixDABase(const MatrixDABase&);
01165 void operator = (const MatrixDABase&);
01166
01167 static REAL * const * build_index(const Matrix& m);
01168
01169 protected:
01170 MatrixDABase(const Matrix& m)
01171 : DimSpec(m),
01172 index( build_index(m) )
01173 {}
01174 inline const REAL& operator () (const int rown, const int coln) const;
01175 public:
01176 ~MatrixDABase(void);
01177 };
01178
01179 class ConstMatrixDA : public Matrix::ConstReference, public MatrixDABase
01180 {
01181 ConstMatrixDA(const ConstMatrixDA&);
01182 void operator = (const ConstMatrixDA&);
01183 public:
01184 ConstMatrixDA(const Matrix& m)
01185 : Matrix::ConstReference(m), MatrixDABase(m)
01186 {}
01187
01188 const REAL operator () (const int rown, const int coln) const
01189 { return MatrixDABase::operator() (rown,coln); }
01190 const REAL operator () (const rowcol rc) const
01191 { return operator () (rc.row,rc.col); }
01192
01193
01194
01195
01196
01197
01198 };
01199
01200 class MatrixDA : public Matrix::Reference, public MatrixDABase
01201 {
01202 MatrixDA(const MatrixDA&);
01203 void operator = (const MatrixDA&);
01204
01205 public:
01206 MatrixDA(Matrix& m)
01207 : Matrix::Reference(m), MatrixDABase(m)
01208 {}
01209
01210 REAL& operator () (const int rown, const int coln) const
01211 { return const_cast<REAL&>(MatrixDABase::operator()(rown,coln)); }
01212 REAL& operator () (const rowcol rc) const
01213 { return operator () (rc.row,rc.col); }
01214
01215
01216
01217
01218
01219 #if defined (__GNUC__)
01220 Matrix& get_container(void) { return Matrix::Reference::operator Matrix& (); }
01221 #else
01222 Matrix& get_container(void) { return *this; }
01223 #endif
01224 };
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234 class ConstMatrixColumn : protected Matrix::ConstReference,
01235 public DimSpec
01236 {
01237 friend class ElementWiseConst;
01238 friend class LAStreamIn;
01239 friend class LAStreamOut;
01240 const REAL * const col_ptr;
01241
01242 ConstMatrixColumn(const ConstMatrixColumn&);
01243 void operator = (const ConstMatrixColumn&);
01244
01245 protected:
01246
01247 inline const REAL& get_ref(const int index) const;
01248
01249 public:
01250 ConstMatrixColumn(const Matrix& m, const int col);
01251
01252 const REAL operator () (const int i) const
01253 { return get_ref(i); }
01254 };
01255
01256 inline ElementWiseConst::ElementWiseConst(const ConstMatrixColumn& mc)
01257 : start_ptr(const_cast<REAL*>(mc.col_ptr)),
01258 end_ptr(const_cast<REAL*>(mc.col_ptr)+mc.nrows)
01259 { }
01260
01261
01262 inline const REAL& ConstMatrixColumn::get_ref(const int index) const
01263 { register const int offset = index - row_lwb;
01264 if( offset >= nrows || offset < 0 )
01265 _error("MatrixColumn index %d is out of row boundaries [%d,%d]",
01266 index,q_row_lwb(),q_row_upb());
01267 return col_ptr[offset];
01268 }
01269
01270
01271 inline ElementWiseConst of_every(const ConstMatrixColumn& mc)
01272 { return mc; }
01273
01274 class MatrixColumn : public ConstMatrixColumn
01275 {
01276 friend class ElementWise;
01277 friend class LAStreamOut;
01278 void operator = (const MatrixColumn&);
01279 public:
01280 MatrixColumn(Matrix& m, const int col) :
01281 ConstMatrixColumn(m,col) {}
01282
01283 REAL& operator () (const int i)
01284 { return const_cast<REAL&>(get_ref(i)); }
01285 };
01286
01287 inline ElementWise::ElementWise(const MatrixColumn& mc) : ElementWiseConst(mc) {}
01288 inline ElementWise to_every(const MatrixColumn& mc) { return mc; }
01289
01290
01291 class ConstMatrixRow : protected Matrix::ConstReference,
01292 public DimSpec
01293 {
01294 friend class ElementWiseStrideConst;
01295 friend class LAStrideStreamIn;
01296 friend class LAStrideStreamOut;
01297 const REAL * const row_ptr;
01298
01299 const int stride;
01300
01301
01302
01303 const REAL * const end_ptr;
01304
01305 ConstMatrixRow(const ConstMatrixRow&);
01306 void operator = (const ConstMatrixRow&);
01307
01308 protected:
01309
01310 inline const REAL& get_ref(const int index) const;
01311
01312 public:
01313 ConstMatrixRow(const Matrix& m, const int row);
01314
01315 const REAL operator () (const int i) const
01316 { return get_ref(i); }
01317 };
01318
01319 inline
01320 ElementWiseStrideConst::ElementWiseStrideConst(const ConstMatrixRow& mr)
01321 : start_ptr(const_cast<REAL*>(mr.row_ptr)),
01322 end_ptr(const_cast<REAL*>(mr.end_ptr)),
01323 stride(mr.stride)
01324 { }
01325
01326
01327
01328 inline ElementWiseStrideConst of_every(const ConstMatrixRow& mr)
01329 { return mr; }
01330
01331
01332
01333 inline const REAL& ConstMatrixRow::get_ref(const int index) const
01334 { register const int offset = index - col_lwb;
01335 if( offset >= ncols || offset < 0 )
01336 _error("MatrixRow index %d is out of column boundaries [%d,%d]",
01337 index,q_col_lwb(),q_col_upb());
01338 return row_ptr[stride*offset];
01339 }
01340
01341
01342 class MatrixRow : public ConstMatrixRow
01343 {
01344 friend class ElementWiseStride;
01345 friend class LAStrideStreamOut;
01346 MatrixRow(const MatrixRow&);
01347 void operator = (const MatrixRow&);
01348 public:
01349 MatrixRow(Matrix& m, const int row) :
01350 ConstMatrixRow(m,row) {}
01351
01352 REAL& operator () (const int i)
01353 { return const_cast<REAL&>(get_ref(i)); }
01354 };
01355
01356 inline ElementWiseStride::ElementWiseStride(const MatrixRow& mr)
01357 : ElementWiseStrideConst(mr) {}
01358
01359 inline ElementWiseStride to_every(const MatrixRow& mr) { return mr; }
01360
01361
01362
01363
01364 class ConstMatrixDiag : protected Matrix::ConstReference,
01365 public DimSpec
01366 {
01367 friend class ElementWiseStrideConst;
01368 friend class LAStrideStreamIn;
01369 friend class LAStrideStreamOut;
01370 const REAL * const start_ptr;
01371 const int stride;
01372
01373
01374
01375 const REAL * const end_ptr;
01376
01377 ConstMatrixDiag(const ConstMatrixDiag&);
01378 void operator = (const ConstMatrixDiag&);
01379
01380 protected:
01381
01382
01383
01384
01385 inline const REAL& get_ref(const int index) const;
01386
01387 public:
01388 explicit ConstMatrixDiag(const Matrix& m);
01389
01390 const REAL operator () (const int i) const
01391 { return get_ref(i); }
01392 };
01393
01394 inline
01395 ElementWiseStrideConst::ElementWiseStrideConst(const ConstMatrixDiag& md)
01396 : start_ptr(const_cast<REAL*>(md.start_ptr)),
01397 end_ptr(const_cast<REAL*>(md.end_ptr)),
01398 stride(md.stride)
01399 { }
01400
01401
01402
01403 inline ElementWiseStrideConst of_every(const ConstMatrixDiag& md)
01404 { return md; }
01405
01406
01407
01408 inline const REAL& ConstMatrixDiag::get_ref(const int index) const
01409 { if( index > ncols || index < 1 )
01410 _error("MatrixDiag index %d is out of diag boundaries [%d,%d]",
01411 index,1,ncols);
01412 return start_ptr[stride*(index-1)];
01413 }
01414
01415
01416 class MatrixDiag : public ConstMatrixDiag
01417 {
01418 friend class ElementWiseStride;
01419 friend class LAStrideStreamOut;
01420 MatrixDiag(const MatrixDiag&);
01421 void operator = (const MatrixDiag&);
01422 public:
01423 explicit MatrixDiag(Matrix& m) : ConstMatrixDiag(m) {}
01424
01425 REAL& operator () (const int i)
01426 { return const_cast<REAL&>(get_ref(i)); }
01427 };
01428
01429 inline ElementWiseStride::ElementWiseStride(const MatrixDiag& md)
01430 : ElementWiseStrideConst(md) {}
01431
01432 inline ElementWiseStride to_every(const MatrixDiag& md) { return md; }
01433
01434
01435
01436
01437
01438
01439
01440 inline Matrix::Matrix(const int no_rows, const int no_cols)
01441 : DimSpec(no_rows,no_cols)
01442 {
01443 allocate();
01444 }
01445
01446 inline Matrix::Matrix(const int row_lwb, const int row_upb,
01447 const int col_lwb, const int col_upb)
01448 : DimSpec(row_lwb,row_upb,col_lwb,col_upb)
01449 {
01450 allocate();
01451 }
01452
01453 inline Matrix::Matrix(const DimSpec& dimension_specs)
01454 : DimSpec(dimension_specs)
01455 {
01456 allocate();
01457 }
01458
01459 inline Matrix::Matrix(const LazyMatrix& lazy_constructor)
01460 : DimSpec(lazy_constructor)
01461 {
01462 allocate();
01463 lazy_constructor.fill_in(*this);
01464 }
01465
01466
01467
01468
01469
01470 inline Matrix& Matrix::operator = (const LazyMatrix& lazy_constructor)
01471 {
01472 is_valid();
01473 if( !(static_cast<const DimSpec&>(lazy_constructor) == *this) )
01474 info(),
01475 _error("The matrix above is incompatible with the assigned "
01476 "Lazy matrix");
01477
01478 lazy_constructor.fill_in(*this);
01479 return *this;
01480 }
01481
01482
01483 inline Matrix::Matrix(const Matrix& another)
01484 : DimSpec(another)
01485 {
01486 another.is_valid();
01487 allocate();
01488 *this = another;
01489 }
01490
01491
01492 inline void Matrix::resize_to(const int row_lwb, const int row_upb,
01493 const int col_lwb, const int col_upb)
01494 {
01495 resize_to(DimSpec(row_lwb,row_upb,col_lwb,col_upb));
01496 }
01497
01498 inline const REAL& MatrixDABase::operator () (const int rown, const int coln) const
01499 {
01500 register int arown = rown-row_lwb;
01501 register int acoln = coln-col_lwb;
01502
01503 if( arown >= nrows || arown < 0 )
01504 _error("Row index %d is out of Matrix boundaries [%d,%d]",
01505 rown,row_lwb,nrows+row_lwb-1);
01506 if( acoln >= ncols || acoln < 0 )
01507 _error("Col index %d is out of Matrix boundaries [%d,%d]",
01508 coln,col_lwb,ncols+col_lwb-1);
01509
01510 return (index[acoln])[arown];
01511 }
01512
01513 inline Matrix& Matrix::clear(void)
01514 {
01515 is_valid();
01516 memset(elements,0,nelems*sizeof(REAL));
01517 return *this;
01518 }
01519
01520 inline void are_compatible(const Matrix& m1, const Matrix& m2)
01521 {
01522 m1.is_valid();
01523 m2.is_valid();
01524
01525 if( !(static_cast<const DimSpec&>(m1) == m2) )
01526 m1.info(), m2.info(), _error("The matrices above are incompatible");
01527 }
01528
01529
01530 inline Matrix& Matrix::operator = (const Matrix& source)
01531 {
01532 are_compatible(*this,source);
01533 memcpy(elements,source.elements,nelems*sizeof(REAL));
01534 return *this;
01535 }
01536
01537
01538
01539 inline ElementWiseConstAction&
01540 ElementWiseConst::apply(ElementWiseConstAction& functor) const
01541 {
01542 register const REAL * ep = start_ptr;
01543 while( ep < end_ptr )
01544 functor(*ep++);
01545 return functor;
01546 }
01547
01548 inline ElementWiseAction&
01549 ElementWise::apply(ElementWiseAction& functor)
01550 {
01551 register REAL * ep = start_ptr;
01552 while( ep < end_ptr )
01553 functor(*ep++);
01554 return functor;
01555 }
01556
01557 inline ElementWiseConstAction&
01558 ElementWiseStrideConst::apply(ElementWiseConstAction& functor) const
01559 {
01560 for(register const REAL * ep=start_ptr; ep < end_ptr; ep += stride )
01561 functor(*ep);
01562 return functor;
01563 }
01564
01565 inline ElementWiseAction&
01566 ElementWiseStride::apply(ElementWiseAction& functor)
01567 {
01568 for(register REAL * ep=start_ptr; ep < end_ptr; ep += stride )
01569 functor(*ep);
01570 return functor;
01571 }
01572
01573
01574
01575 inline Vector::Vector(const int n)
01576 : Matrix(n,1) {}
01577
01578
01579 inline Vector::Vector(const int lwb, const int upb) : Matrix(lwb,upb,1,1) {}
01580
01581
01582 inline Vector& Vector::operator = (const LazyMatrix& source)
01583 { Matrix::operator =(source); return *this; }
01584
01585
01586
01587
01588
01589
01590 inline void Vector::resize_to(const int n) { Vector::resize_to(1,n); }
01591
01592 inline void Vector::resize_to(const Vector& v)
01593 {
01594 Vector::resize_to(v.q_lwb(),v.q_upb());
01595 }
01596
01597
01598 inline REAL& Vector::operator () (const int ind)
01599 {
01600 is_valid();
01601 register const int aind = ind - row_lwb;
01602 if( aind >= nelems || aind < 0 )
01603 _error("Requested element %d is out of Vector boundaries [%d,%d]",
01604 ind,row_lwb,nrows+row_lwb-1);
01605
01606 return elements[aind];
01607 }
01608
01609
01610 inline Vector::Vector(const LazyMatrix& lazy_constructor)
01611 : Matrix(lazy_constructor)
01612 {
01613 assure(ncols == 1,
01614 "cannot make a vector from a promise of a full matrix");
01615 }
01616
01617
01618
01619 inline Matrix& Matrix::operator = (const REAL val)
01620 { to_every(*this) = val; return *this; }
01621
01622 inline Matrix& operator -= (Matrix& m, const double val)
01623 { to_every(m) -= val; return m; }
01624 inline Matrix& operator += (Matrix& m, const double val)
01625 { to_every(m) += val; return m; }
01626 inline Matrix& operator *= (Matrix& m, const double val)
01627 { to_every(m) *= val; return m; }
01628
01629 inline Matrix::Reference& operator -= (Matrix::Reference& m, const double val)
01630 { to_every(m) -= val; return m; }
01631 inline Matrix::Reference& operator += (Matrix::Reference& m, const double val)
01632 { to_every(m) += val; return m; }
01633 inline Matrix::Reference& operator *= (Matrix::Reference& m, const double val)
01634 { to_every(m) *= val; return m; }
01635
01636
01637 inline bool operator == (const Matrix& m, const REAL val)
01638 { return of_every(m) == val; }
01639 inline bool operator != (const Matrix& m, const REAL val)
01640 { return of_every(m) != val; }
01641 inline bool operator == (const Matrix::ConstReference& m, const REAL val)
01642 { return of_every(m) == val; }
01643 inline bool operator != (const Matrix::ConstReference& m, const REAL val)
01644 { return of_every(m) != val; }
01645 }
01646
01647
01648 #endif
01649