Changeset 9434 in ntrip


Ignore:
Timestamp:
May 19, 2021, 1:32:38 PM (3 years ago)
Author:
stuerze
Message:

newmat is reverted to its stable version 10 because version 11beta was not working with mac compilers

Location:
trunk/BNC/newmat
Files:
4 deleted
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/BNC/newmat/hholder.cpp

    r8901 r9434  
    335335}
    336336
    337 // Following previous transformation,
    338 // now apply the same orthogonal transformation to (MX & MU)
    339 // Need the X Matrix but not the U.
    340 // Not optimised for accessing consecutive memory
    341 
    342 void updateQRZ(const Matrix& X, Matrix& MX, Matrix& MU)
    343 {
    344    REPORT
    345    Tracer et("updateQRZ(2)");
    346    int s = X.Ncols(); int n = X.Nrows();
    347    if (n != MX.Nrows())
    348       Throw(ProgramException("Incompatible dimensions",X,MX));
    349    if (s != MU.Nrows())
    350       Throw(ProgramException("Incompatible dimensions",X,MU));
    351    int t = MX.Ncols();
    352    if (t != MU.Ncols())
    353       Throw(ProgramException("Incompatible dimensions",MX,MU));
    354    
    355    if (s == 0) return;
    356    
    357    const Real* xi0 = X.data(); Real* mx = MX.data(); Real* muj = MU.data();
    358    for (int i=1; i<=s; ++i)
    359    {
    360                   Real sum = 0.0;
    361                         {
    362          const Real* xi=xi0; int k=n;
    363                            while(k--) { sum += square(*xi); xi+= s;}
    364                         }
    365       Real a0 = sqrt(2.0 - sum); Real* mxj0 = mx;
    366       for (int j=1; j<=t; ++j)
    367       {
    368          Real sum = 0.0;
    369          const Real* xi=xi0; Real* mxj=mxj0; int k=n;
    370          while(--k) { sum += *xi * *mxj; xi += s; mxj += t; }
    371          sum += *xi * *mxj;    // last line of loop
    372          sum += a0 * *muj;
    373          xi=xi0; mxj=mxj0; k=n;
    374          while(--k) { *mxj -= sum * *xi; xi += s; mxj += t; }
    375          *mxj -= sum * *xi;    // last line of loop
    376          *muj -= sum * a0; ++mxj0; ++muj;
    377       }
    378                         ++xi0;
    379    }
    380 }
    381 
    382 
    383 
    384 // same thing as updateQRZ(Matrix& X, UpperTriangularMatrix& U)
    385 // except that X is upper triangular
    386 // contents of X are destroyed - results are in U
    387 // assume we can access efficiently by columns
    388 // e.g. X and U will fit in cache memory
    389 
    390 void updateQRZ(UpperTriangularMatrix& X, UpperTriangularMatrix& U)
    391 {
    392    REPORT
    393    Tracer et("updateQRZ(3)");
    394    int s = X.Ncols();
    395    if (s != U.Ncols())
    396       Throw(ProgramException("Incompatible dimensions",X,U));
    397    if (s == 0) return;
    398    Real* xi0 = X.data(); Real* u = U.data();
    399    for (int i=1; i<=s; ++i)
    400    {
    401                   Real r = *u; Real sum = 0.0;
    402                         {
    403          Real* xi=xi0; int k=i; int l=s;
    404                            while(k--) { sum += square(*xi); xi+= --l;}
    405                         }
    406       sum = sqrt(sum + square(r));
    407       if (sum == 0.0) { REPORT X.column(i) = 0.0; *u = 0.0; }
    408       else
    409       {
    410          Real frs = fabs(r) + sum;
    411          Real a0 = sqrt(frs / sum); Real alpha = a0 / frs;
    412          if (r <= 0) { REPORT *u = sum; alpha = -alpha; }
    413          else { REPORT *u = -sum; }
    414          {
    415             Real* xj0=xi0; int k=i; int l=s;
    416             while(k--) { *xj0 *= alpha; --l; xj0 += l;}
    417          }
    418          Real* xj0=xi0; Real* uj=u;
    419          for (int j=i+1; j<=s; ++j)
    420          {
    421             Real sum = 0.0; ++xj0; ++uj;
    422             Real* xi=xi0; Real* xj=xj0; int k=i; int l=s;
    423             while(k--) { sum += *xi * *xj; --l; xi += l; xj += l; }
    424             sum += a0 * *uj;
    425             xi=xi0; xj=xj0; k=i; l=s;
    426             while(k--) { *xj -= sum * *xi; --l; xi += l; xj += l; }
    427             *uj -= sum * a0;
    428          }
    429       }
    430                         ++xi0; u += s-i+1;
    431    }
    432 }
    433 
    434 // Following previous transformation,
    435 // now apply the same orthogonal transformation to (MX & MU)
    436 // Need the X UpperTriangularMatrix but not the U.
    437 
    438 void updateQRZ(const UpperTriangularMatrix& X, Matrix& MX, Matrix& MU)
    439 {
    440    REPORT
    441    Tracer et("updateQRZ(4)");
    442    int s = X.Ncols();
    443    if (s != MX.Nrows())
    444       Throw(ProgramException("Incompatible dimensions",X,MX));
    445    if (s != MU.Nrows())
    446       Throw(ProgramException("Incompatible dimensions",X,MU));
    447    int t = MX.Ncols();
    448    if (t != MU.Ncols())
    449       Throw(ProgramException("Incompatible dimensions",MX,MU));
    450    if (s == 0) return;
    451    
    452    const Real* xi0 = X.data(); Real* mx = MX.data(); Real* muj = MU.data();
    453    for (int i=1; i<=s; ++i)
    454    {
    455                   Real sum = 0.0;
    456                         {
    457          const Real* xi=xi0; int k=i; int l=s;
    458                            while(k--) { sum += square(*xi); xi+= --l;}
    459                         }
    460       Real a0 = sqrt(2.0 - sum); Real* mxj0 = mx;
    461       for (int j=1; j<=t; ++j)
    462       {
    463          Real sum = 0.0;
    464          const Real* xi=xi0; Real* mxj=mxj0; int k=i; int l=s;
    465          while(--k) { sum += *xi * *mxj; --l; xi += l; mxj += t; }
    466          sum += *xi * *mxj;    // last line of loop
    467          sum += a0 * *muj;
    468          xi=xi0; mxj=mxj0; k=i; l=s;
    469          while(--k) { *mxj -= sum * *xi; --l; xi += l; mxj += t; }
    470          *mxj -= sum * *xi;    // last line of loop
    471          *muj -= sum * a0; ++mxj0; ++muj;
    472       }
    473                         ++xi0;
    474    }
    475 }
    476 
    477 
    478 
    479 
    480 
    481337// Matrix A's first n columns are orthonormal
    482338// so A.Columns(1,n).t() * A.Columns(1,n) is the identity matrix.

  • trunk/BNC/newmat/include.h

    r8901 r9434  
    1 /// \defgroup rbd_common RBD common library
     1/// \defgroup rbd_common RBD common library 
    22///@{
    33
    … …  
    4444//#define HAS_INT64                     // if unsigned _int64 is recognised
    4545                                        // used by newran03
    46 
    47 //#define set_unix_options              // set if running UNIX or LINUX
    48 
     46                                       
    4947// comment out next line if Exception causes a problem
    5048#define TypeDefException
    … …  
    6563#endif
    6664
    67 // for Intel C++ for Windows
    68 #if defined __ICL
    69    #define _STANDARD_                   // use standard library
    70    #define ios_format_flags ios::fmtflags
    71 #endif
    72 
    7365// for Microsoft Visual C++ 7 and above (and Intel simulating these)
    7466#if defined _MSC_VER && _MSC_VER >= 1300
    7567   #define _STANDARD_                   // use standard library
    76 #endif
    77 
    78 // for Borland Builder C++ 2006 and above
    79 #if defined __BCPLUSPLUS__ && __BCPLUSPLUS__ >= 0x0570
    80    #define _STANDARD_                   // use standard library
    81    #define ios_format_flags ios::fmtflags
    8268#endif
    8369
    … …  
    10490      #include <fstream>
    10591   #endif
    106    using namespace std;
     92   ////using namespace std;
     93   #define USE_STD_NAMESPACE
    10794#else
    10895

  • trunk/BNC/newmat/myexcept.cpp

    r8901 r9434  
    2626#endif
    2727
     28#ifdef USE_STD_NAMESPACE
     29using namespace std;
     30#endif
    2831
    2932//#define REG_DEREG                    // for print out uses of new/delete

  • trunk/BNC/newmat/myexcept.h

    r8901 r9434  
    7575   static void AddTrace();               // insert trace in exception record
    7676   static Tracer* last;                  // points to Tracer list
    77    static void clear() {}                // for compatibility
    7877   friend class BaseException;
    7978};
    … …  
    9493   static const char* what() { return what_error; }
    9594                                         // for getting error message
    96    static void clear() {}                // for compatibility
    9795};
    9896

  • trunk/BNC/newmat/newmat.h

    r8901 r9434  
    447447class GeneralMatrix : public BaseMatrix         // declarable matrix types
    448448{
     449   virtual GeneralMatrix* Image() const;        // copy of matrix
    449450protected:
    450451   int tag_val;                                 // shows whether can reuse
    … …  
    464465   void ReverseElements();                      // internal reverse of elements
    465466   void ReverseElements(GeneralMatrix*);        // reverse order of elements
     467   void operator=(Real);                        // set matrix to constant
    466468   Real* GetStore();                            // get store or copy
    467469   GeneralMatrix* BorrowStore(GeneralMatrix*, MatrixType);
    468470                                                // temporarily access store
    469471   void GetMatrix(const GeneralMatrix*);        // used by = and initialise
     472   void Eq(const BaseMatrix&, MatrixType);      // used by =
     473   void Eq(const GeneralMatrix&);               // version with no conversion
     474   void Eq(const BaseMatrix&, MatrixType, bool);// used by <<
     475   void Eq2(const BaseMatrix&, MatrixType);     // cut down version of Eq
    470476   int search(const BaseMatrix*) const;
    471477   virtual GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
    … …  
    478484             // CleanUp when the data array has already been deleted
    479485   void PlusEqual(const GeneralMatrix& gm);
    480    void SP_Equal(const GeneralMatrix& gm);
    481486   void MinusEqual(const GeneralMatrix& gm);
    482487   void PlusEqual(Real f);
    … …  
    485490public:
    486491   GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
    487    void Eq(const BaseMatrix&, MatrixType);      // used by =
    488    void Eq(const GeneralMatrix&);               // version with no conversion
    489    void Eq(const BaseMatrix&, MatrixType, bool);// used by <<
    490    void Eq2(const BaseMatrix&, MatrixType);     // cut down version of Eq
    491492   virtual MatrixType type() const = 0;         // type of a matrix
    492493   MatrixType Type() const { return type(); }
    … …  
    502503   const Real* data() const { return store; }
    503504   const Real* const_data() const { return store; }
    504    void operator=(Real);                        // set matrix to constant
    505505   virtual ~GeneralMatrix();                    // delete store if set
    506506   void tDelete();                              // delete if tag_val permits
    … …  
    526526   void Inject(const GeneralMatrix& GM) { inject(GM); }
    527527   void operator+=(const BaseMatrix&);
    528    void SP_eq(const BaseMatrix&);
    529528   void operator-=(const BaseMatrix&);
    530529   void operator*=(const BaseMatrix&);
    … …  
    580579//   ReturnMatrix Reverse() const;                // reverse order of elements
    581580   void cleanup();                              // to clear store
    582    virtual GeneralMatrix* Image() const;        // copy of matrix
    583581
    584582   friend class Matrix;
    … …  
    627625class Matrix : public GeneralMatrix
    628626{
     627   GeneralMatrix* Image() const;                // copy of matrix
    629628public:
    630629   Matrix() {}
    … …  
    669668   Real minimum2(int& i, int& j) const;
    670669   void operator+=(const Matrix& M) { PlusEqual(M); }
    671    void SP_eq(const Matrix& M) { SP_Equal(M); }
    672670   void operator-=(const Matrix& M) { MinusEqual(M); }
    673    void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
    674    void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
    675    void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    676671   void operator+=(Real f) { GeneralMatrix::Add(f); }
    677672   void operator-=(Real f) { GeneralMatrix::Add(-f); }
    678673   void swap(Matrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
    679    GeneralMatrix* Image() const;                // copy of matrix
    680674   friend Real dotproduct(const Matrix& A, const Matrix& B);
    681675   NEW_DELETE(Matrix)
    … …  
    685679class SquareMatrix : public Matrix
    686680{
     681   GeneralMatrix* Image() const;                // copy of matrix
    687682public:
    688683   SquareMatrix() {}
    … …  
    706701   void ReSize(const GeneralMatrix& A) { resize(A); }
    707702   void operator+=(const Matrix& M) { PlusEqual(M); }
    708    void SP_eq(const Matrix& M) { SP_Equal(M); }
    709703   void operator-=(const Matrix& M) { MinusEqual(M); }
    710    void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
    711    void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
    712    void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    713704   void operator+=(Real f) { GeneralMatrix::Add(f); }
    714705   void operator-=(Real f) { GeneralMatrix::Add(-f); }
    715706   void swap(SquareMatrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
    716    GeneralMatrix* Image() const;                // copy of matrix
    717707   NEW_DELETE(SquareMatrix)
    718708};
    … …  
    721711class nricMatrix : public Matrix
    722712{
     713   GeneralMatrix* Image() const;                // copy of matrix
    723714   Real** row_pointer;                          // points to rows
    724715   void MakeRowPointer();                       // build rowpointer
    … …  
    752743   void MiniCleanUp();
    753744   void operator+=(const Matrix& M) { PlusEqual(M); }
    754    void SP_eq(const Matrix& M) { SP_Equal(M); }
    755745   void operator-=(const Matrix& M) { MinusEqual(M); }
    756    void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
    757    void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
    758    void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    759746   void operator+=(Real f) { GeneralMatrix::Add(f); }
    760747   void operator-=(Real f) { GeneralMatrix::Add(-f); }
    761748   void swap(nricMatrix& gm);
    762    GeneralMatrix* Image() const;                // copy of matrix
    763749   NEW_DELETE(nricMatrix)
    764750};
    … …  
    767753class SymmetricMatrix : public GeneralMatrix
    768754{
     755   GeneralMatrix* Image() const;                // copy of matrix
    769756public:
    770757   SymmetricMatrix() {}
    … …  
    802789   void ReSize(const GeneralMatrix& A) { resize(A); }
    803790   void operator+=(const SymmetricMatrix& M) { PlusEqual(M); }
    804    void SP_eq(const SymmetricMatrix& M) { SP_Equal(M); }
    805791   void operator-=(const SymmetricMatrix& M) { MinusEqual(M); }
    806    void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
    807    void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
    808    void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    809792   void operator+=(Real f) { GeneralMatrix::Add(f); }
    810793   void operator-=(Real f) { GeneralMatrix::Add(-f); }
    811794   void swap(SymmetricMatrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
    812    GeneralMatrix* Image() const;                // copy of matrix
    813795   NEW_DELETE(SymmetricMatrix)
    814796};
    … …  
    817799class UpperTriangularMatrix : public GeneralMatrix
    818800{
     801   GeneralMatrix* Image() const;                // copy of matrix
    819802public:
    820803   UpperTriangularMatrix() {}
    … …  
    854837   MatrixBandWidth bandwidth() const;
    855838   void operator+=(const UpperTriangularMatrix& M) { PlusEqual(M); }
    856    void SP_eq(const UpperTriangularMatrix& M) { SP_Equal(M); }
    857839   void operator-=(const UpperTriangularMatrix& M) { MinusEqual(M); }
    858    void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
    859    void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
    860    void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    861840   void operator+=(Real f) { GeneralMatrix::operator+=(f); }
    862841   void operator-=(Real f) { GeneralMatrix::operator-=(f); }
    863842   void swap(UpperTriangularMatrix& gm)
    864843      { GeneralMatrix::swap((GeneralMatrix&)gm); }
    865    GeneralMatrix* Image() const;                // copy of matrix
    866844   NEW_DELETE(UpperTriangularMatrix)
    867845};
    … …  
    870848class LowerTriangularMatrix : public GeneralMatrix
    871849{
     850   GeneralMatrix* Image() const;                // copy of matrix
    872851public:
    873852   LowerTriangularMatrix() {}
    … …  
    906885   MatrixBandWidth bandwidth() const;
    907886   void operator+=(const LowerTriangularMatrix& M) { PlusEqual(M); }
    908    void SP_eq(const LowerTriangularMatrix& M) { SP_Equal(M); }
    909887   void operator-=(const LowerTriangularMatrix& M) { MinusEqual(M); }
    910    void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
    911    void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
    912    void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    913888   void operator+=(Real f) { GeneralMatrix::operator+=(f); }
    914889   void operator-=(Real f) { GeneralMatrix::operator-=(f); }
    915890   void swap(LowerTriangularMatrix& gm)
    916891      { GeneralMatrix::swap((GeneralMatrix&)gm); }
    917    GeneralMatrix* Image() const;                // copy of matrix
    918892   NEW_DELETE(LowerTriangularMatrix)
    919893};
    … …  
    922896class DiagonalMatrix : public GeneralMatrix
    923897{
     898   GeneralMatrix* Image() const;                // copy of matrix
    924899public:
    925900   DiagonalMatrix() {}
    … …  
    967942//   ReturnMatrix Reverse() const;                // reverse order of elements
    968943   void operator+=(const DiagonalMatrix& M) { PlusEqual(M); }
    969    void SP_eq(const DiagonalMatrix& M) { SP_Equal(M); }
    970944   void operator-=(const DiagonalMatrix& M) { MinusEqual(M); }
    971    void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
    972    void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
    973    void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    974945   void operator+=(Real f) { GeneralMatrix::operator+=(f); }
    975946   void operator-=(Real f) { GeneralMatrix::operator-=(f); }
    976947   void swap(DiagonalMatrix& gm)
    977948      { GeneralMatrix::swap((GeneralMatrix&)gm); }
    978    GeneralMatrix* Image() const;                // copy of matrix
    979949   NEW_DELETE(DiagonalMatrix)
    980950};
    … …  
    983953class RowVector : public Matrix
    984954{
     955   GeneralMatrix* Image() const;                // copy of matrix
    985956public:
    986957   RowVector() { nrows_val = 1; }
    … …  
    1026997   // friend ReturnMatrix GetMatrixRow(Matrix& A, int row);
    1027998   void operator+=(const Matrix& M) { PlusEqual(M); }
    1028    void SP_eq(const Matrix& M) { SP_Equal(M); }
    1029999   void operator-=(const Matrix& M) { MinusEqual(M); }
    1030    void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
    1031    void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
    1032    void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    10331000   void operator+=(Real f) { GeneralMatrix::Add(f); }
    10341001   void operator-=(Real f) { GeneralMatrix::Add(-f); }
    10351002   void swap(RowVector& gm)
    10361003      { GeneralMatrix::swap((GeneralMatrix&)gm); }
    1037    GeneralMatrix* Image() const;                // copy of matrix
    10381004   NEW_DELETE(RowVector)
    10391005};
    … …  
    10421008class ColumnVector : public Matrix
    10431009{
     1010   GeneralMatrix* Image() const;                // copy of matrix
    10441011public:
    10451012   ColumnVector() { ncols_val = 1; }
    … …  
    10791046//   ReturnMatrix Reverse() const;                // reverse order of elements
    10801047   void operator+=(const Matrix& M) { PlusEqual(M); }
    1081    void SP_eq(const Matrix& M) { SP_Equal(M); }
    10821048   void operator-=(const Matrix& M) { MinusEqual(M); }
    1083    void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
    1084    void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
    1085    void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    10861049   void operator+=(Real f) { GeneralMatrix::Add(f); }
    10871050   void operator-=(Real f) { GeneralMatrix::Add(-f); }
    10881051   void swap(ColumnVector& gm)
    10891052      { GeneralMatrix::swap((GeneralMatrix&)gm); }
    1090    GeneralMatrix* Image() const;                // copy of matrix
    10911053   NEW_DELETE(ColumnVector)
    10921054};
    … …  
    11021064   void ludcmp();
    11031065   void get_aux(CroutMatrix&);                  // for copying indx[] etc
     1066   GeneralMatrix* Image() const;                // copy of matrix
    11041067public:
    11051068   CroutMatrix(const BaseMatrix&);
    … …  
    11251088   bool even_exchanges() const { return d; }
    11261089   void swap(CroutMatrix& gm);
    1127    GeneralMatrix* Image() const;                // copy of matrix
    11281090   NEW_DELETE(CroutMatrix)
    11291091};
    … …  
    11341096class BandMatrix : public GeneralMatrix
    11351097{
     1098   GeneralMatrix* Image() const;                // copy of matrix
    11361099protected:
    11371100   void CornerClear() const;                    // set unused elements to zero
    … …  
    11981161   void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
    11991162   void swap(BandMatrix& gm);
    1200    GeneralMatrix* Image() const;                // copy of matrix
    12011163   NEW_DELETE(BandMatrix)
    12021164};
    … …  
    12051167class UpperBandMatrix : public BandMatrix
    12061168{
     1169   GeneralMatrix* Image() const;                // copy of matrix
    12071170public:
    12081171   UpperBandMatrix() {}
    … …  
    12371200   void swap(UpperBandMatrix& gm)
    12381201      { BandMatrix::swap((BandMatrix&)gm); }
    1239    GeneralMatrix* Image() const;                // copy of matrix
    12401202   NEW_DELETE(UpperBandMatrix)
    12411203};
    … …  
    12441206class LowerBandMatrix : public BandMatrix
    12451207{
     1208   GeneralMatrix* Image() const;                // copy of matrix
    12461209public:
    12471210   LowerBandMatrix() {}
    … …  
    12761239   void swap(LowerBandMatrix& gm)
    12771240      { BandMatrix::swap((BandMatrix&)gm); }
    1278    GeneralMatrix* Image() const;                // copy of matrix
    12791241   NEW_DELETE(LowerBandMatrix)
    12801242};
    … …  
    12831245class SymmetricBandMatrix : public GeneralMatrix
    12841246{
     1247   GeneralMatrix* Image() const;                // copy of matrix
    12851248   void CornerClear() const;                    // set unused elements to zero
    12861249   short SimpleAddOK(const GeneralMatrix* gm);
    … …  
    13371300   void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
    13381301   void swap(SymmetricBandMatrix& gm);
    1339    GeneralMatrix* Image() const;                // copy of matrix
    13401302   NEW_DELETE(SymmetricBandMatrix)
    13411303};
    … …  
    13521314   void ludcmp();
    13531315   void get_aux(BandLUMatrix&);                 // for copying indx[] etc
     1316   GeneralMatrix* Image() const;                // copy of matrix
    13541317public:
    13551318   BandLUMatrix(const BaseMatrix&);
    … …  
    13791342   MatrixBandWidth bandwidth() const;
    13801343   void swap(BandLUMatrix& gm);
    1381    GeneralMatrix* Image() const;                // copy of matrix
    13821344   NEW_DELETE(BandLUMatrix)
    13831345};
    … …  
    13881350class IdentityMatrix : public GeneralMatrix
    13891351{
     1352   GeneralMatrix* Image() const;          // copy of matrix
    13901353public:
    13911354   IdentityMatrix() {}
    … …  
    14231386   void swap(IdentityMatrix& gm)
    14241387      { GeneralMatrix::swap((GeneralMatrix&)gm); }
    1425    GeneralMatrix* Image() const;          // copy of matrix
    14261388   NEW_DELETE(IdentityMatrix)
    14271389};
    … …  
    14471409   void operator=(const BaseMatrix&);
    14481410   void operator+=(const BaseMatrix&);
    1449    void SP_eq(const BaseMatrix&);
    14501411   void operator-=(const BaseMatrix&);
    14511412   void operator*=(const BaseMatrix&);
    … …  
    18201781   void operator=(const BaseMatrix&);
    18211782   void operator+=(const BaseMatrix&);
    1822    void SP_eq(const BaseMatrix&);
    18231783   void operator-=(const BaseMatrix&);
    18241784   void operator=(const GetSubMatrix& m) { operator=((const BaseMatrix&)m); }

  • trunk/BNC/newmat/newmat.pro

    r8901 r9434  
    1010MOC_DIR     = .moc
    1111
    12 HEADERS += boolean.h controlw.h include.h myexcept.h  \
    13            newmatap.h newmat.h newmatio.h newmatnl.h  \
    14            newmatrc.h newmatrm.h precisio.h solution.h
     12HEADERS += controlw.h include.h myexcept.h  \
     13           newmatap.h newmat.h newmatio.h   \
     14           newmatrc.h newmatrm.h precisio.h
    1515
    16 SOURCES += bandmat.cpp cholesky.cpp evalue.cpp     \
    17            fft.cpp hholder.cpp jacobi.cpp          \
    18            myexcept.cpp newfft.cpp newmat1.cpp     \
    19            newmat2.cpp newmat3.cpp newmat4.cpp     \
    20            newmat5.cpp newmat6.cpp newmat7.cpp     \
    21            newmat8.cpp newmat9.cpp newmatex.cpp    \
    22            newmatrm.cpp newmatnl.cpp nm_misc.cpp   \
    23            sort.cpp submat.cpp svd.cpp solution.cpp
     16SOURCES += bandmat.cpp cholesky.cpp evalue.cpp  \
     17           fft.cpp hholder.cpp jacobi.cpp       \
     18           myexcept.cpp newfft.cpp newmat1.cpp  \
     19           newmat2.cpp newmat3.cpp newmat4.cpp  \
     20           newmat5.cpp newmat6.cpp newmat7.cpp  \
     21           newmat8.cpp newmat9.cpp newmatex.cpp \
     22           newmatrm.cpp nm_misc.cpp sort.cpp    \
     23           submat.cpp svd.cpp
    2424

  • trunk/BNC/newmat/newmat6.cpp

    r8901 r9434  
    495495   Tracer tr("GeneralMatrix::operator+=");
    496496   // MatrixConversionCheck mcc;
    497    Protect();  // so it cannot get deleted during Evaluate
     497   Protect();                                   // so it cannot get deleted
     498                                                // during Evaluate
    498499   GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
    499500   AddedMatrix am(this,gm);
    500501   if (gm==this) Release(2); else Release();
    501502   Eq2(am,type());
    502 }
    503 
    504 // GeneralMatrix operators
    505 
    506 void GeneralMatrix::SP_eq(const BaseMatrix& X)
    507 {
    508    REPORT
    509    Tracer tr("GeneralMatrix::SP_eq");
    510    // MatrixConversionCheck mcc;
    511    Protect();  // so it cannot get deleted during Evaluate
    512    GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
    513    SPMatrix spm(this,gm);
    514    if (gm==this) Release(2); else Release();
    515    Eq2(spm,type());
    516503}
    517504
    … …  
    599586   if (gmx==gm) gm->Release(2); else gm->Release();
    600587   GeneralMatrix* gmy = am.Evaluate();
    601    if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
    602    else { REPORT }
    603    gm->Protect();
    604 }
    605 
    606 void GenericMatrix::SP_eq(const BaseMatrix& X)
    607 {
    608    REPORT
    609    Tracer tr("GenericMatrix::SP_eq");
    610    if (!gm) Throw(ProgramException("GenericMatrix is null"));
    611    gm->Protect();            // so it cannot get deleted during Evaluate
    612    GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
    613    SPMatrix spm(gm,gmx);
    614    if (gmx==gm) gm->Release(2); else gm->Release();
    615    GeneralMatrix* gmy = spm.Evaluate();
    616588   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
    617589   else { REPORT }

  • trunk/BNC/newmat/newmat7.cpp

    r8901 r9434  
    214214}
    215215
    216 static void SP(GeneralMatrix* gm, const GeneralMatrix* gm2)
    217 {
    218    REPORT
    219    const Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
     216static void SP(GeneralMatrix* gm, GeneralMatrix* gm2)
     217{
     218   REPORT
     219   Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
    220220   while (i--)
    221221   { *s++ *= *s2++; *s++ *= *s2++; *s++ *= *s2++; *s++ *= *s2++; }
    222222   i=gm->Storage() & 3; while (i--) *s++ *= *s2++;
    223 }
    224 
    225 void GeneralMatrix::SP_Equal(const GeneralMatrix& gm)
    226 {
    227    REPORT
    228    if (nrows_val != gm.nrows_val || ncols_val != gm.ncols_val)
    229       Throw(IncompatibleDimensionsException(*this, gm));
    230    SP(this, &gm);
    231223}
    232224

  • trunk/BNC/newmat/newmat9.cpp

    r8901 r9434  
    4444   MatrixRow mr((GeneralMatrix*)&X, LoadOnEntry);
    4545   int w = s.width();  int nr = X.Nrows();  ios_format_flags f = s.flags();
    46    if (f & ios::scientific) s.setf(ios::scientific, ios::floatfield);
    47    else s.setf(ios::fixed, ios::floatfield);
     46   s.setf(ios::fixed, ios::floatfield);
    4847   for (int i=1; i<=nr; i++)
    4948   {

  • trunk/BNC/newmat/newmatap.h

    r8901 r9434  
    2626void QRZ(Matrix&, UpperTriangularMatrix&);
    2727
    28 void QRZ(const Matrix&, Matrix&, Matrix&);
    29 
    30 inline void QRZT(Matrix& X, Matrix& Y, LowerTriangularMatrix& L, Matrix& M)
    31    { QRZT(X, L); QRZT(X, Y, M); }
    32 
    33 inline void QRZ(Matrix& X, Matrix& Y, UpperTriangularMatrix& U, Matrix& M)
    34    { QRZ(X, U); QRZ(X, Y, M); }
     28void QRZ(const Matrix&, Matrix&, Matrix&);
    3529
    3630inline void HHDecompose(Matrix& X, LowerTriangularMatrix& L)
    … …  
    4438void updateQRZ(Matrix& X, UpperTriangularMatrix& U);
    4539
    46 void updateQRZ(const Matrix& X, Matrix& MX, Matrix& MU);
    47 
    48 void updateQRZ(UpperTriangularMatrix& X, UpperTriangularMatrix& U);
    49 
    50 void updateQRZ(const UpperTriangularMatrix& X, Matrix& MX, Matrix& MU);
    51 
    5240inline void UpdateQRZT(Matrix& X, LowerTriangularMatrix& L)
    5341   { updateQRZT(X, L); }
    … …  
    5543inline void UpdateQRZ(Matrix& X, UpperTriangularMatrix& U)
    5644   { updateQRZ(X, U); }
    57 
    58 inline void UpdateQRZ(UpperTriangularMatrix& X, UpperTriangularMatrix& U)
    59    { updateQRZ(X, U); }
    60 
    61 inline void UpdateQRZ(const UpperTriangularMatrix& X, Matrix& MX, Matrix& MU)
    62    { updateQRZ(X, MX, MU); }
    63    
    64 inline void UpdateQRZ(const Matrix& X, Matrix& MX, Matrix& MU)
    65    { updateQRZ(X, MX, MU); }
    66 
    6745
    6846// Matrix A's first n columns are orthonormal
    … …  
    7856
    7957
    80 // produces the Cholesky decomposition of A + x.t() * x
    81 // where A = chol.t() * chol and x is a RowVector
     58// produces the Cholesky decomposition of A + x.t() * x where A = chol.t() * chol
     59// and x is a RowVector
    8260void update_Cholesky(UpperTriangularMatrix& chol, RowVector x);
    8361inline void UpdateCholesky(UpperTriangularMatrix& chol, const RowVector& x)
    8462   { update_Cholesky(chol, x); }
    8563
    86 // produces the Cholesky decomposition of A - x.t() * x
    87 // where A = chol.t() * chol and x is a RowVector
     64// produces the Cholesky decomposition of A - x.t() * x where A = chol.t() * chol
     65// and x is a RowVector
    8866void downdate_Cholesky(UpperTriangularMatrix &chol, RowVector x);
    8967inline void DowndateCholesky(UpperTriangularMatrix &chol, const RowVector& x)

  • trunk/BNC/newmat/newmatio.h

    r8901 r9434  
    1616#include "newmat.h"
    1717
    18 #ifdef use_namespace
    19 namespace NEWMAT {
     18#ifdef USE_STD_NAMESPACE
     19using namespace std;
    2020#endif
    21 
    22 
    2321
    2422// **************************** input/output *****************************/

  • trunk/BNC/newmat/precisio.h

    r8901 r9434  
    2525#endif
    2626
    27 using namespace std;
    28        
    2927/// Floating point precision.
    3028class FloatingPointPrecision
    … …  
    3230public:
    3331   static int Dig()              // number of decimal digits or precision
    34       { return numeric_limits<Real>::digits10 ; }
     32      { return std::numeric_limits<Real>::digits10 ; }
    3533
    3634   static Real Epsilon()         // smallest number such that 1+Eps!=Eps
    37       { return numeric_limits<Real>::epsilon(); }
     35      { return std::numeric_limits<Real>::epsilon(); }
    3836
    3937   static int Mantissa()         // bits in mantisa
    40       { return numeric_limits<Real>::digits; }
     38      { return std::numeric_limits<Real>::digits; }
    4139
    4240   static Real Maximum()         // maximum value
    43       { return numeric_limits<Real>::max(); }
     41      { return std::numeric_limits<Real>::max(); }
    4442
    4543   static int MaximumDecimalExponent()  // maximum decimal exponent
    46       { return numeric_limits<Real>::max_exponent10; }
     44      { return std::numeric_limits<Real>::max_exponent10; }
    4745
    4846   static int MaximumExponent()  // maximum binary exponent
    49       { return numeric_limits<Real>::max_exponent; }
     47      { return std::numeric_limits<Real>::max_exponent; }
    5048
    5149   static Real LnMaximum()       // natural log of maximum
    … …  
    5351
    5452   static Real Minimum()         // minimum positive value
    55       { return numeric_limits<Real>::min(); }
     53      { return std::numeric_limits<Real>::min(); }
    5654
    5755   static int MinimumDecimalExponent() // minimum decimal exponent
    58       { return numeric_limits<Real>::min_exponent10; }
     56      { return std::numeric_limits<Real>::min_exponent10; }
    5957
    6058   static int MinimumExponent()  // minimum binary exponent
    61       { return numeric_limits<Real>::min_exponent; }
     59      { return std::numeric_limits<Real>::min_exponent; }
    6260
    6361   static Real LnMinimum()       // natural log of minimum
    … …  
    6563
    6664   static int Radix()            // exponent radix
    67       { return numeric_limits<Real>::radix; }
     65      { return std::numeric_limits<Real>::radix; }
    6866
    6967   static int Rounds()           // addition rounding (1 = does round)
    7068   {
    71           return numeric_limits<Real>::round_style ==
    72                  round_to_nearest ? 1 : 0;
     69          return std::numeric_limits<Real>::round_style ==
     70                 std::round_to_nearest ? 1 : 0;
    7371   }
    7472

  • trunk/BNC/newmat/submat.cpp

    r8901 r9434  
    257257      if (row_number != gmx->Nrows() || col_number != gmx->Ncols())
    258258         Throw(IncompatibleDimensionsException());
    259       if (gm->type().is_symmetric() &&
    260          ( ! gmx->type().is_symmetric() || row_skip != col_skip) )
    261          Throw(ProgramException("Illegal operation on symmetric"));
    262259      MatrixRow mrx(gmx, LoadOnEntry);
    263260      MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
    … …  
    280277}
    281278
    282 void GetSubMatrix::SP_eq(const BaseMatrix& bmx)
    283 {
    284    REPORT
    285    Tracer tr("SubMatrix(SP_eq)"); GeneralMatrix* gmx = 0;
     279void GetSubMatrix::operator-=(const BaseMatrix& bmx)
     280{
     281   REPORT
     282   Tracer tr("SubMatrix(-=)"); GeneralMatrix* gmx = 0;
    286283   // MatrixConversionCheck mcc;         // Check for loss of info
    287284   Try
    … …  
    290287      if (row_number != gmx->Nrows() || col_number != gmx->Ncols())
    291288         Throw(IncompatibleDimensionsException());
    292       if (gm->type().is_symmetric() &&
    293          ( ! gmx->type().is_symmetric() || row_skip != col_skip) )
    294          Throw(ProgramException("Illegal operation on symmetric"));
    295       MatrixRow mrx(gmx, LoadOnEntry);
    296       MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
    297                                      // do need LoadOnEntry
    298       MatrixRowCol sub; int i = row_number;
    299       while (i--)
    300       {
    301          mr.SubRowCol(sub, col_skip, col_number);   // put values in sub
    302          sub.Multiply(mrx); mr.Next(); mrx.Next();
    303       }
    304       gmx->tDelete();
    305    }
    306 
    307    CatchAll
    308    {
    309       if (gmx) gmx->tDelete();
    310       ReThrow;
    311    }
    312 }
    313 
    314 void GetSubMatrix::operator-=(const BaseMatrix& bmx)
    315 {
    316    REPORT
    317    Tracer tr("SubMatrix(-=)"); GeneralMatrix* gmx = 0;
    318    // MatrixConversionCheck mcc;         // Check for loss of info
    319    Try
    320    {
    321       SetUpLHS(); gmx = ((BaseMatrix&)bmx).Evaluate();
    322       if (row_number != gmx->Nrows() || col_number != gmx->Ncols())
    323          Throw(IncompatibleDimensionsException());
    324       if (gm->type().is_symmetric() &&
    325          ( ! gmx->type().is_symmetric() || row_skip != col_skip) )
    326          Throw(ProgramException("Illegal operation on symmetric"));
    327289      MatrixRow mrx(gmx, LoadOnEntry);
    328290      MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
Note: See TracChangeset for help on using the changeset viewer.