| 1 | /// \ingroup newmat
 | 
|---|
| 2 | ///@{
 | 
|---|
| 3 | 
 | 
|---|
| 4 | /// \file newmatrm.h
 | 
|---|
| 5 | /// Rectangular matrix operations.
 | 
|---|
| 6 | 
 | 
|---|
| 7 | // Copyright (C) 1991,2,3,4: R B Davies
 | 
|---|
| 8 | 
 | 
|---|
| 9 | #ifndef NEWMATRM_LIB
 | 
|---|
| 10 | #define NEWMATRM_LIB 0
 | 
|---|
| 11 | 
 | 
|---|
| 12 | #ifdef use_namespace
 | 
|---|
| 13 | namespace NEWMAT {
 | 
|---|
| 14 | #endif
 | 
|---|
| 15 | 
 | 
|---|
| 16 | 
 | 
|---|
| 17 | class RectMatrixCol;
 | 
|---|
| 18 | 
 | 
|---|
| 19 | /// Access rows and columns of a rectangular matrix.
 | 
|---|
| 20 | /// \internal
 | 
|---|
| 21 | class RectMatrixRowCol
 | 
|---|
| 22 | {
 | 
|---|
| 23 | protected:
 | 
|---|
| 24 | #ifdef use_namespace              // to make namespace work
 | 
|---|
| 25 | public:
 | 
|---|
| 26 | #endif
 | 
|---|
| 27 |    Real* store;                   // pointer to storage
 | 
|---|
| 28 |    int n;                         // number of elements
 | 
|---|
| 29 |    int spacing;                   // space between elements
 | 
|---|
| 30 |    int shift;                     // space between cols or rows
 | 
|---|
| 31 |    RectMatrixRowCol(Real* st, int nx, int sp, int sh)
 | 
|---|
| 32 |       : store(st), n(nx), spacing(sp), shift(sh) {}
 | 
|---|
| 33 |    void Reset(Real* st, int nx, int sp, int sh)
 | 
|---|
| 34 |       { store=st; n=nx; spacing=sp; shift=sh; }
 | 
|---|
| 35 | public:
 | 
|---|
| 36 |    Real operator*(const RectMatrixRowCol&) const;         // dot product
 | 
|---|
| 37 |    void AddScaled(const RectMatrixRowCol&, Real);         // add scaled
 | 
|---|
| 38 |    void Divide(const RectMatrixRowCol&, Real);            // scaling
 | 
|---|
| 39 |    void Divide(Real);                                     // scaling
 | 
|---|
| 40 |    void Negate();                                         // change sign
 | 
|---|
| 41 |    void Zero();                                           // zero row col
 | 
|---|
| 42 |    Real& operator[](int i) { return *(store+i*spacing); } // element
 | 
|---|
| 43 |    Real SumSquare() const;                                // sum of squares
 | 
|---|
| 44 |    Real& First() { return *store; }                       // get first element
 | 
|---|
| 45 |    void DownDiag() { store += (shift+spacing); n--; }
 | 
|---|
| 46 |    void UpDiag() { store -= (shift+spacing); n++; }
 | 
|---|
| 47 |    friend void ComplexScale(RectMatrixCol&, RectMatrixCol&, Real, Real);
 | 
|---|
| 48 |    friend void Rotate(RectMatrixCol&, RectMatrixCol&, Real, Real);
 | 
|---|
| 49 |    FREE_CHECK(RectMatrixRowCol)
 | 
|---|
| 50 | };
 | 
|---|
| 51 | 
 | 
|---|
| 52 | /// Access rows of a rectangular matrix.
 | 
|---|
| 53 | /// \internal
 | 
|---|
| 54 | class RectMatrixRow : public RectMatrixRowCol
 | 
|---|
| 55 | {
 | 
|---|
| 56 | public:
 | 
|---|
| 57 |    RectMatrixRow(const Matrix&, int, int, int);
 | 
|---|
| 58 |    RectMatrixRow(const Matrix&, int);
 | 
|---|
| 59 |    void Reset(const Matrix&, int, int, int);
 | 
|---|
| 60 |    void Reset(const Matrix&, int);
 | 
|---|
| 61 |    Real& operator[](int i) { return *(store+i); }
 | 
|---|
| 62 |    void Down() { store += shift; }
 | 
|---|
| 63 |    void Right() { store++; n--; }
 | 
|---|
| 64 |    void Up() { store -= shift; }
 | 
|---|
| 65 |    void Left() { store--; n++; }
 | 
|---|
| 66 |    FREE_CHECK(RectMatrixRow)
 | 
|---|
| 67 | };
 | 
|---|
| 68 | 
 | 
|---|
| 69 | /// Access columns of a rectangular matrix.
 | 
|---|
| 70 | /// \internal
 | 
|---|
| 71 | class RectMatrixCol : public RectMatrixRowCol
 | 
|---|
| 72 | {
 | 
|---|
| 73 | public:
 | 
|---|
| 74 |    RectMatrixCol(const Matrix&, int, int, int);
 | 
|---|
| 75 |    RectMatrixCol(const Matrix&, int);
 | 
|---|
| 76 |    void Reset(const Matrix&, int, int, int);
 | 
|---|
| 77 |    void Reset(const Matrix&, int);
 | 
|---|
| 78 |    void Down() { store += spacing; n--; }
 | 
|---|
| 79 |    void Right() { store++; }
 | 
|---|
| 80 |    void Up() { store -= spacing; n++; }
 | 
|---|
| 81 |    void Left() { store--; }
 | 
|---|
| 82 |    friend void ComplexScale(RectMatrixCol&, RectMatrixCol&, Real, Real);
 | 
|---|
| 83 |    friend void Rotate(RectMatrixCol&, RectMatrixCol&, Real, Real);
 | 
|---|
| 84 |    FREE_CHECK(RectMatrixCol)
 | 
|---|
| 85 | };
 | 
|---|
| 86 | 
 | 
|---|
| 87 | /// Access diagonal of a rectangular matrix.
 | 
|---|
| 88 | /// \internal
 | 
|---|
| 89 | class RectMatrixDiag : public RectMatrixRowCol
 | 
|---|
| 90 | {
 | 
|---|
| 91 | public:
 | 
|---|
| 92 |    RectMatrixDiag(const DiagonalMatrix& D)
 | 
|---|
| 93 |       : RectMatrixRowCol(D.Store(), D.Nrows(), 1, 1) {}
 | 
|---|
| 94 |    Real& operator[](int i) { return *(store+i); }
 | 
|---|
| 95 |    void DownDiag() { store++; n--; }
 | 
|---|
| 96 |    void UpDiag() { store--; n++; }
 | 
|---|
| 97 |    FREE_CHECK(RectMatrixDiag)
 | 
|---|
| 98 | };
 | 
|---|
| 99 | 
 | 
|---|
| 100 | 
 | 
|---|
| 101 | 
 | 
|---|
| 102 | 
 | 
|---|
| 103 | inline RectMatrixRow::RectMatrixRow
 | 
|---|
| 104 |    (const Matrix& M, int row, int skip, int length)
 | 
|---|
| 105 |    : RectMatrixRowCol( M.Store()+row*M.Ncols()+skip, length, 1, M.Ncols() ) {}
 | 
|---|
| 106 | 
 | 
|---|
| 107 | inline RectMatrixRow::RectMatrixRow (const Matrix& M, int row)
 | 
|---|
| 108 |    : RectMatrixRowCol( M.Store()+row*M.Ncols(), M.Ncols(), 1, M.Ncols() ) {}
 | 
|---|
| 109 | 
 | 
|---|
| 110 | inline RectMatrixCol::RectMatrixCol
 | 
|---|
| 111 |    (const Matrix& M, int skip, int col, int length)
 | 
|---|
| 112 |    : RectMatrixRowCol( M.Store()+col+skip*M.Ncols(), length, M.Ncols(), 1 ) {}
 | 
|---|
| 113 | 
 | 
|---|
| 114 | inline RectMatrixCol::RectMatrixCol (const Matrix& M, int col)
 | 
|---|
| 115 |    : RectMatrixRowCol( M.Store()+col, M.Nrows(), M.Ncols(), 1 ) {}
 | 
|---|
| 116 | 
 | 
|---|
| 117 | inline Real square(Real x) { return x*x; }
 | 
|---|
| 118 | inline Real sign(Real x, Real y)
 | 
|---|
| 119 |    { return (y>=0) ? x : -x; }                    // assume x >=0
 | 
|---|
| 120 | 
 | 
|---|
| 121 | 
 | 
|---|
| 122 | // Misc numerical things
 | 
|---|
| 123 | 
 | 
|---|
| 124 | Real pythag(Real f, Real g, Real& c, Real& s);
 | 
|---|
| 125 | 
 | 
|---|
| 126 | inline void GivensRotation(Real cGivens, Real sGivens, Real& x, Real& y)
 | 
|---|
| 127 | {
 | 
|---|
| 128 |    // allow for possibility &x = &y
 | 
|---|
| 129 |    Real tmp0 = cGivens * x + sGivens * y;
 | 
|---|
| 130 |    Real tmp1 = -sGivens * x + cGivens * y;
 | 
|---|
| 131 |    x = tmp0; y = tmp1;
 | 
|---|
| 132 | }
 | 
|---|
| 133 |    
 | 
|---|
| 134 | inline void GivensRotationR(Real cGivens, Real sGivens, Real& x, Real& y)
 | 
|---|
| 135 | {
 | 
|---|
| 136 |    // also change sign of y
 | 
|---|
| 137 |    // allow for possibility &x = &y
 | 
|---|
| 138 |    Real tmp0 = cGivens * x + sGivens * y;
 | 
|---|
| 139 |    Real tmp1 = sGivens * x - cGivens * y;
 | 
|---|
| 140 |    x = tmp0; y = tmp1;
 | 
|---|
| 141 | }   
 | 
|---|
| 142 | 
 | 
|---|
| 143 | 
 | 
|---|
| 144 | 
 | 
|---|
| 145 | 
 | 
|---|
| 146 | 
 | 
|---|
| 147 | #ifdef use_namespace
 | 
|---|
| 148 | }
 | 
|---|
| 149 | #endif
 | 
|---|
| 150 | 
 | 
|---|
| 151 | #endif
 | 
|---|
| 152 | 
 | 
|---|
| 153 | // body file: newmatrm.cpp
 | 
|---|
| 154 | 
 | 
|---|
| 155 | 
 | 
|---|
| 156 | ///@}
 | 
|---|