Changeset 8901 in ntrip


Ignore:
Timestamp:
Mar 18, 2020, 11:06:13 AM (5 years ago)
Author:
stuerze
Message:

upgrade to newmat11 library

Location:
trunk/BNC
Files:
4 added
16 edited

Legend:

Unmodified
Added
Removed

  • trunk/BNC/newmat/hholder.cpp

    r2013 r8901  
    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
     342void 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
     390void 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
     438void 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
    337481// Matrix A's first n columns are orthonormal
    338482// so A.Columns(1,n).t() * A.Columns(1,n) is the identity matrix.

  • trunk/BNC/newmat/include.h

    r2013 r8901  
    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                                        
     46
     47//#define set_unix_options              // set if running UNIX or LINUX
     48
    4749// comment out next line if Exception causes a problem
    4850#define TypeDefException
    … …  
    6365#endif
    6466
     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
    6573// for Microsoft Visual C++ 7 and above (and Intel simulating these)
    6674#if defined _MSC_VER && _MSC_VER >= 1300
    6775   #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
    6882#endif
    6983
    … …  
    90104      #include <fstream>
    91105   #endif
    92    ////using namespace std;
    93    #define USE_STD_NAMESPACE
     106   using namespace std;
    94107#else
    95108

  • trunk/BNC/newmat/myexcept.cpp

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

  • trunk/BNC/newmat/myexcept.h

    r2013 r8901  
    7575   static void AddTrace();               // insert trace in exception record
    7676   static Tracer* last;                  // points to Tracer list
     77   static void clear() {}                // for compatibility
    7778   friend class BaseException;
    7879};
    … …  
    9394   static const char* what() { return what_error; }
    9495                                         // for getting error message
     96   static void clear() {}                // for compatibility
    9597};
    9698

  • trunk/BNC/newmat/newmat.h

    r2013 r8901  
    447447class GeneralMatrix : public BaseMatrix         // declarable matrix types
    448448{
    449    virtual GeneralMatrix* Image() const;        // copy of matrix
    450449protected:
    451450   int tag_val;                                 // shows whether can reuse
    … …  
    465464   void ReverseElements();                      // internal reverse of elements
    466465   void ReverseElements(GeneralMatrix*);        // reverse order of elements
    467    void operator=(Real);                        // set matrix to constant
    468466   Real* GetStore();                            // get store or copy
    469467   GeneralMatrix* BorrowStore(GeneralMatrix*, MatrixType);
    470468                                                // temporarily access store
    471469   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
    476470   int search(const BaseMatrix*) const;
    477471   virtual GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
    … …  
    484478             // CleanUp when the data array has already been deleted
    485479   void PlusEqual(const GeneralMatrix& gm);
     480   void SP_Equal(const GeneralMatrix& gm);
    486481   void MinusEqual(const GeneralMatrix& gm);
    487482   void PlusEqual(Real f);
    … …  
    490485public:
    491486   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
    492491   virtual MatrixType type() const = 0;         // type of a matrix
    493492   MatrixType Type() const { return type(); }
    … …  
    503502   const Real* data() const { return store; }
    504503   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&);
    528529   void operator-=(const BaseMatrix&);
    529530   void operator*=(const BaseMatrix&);
    … …  
    579580//   ReturnMatrix Reverse() const;                // reverse order of elements
    580581   void cleanup();                              // to clear store
     582   virtual GeneralMatrix* Image() const;        // copy of matrix
    581583
    582584   friend class Matrix;
    … …  
    625627class Matrix : public GeneralMatrix
    626628{
    627    GeneralMatrix* Image() const;                // copy of matrix
    628629public:
    629630   Matrix() {}
    … …  
    668669   Real minimum2(int& i, int& j) const;
    669670   void operator+=(const Matrix& M) { PlusEqual(M); }
     671   void SP_eq(const Matrix& M) { SP_Equal(M); }
    670672   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); }
    671676   void operator+=(Real f) { GeneralMatrix::Add(f); }
    672677   void operator-=(Real f) { GeneralMatrix::Add(-f); }
    673678   void swap(Matrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
     679   GeneralMatrix* Image() const;                // copy of matrix
    674680   friend Real dotproduct(const Matrix& A, const Matrix& B);
    675681   NEW_DELETE(Matrix)
    … …  
    679685class SquareMatrix : public Matrix
    680686{
    681    GeneralMatrix* Image() const;                // copy of matrix
    682687public:
    683688   SquareMatrix() {}
    … …  
    701706   void ReSize(const GeneralMatrix& A) { resize(A); }
    702707   void operator+=(const Matrix& M) { PlusEqual(M); }
     708   void SP_eq(const Matrix& M) { SP_Equal(M); }
    703709   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); }
    704713   void operator+=(Real f) { GeneralMatrix::Add(f); }
    705714   void operator-=(Real f) { GeneralMatrix::Add(-f); }
    706715   void swap(SquareMatrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
     716   GeneralMatrix* Image() const;                // copy of matrix
    707717   NEW_DELETE(SquareMatrix)
    708718};
    … …  
    711721class nricMatrix : public Matrix
    712722{
    713    GeneralMatrix* Image() const;                // copy of matrix
    714723   Real** row_pointer;                          // points to rows
    715724   void MakeRowPointer();                       // build rowpointer
    … …  
    743752   void MiniCleanUp();
    744753   void operator+=(const Matrix& M) { PlusEqual(M); }
     754   void SP_eq(const Matrix& M) { SP_Equal(M); }
    745755   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); }
    746759   void operator+=(Real f) { GeneralMatrix::Add(f); }
    747760   void operator-=(Real f) { GeneralMatrix::Add(-f); }
    748761   void swap(nricMatrix& gm);
     762   GeneralMatrix* Image() const;                // copy of matrix
    749763   NEW_DELETE(nricMatrix)
    750764};
    … …  
    753767class SymmetricMatrix : public GeneralMatrix
    754768{
    755    GeneralMatrix* Image() const;                // copy of matrix
    756769public:
    757770   SymmetricMatrix() {}
    … …  
    789802   void ReSize(const GeneralMatrix& A) { resize(A); }
    790803   void operator+=(const SymmetricMatrix& M) { PlusEqual(M); }
     804   void SP_eq(const SymmetricMatrix& M) { SP_Equal(M); }
    791805   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); }
    792809   void operator+=(Real f) { GeneralMatrix::Add(f); }
    793810   void operator-=(Real f) { GeneralMatrix::Add(-f); }
    794811   void swap(SymmetricMatrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
     812   GeneralMatrix* Image() const;                // copy of matrix
    795813   NEW_DELETE(SymmetricMatrix)
    796814};
    … …  
    799817class UpperTriangularMatrix : public GeneralMatrix
    800818{
    801    GeneralMatrix* Image() const;                // copy of matrix
    802819public:
    803820   UpperTriangularMatrix() {}
    … …  
    837854   MatrixBandWidth bandwidth() const;
    838855   void operator+=(const UpperTriangularMatrix& M) { PlusEqual(M); }
     856   void SP_eq(const UpperTriangularMatrix& M) { SP_Equal(M); }
    839857   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); }
    840861   void operator+=(Real f) { GeneralMatrix::operator+=(f); }
    841862   void operator-=(Real f) { GeneralMatrix::operator-=(f); }
    842863   void swap(UpperTriangularMatrix& gm)
    843864      { GeneralMatrix::swap((GeneralMatrix&)gm); }
     865   GeneralMatrix* Image() const;                // copy of matrix
    844866   NEW_DELETE(UpperTriangularMatrix)
    845867};
    … …  
    848870class LowerTriangularMatrix : public GeneralMatrix
    849871{
    850    GeneralMatrix* Image() const;                // copy of matrix
    851872public:
    852873   LowerTriangularMatrix() {}
    … …  
    885906   MatrixBandWidth bandwidth() const;
    886907   void operator+=(const LowerTriangularMatrix& M) { PlusEqual(M); }
     908   void SP_eq(const LowerTriangularMatrix& M) { SP_Equal(M); }
    887909   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); }
    888913   void operator+=(Real f) { GeneralMatrix::operator+=(f); }
    889914   void operator-=(Real f) { GeneralMatrix::operator-=(f); }
    890915   void swap(LowerTriangularMatrix& gm)
    891916      { GeneralMatrix::swap((GeneralMatrix&)gm); }
     917   GeneralMatrix* Image() const;                // copy of matrix
    892918   NEW_DELETE(LowerTriangularMatrix)
    893919};
    … …  
    896922class DiagonalMatrix : public GeneralMatrix
    897923{
    898    GeneralMatrix* Image() const;                // copy of matrix
    899924public:
    900925   DiagonalMatrix() {}
    … …  
    942967//   ReturnMatrix Reverse() const;                // reverse order of elements
    943968   void operator+=(const DiagonalMatrix& M) { PlusEqual(M); }
     969   void SP_eq(const DiagonalMatrix& M) { SP_Equal(M); }
    944970   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); }
    945974   void operator+=(Real f) { GeneralMatrix::operator+=(f); }
    946975   void operator-=(Real f) { GeneralMatrix::operator-=(f); }
    947976   void swap(DiagonalMatrix& gm)
    948977      { GeneralMatrix::swap((GeneralMatrix&)gm); }
     978   GeneralMatrix* Image() const;                // copy of matrix
    949979   NEW_DELETE(DiagonalMatrix)
    950980};
    … …  
    953983class RowVector : public Matrix
    954984{
    955    GeneralMatrix* Image() const;                // copy of matrix
    956985public:
    957986   RowVector() { nrows_val = 1; }
    … …  
    9971026   // friend ReturnMatrix GetMatrixRow(Matrix& A, int row);
    9981027   void operator+=(const Matrix& M) { PlusEqual(M); }
     1028   void SP_eq(const Matrix& M) { SP_Equal(M); }
    9991029   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); }
    10001033   void operator+=(Real f) { GeneralMatrix::Add(f); }
    10011034   void operator-=(Real f) { GeneralMatrix::Add(-f); }
    10021035   void swap(RowVector& gm)
    10031036      { GeneralMatrix::swap((GeneralMatrix&)gm); }
     1037   GeneralMatrix* Image() const;                // copy of matrix
    10041038   NEW_DELETE(RowVector)
    10051039};
    … …  
    10081042class ColumnVector : public Matrix
    10091043{
    1010    GeneralMatrix* Image() const;                // copy of matrix
    10111044public:
    10121045   ColumnVector() { ncols_val = 1; }
    … …  
    10461079//   ReturnMatrix Reverse() const;                // reverse order of elements
    10471080   void operator+=(const Matrix& M) { PlusEqual(M); }
     1081   void SP_eq(const Matrix& M) { SP_Equal(M); }
    10481082   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); }
    10491086   void operator+=(Real f) { GeneralMatrix::Add(f); }
    10501087   void operator-=(Real f) { GeneralMatrix::Add(-f); }
    10511088   void swap(ColumnVector& gm)
    10521089      { GeneralMatrix::swap((GeneralMatrix&)gm); }
     1090   GeneralMatrix* Image() const;                // copy of matrix
    10531091   NEW_DELETE(ColumnVector)
    10541092};
    … …  
    10641102   void ludcmp();
    10651103   void get_aux(CroutMatrix&);                  // for copying indx[] etc
    1066    GeneralMatrix* Image() const;                // copy of matrix
    10671104public:
    10681105   CroutMatrix(const BaseMatrix&);
    … …  
    10881125   bool even_exchanges() const { return d; }
    10891126   void swap(CroutMatrix& gm);
     1127   GeneralMatrix* Image() const;                // copy of matrix
    10901128   NEW_DELETE(CroutMatrix)
    10911129};
    … …  
    10961134class BandMatrix : public GeneralMatrix
    10971135{
    1098    GeneralMatrix* Image() const;                // copy of matrix
    10991136protected:
    11001137   void CornerClear() const;                    // set unused elements to zero
    … …  
    11611198   void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
    11621199   void swap(BandMatrix& gm);
     1200   GeneralMatrix* Image() const;                // copy of matrix
    11631201   NEW_DELETE(BandMatrix)
    11641202};
    … …  
    11671205class UpperBandMatrix : public BandMatrix
    11681206{
    1169    GeneralMatrix* Image() const;                // copy of matrix
    11701207public:
    11711208   UpperBandMatrix() {}
    … …  
    12001237   void swap(UpperBandMatrix& gm)
    12011238      { BandMatrix::swap((BandMatrix&)gm); }
     1239   GeneralMatrix* Image() const;                // copy of matrix
    12021240   NEW_DELETE(UpperBandMatrix)
    12031241};
    … …  
    12061244class LowerBandMatrix : public BandMatrix
    12071245{
    1208    GeneralMatrix* Image() const;                // copy of matrix
    12091246public:
    12101247   LowerBandMatrix() {}
    … …  
    12391276   void swap(LowerBandMatrix& gm)
    12401277      { BandMatrix::swap((BandMatrix&)gm); }
     1278   GeneralMatrix* Image() const;                // copy of matrix
    12411279   NEW_DELETE(LowerBandMatrix)
    12421280};
    … …  
    12451283class SymmetricBandMatrix : public GeneralMatrix
    12461284{
    1247    GeneralMatrix* Image() const;                // copy of matrix
    12481285   void CornerClear() const;                    // set unused elements to zero
    12491286   short SimpleAddOK(const GeneralMatrix* gm);
    … …  
    13001337   void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
    13011338   void swap(SymmetricBandMatrix& gm);
     1339   GeneralMatrix* Image() const;                // copy of matrix
    13021340   NEW_DELETE(SymmetricBandMatrix)
    13031341};
    … …  
    13141352   void ludcmp();
    13151353   void get_aux(BandLUMatrix&);                 // for copying indx[] etc
    1316    GeneralMatrix* Image() const;                // copy of matrix
    13171354public:
    13181355   BandLUMatrix(const BaseMatrix&);
    … …  
    13421379   MatrixBandWidth bandwidth() const;
    13431380   void swap(BandLUMatrix& gm);
     1381   GeneralMatrix* Image() const;                // copy of matrix
    13441382   NEW_DELETE(BandLUMatrix)
    13451383};
    … …  
    13501388class IdentityMatrix : public GeneralMatrix
    13511389{
    1352    GeneralMatrix* Image() const;          // copy of matrix
    13531390public:
    13541391   IdentityMatrix() {}
    … …  
    13861423   void swap(IdentityMatrix& gm)
    13871424      { GeneralMatrix::swap((GeneralMatrix&)gm); }
     1425   GeneralMatrix* Image() const;          // copy of matrix
    13881426   NEW_DELETE(IdentityMatrix)
    13891427};
    … …  
    14091447   void operator=(const BaseMatrix&);
    14101448   void operator+=(const BaseMatrix&);
     1449   void SP_eq(const BaseMatrix&);
    14111450   void operator-=(const BaseMatrix&);
    14121451   void operator*=(const BaseMatrix&);
    … …  
    17811820   void operator=(const BaseMatrix&);
    17821821   void operator+=(const BaseMatrix&);
     1822   void SP_eq(const BaseMatrix&);
    17831823   void operator-=(const BaseMatrix&);
    17841824   void operator=(const GetSubMatrix& m) { operator=((const BaseMatrix&)m); }

  • trunk/BNC/newmat/newmat.pro

    r5182 r8901  
    1010MOC_DIR     = .moc
    1111
    12 HEADERS += controlw.h include.h myexcept.h  \
    13            newmatap.h newmat.h newmatio.h   \
    14            newmatrc.h newmatrm.h precisio.h
     12HEADERS += 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
    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 nm_misc.cpp sort.cpp    \
    23            submat.cpp svd.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 newmatnl.cpp nm_misc.cpp   \
     23           sort.cpp submat.cpp svd.cpp solution.cpp
    2424

  • trunk/BNC/newmat/newmat6.cpp

    r2013 r8901  
    495495   Tracer tr("GeneralMatrix::operator+=");
    496496   // MatrixConversionCheck mcc;
    497    Protect();                                   // so it cannot get deleted
    498                                                 // during Evaluate
     497   Protect();  // so it cannot get deleted during Evaluate
    499498   GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
    500499   AddedMatrix am(this,gm);
    501500   if (gm==this) Release(2); else Release();
    502501   Eq2(am,type());
     502}
     503
     504// GeneralMatrix operators
     505
     506void 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());
    503516}
    504517
    … …  
    586599   if (gmx==gm) gm->Release(2); else gm->Release();
    587600   GeneralMatrix* gmy = am.Evaluate();
     601   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
     602   else { REPORT }
     603   gm->Protect();
     604}
     605
     606void 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();
    588616   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
    589617   else { REPORT }

  • trunk/BNC/newmat/newmat7.cpp

    r2013 r8901  
    214214}
    215215
    216 static void SP(GeneralMatrix* gm, GeneralMatrix* gm2)
    217 {
    218    REPORT
    219    Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
     216static 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;
    220220   while (i--)
    221221   { *s++ *= *s2++; *s++ *= *s2++; *s++ *= *s2++; *s++ *= *s2++; }
    222222   i=gm->Storage() & 3; while (i--) *s++ *= *s2++;
     223}
     224
     225void 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);
    223231}
    224232

  • trunk/BNC/newmat/newmat9.cpp

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

  • trunk/BNC/newmat/newmatap.h

    r2013 r8901  
    2626void QRZ(Matrix&, UpperTriangularMatrix&);
    2727
    28 void QRZ(const Matrix&, Matrix&, Matrix&);
     28void QRZ(const Matrix&, Matrix&, Matrix&);
     29
     30inline void QRZT(Matrix& X, Matrix& Y, LowerTriangularMatrix& L, Matrix& M)
     31   { QRZT(X, L); QRZT(X, Y, M); }
     32
     33inline void QRZ(Matrix& X, Matrix& Y, UpperTriangularMatrix& U, Matrix& M)
     34   { QRZ(X, U); QRZ(X, Y, M); }
    2935
    3036inline void HHDecompose(Matrix& X, LowerTriangularMatrix& L)
    … …  
    3844void updateQRZ(Matrix& X, UpperTriangularMatrix& U);
    3945
     46void updateQRZ(const Matrix& X, Matrix& MX, Matrix& MU);
     47
     48void updateQRZ(UpperTriangularMatrix& X, UpperTriangularMatrix& U);
     49
     50void updateQRZ(const UpperTriangularMatrix& X, Matrix& MX, Matrix& MU);
     51
    4052inline void UpdateQRZT(Matrix& X, LowerTriangularMatrix& L)
    4153   { updateQRZT(X, L); }
    … …  
    4355inline void UpdateQRZ(Matrix& X, UpperTriangularMatrix& U)
    4456   { updateQRZ(X, U); }
     57
     58inline void UpdateQRZ(UpperTriangularMatrix& X, UpperTriangularMatrix& U)
     59   { updateQRZ(X, U); }
     60
     61inline void UpdateQRZ(const UpperTriangularMatrix& X, Matrix& MX, Matrix& MU)
     62   { updateQRZ(X, MX, MU); }
     63   
     64inline void UpdateQRZ(const Matrix& X, Matrix& MX, Matrix& MU)
     65   { updateQRZ(X, MX, MU); }
     66
    4567
    4668// Matrix A's first n columns are orthonormal
    … …  
    5678
    5779
    58 // produces the Cholesky decomposition of A + x.t() * x where A = chol.t() * chol
    59 // and x is a RowVector
     80// produces the Cholesky decomposition of A + x.t() * x
     81// where A = chol.t() * chol and x is a RowVector
    6082void update_Cholesky(UpperTriangularMatrix& chol, RowVector x);
    6183inline void UpdateCholesky(UpperTriangularMatrix& chol, const RowVector& x)
    6284   { update_Cholesky(chol, x); }
    6385
    64 // produces the Cholesky decomposition of A - x.t() * x where A = chol.t() * chol
    65 // and x is a RowVector
     86// produces the Cholesky decomposition of A - x.t() * x
     87// where A = chol.t() * chol and x is a RowVector
    6688void downdate_Cholesky(UpperTriangularMatrix &chol, RowVector x);
    6789inline void DowndateCholesky(UpperTriangularMatrix &chol, const RowVector& x)

  • trunk/BNC/newmat/newmatio.h

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

  • trunk/BNC/newmat/precisio.h

    r2013 r8901  
    2525#endif
    2626
     27using namespace std;
     28       
    2729/// Floating point precision.
    2830class FloatingPointPrecision
    … …  
    3032public:
    3133   static int Dig()              // number of decimal digits or precision
    32       { return std::numeric_limits<Real>::digits10 ; }
     34      { return numeric_limits<Real>::digits10 ; }
    3335
    3436   static Real Epsilon()         // smallest number such that 1+Eps!=Eps
    35       { return std::numeric_limits<Real>::epsilon(); }
     37      { return numeric_limits<Real>::epsilon(); }
    3638
    3739   static int Mantissa()         // bits in mantisa
    38       { return std::numeric_limits<Real>::digits; }
     40      { return numeric_limits<Real>::digits; }
    3941
    4042   static Real Maximum()         // maximum value
    41       { return std::numeric_limits<Real>::max(); }
     43      { return numeric_limits<Real>::max(); }
    4244
    4345   static int MaximumDecimalExponent()  // maximum decimal exponent
    44       { return std::numeric_limits<Real>::max_exponent10; }
     46      { return numeric_limits<Real>::max_exponent10; }
    4547
    4648   static int MaximumExponent()  // maximum binary exponent
    47       { return std::numeric_limits<Real>::max_exponent; }
     49      { return numeric_limits<Real>::max_exponent; }
    4850
    4951   static Real LnMaximum()       // natural log of maximum
    … …  
    5153
    5254   static Real Minimum()         // minimum positive value
    53       { return std::numeric_limits<Real>::min(); }
     55      { return numeric_limits<Real>::min(); }
    5456
    5557   static int MinimumDecimalExponent() // minimum decimal exponent
    56       { return std::numeric_limits<Real>::min_exponent10; }
     58      { return numeric_limits<Real>::min_exponent10; }
    5759
    5860   static int MinimumExponent()  // minimum binary exponent
    59       { return std::numeric_limits<Real>::min_exponent; }
     61      { return numeric_limits<Real>::min_exponent; }
    6062
    6163   static Real LnMinimum()       // natural log of minimum
    … …  
    6365
    6466   static int Radix()            // exponent radix
    65       { return std::numeric_limits<Real>::radix; }
     67      { return numeric_limits<Real>::radix; }
    6668
    6769   static int Rounds()           // addition rounding (1 = does round)
    6870   {
    69           return std::numeric_limits<Real>::round_style ==
    70                  std::round_to_nearest ? 1 : 0;
     71          return numeric_limits<Real>::round_style ==
     72                 round_to_nearest ? 1 : 0;
    7173   }
    7274

  • trunk/BNC/newmat/submat.cpp

    r2013 r8901  
    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"));
    259262      MatrixRow mrx(gmx, LoadOnEntry);
    260263      MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
    … …  
    277280}
    278281
    279 void GetSubMatrix::operator-=(const BaseMatrix& bmx)
    280 {
    281    REPORT
    282    Tracer tr("SubMatrix(-=)"); GeneralMatrix* gmx = 0;
     282void GetSubMatrix::SP_eq(const BaseMatrix& bmx)
     283{
     284   REPORT
     285   Tracer tr("SubMatrix(SP_eq)"); GeneralMatrix* gmx = 0;
    283286   // MatrixConversionCheck mcc;         // Check for loss of info
    284287   Try
    … …  
    287290      if (row_number != gmx->Nrows() || col_number != gmx->Ncols())
    288291         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
     314void 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"));
    289327      MatrixRow mrx(gmx, LoadOnEntry);
    290328      MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);

  • trunk/BNC/src/PPP_SSR_I/pppFilter.cpp

    r8253 r8901  
    332332  xRec(3) = z();
    333333
    334   double rho0 = (satData->xx - xRec).norm_Frobenius();
     334  double rho0 = (satData->xx - xRec).NormFrobenius();
    335335  double dPhi = t_CST::omega * rho0 / t_CST::c;
    336336
    … …  
    341341  xRec += _tides->displacement(_time, xRec);
    342342
    343   satData->rho = (satData->xx - xRec).norm_Frobenius();
     343  satData->rho = (satData->xx - xRec).NormFrobenius();
    344344
    345345  double tropDelay = delay_saast(satData->eleSat) +
    … …  
    796796    // --------------------------------------
    797797    ColumnVector rho = rRec - rSat;
    798     rho /= rho.norm_Frobenius();
     798    rho /= rho.NormFrobenius();
    799799
    800800    // GPS Satellite unit Vectors sz, sy, sx
    801801    // -------------------------------------
    802     ColumnVector sz = -rSat / rSat.norm_Frobenius();
     802    ColumnVector sz = -rSat / rSat.NormFrobenius();
    803803
    804804    ColumnVector xSun = t_astro::Sun(Mjd);
    805     xSun /= xSun.norm_Frobenius();
     805    xSun /= xSun.NormFrobenius();
    806806
    807807    ColumnVector sy = crossproduct(sz, xSun);
    … …  
    839839    // ----------------
    840840    double alpha = DotProduct(dipSat,dipRec) /
    841                       (dipSat.norm_Frobenius() * dipRec.norm_Frobenius());
     841                      (dipSat.NormFrobenius() * dipRec.NormFrobenius());
    842842
    843843    if (alpha >  1.0) alpha =  1.0;
    … …  
    861861  Tracer tracer("t_pppFilter::cmpEle");
    862862  ColumnVector rr = satData->xx - _xcBanc.Rows(1,3);
    863   double       rho = rr.norm_Frobenius();
     863  double       rho = rr.NormFrobenius();
    864864
    865865  double neu[3];

  • trunk/BNC/src/combination/bnccomb.cpp

    r8736 r8901  
    11411141      }
    11421142      else {
    1143         double normMax = maxDiff[prn]->_diffRao.norm_Frobenius();
    1144         double norm    = corr->_diffRao.norm_Frobenius();
     1143        double normMax = maxDiff[prn]->_diffRao.NormFrobenius();
     1144        double norm    = corr->_diffRao.NormFrobenius();
    11451145        if (norm > normMax) {
    11461146          maxDiff[prn] = corr;
    … …  
    11651165      }
    11661166      else if (corr == maxDiff[prn]) {
    1167         double norm = corr->_diffRao.norm_Frobenius();
     1167        double norm = corr->_diffRao.NormFrobenius();
    11681168        if (norm > MAX_DISPLACEMENT) {
    11691169          out << _resTime.datestr().c_str()    << " "

  • trunk/BNC/src/rinex/reqcanalyze.cpp

    r8756 r8901  
    278278        ++nSatUsed;
    279279        ColumnVector dx = xSat.Rows(1,3) - xyzSta;
    280         double rho = dx.norm_Frobenius();
     280        double rho = dx.NormFrobenius();
    281281        AA(nSatUsed,1) = dx(1) / rho;
    282282        AA(nSatUsed,2) = dx(2) / rho;
Note: See TracChangeset for help on using the changeset viewer.