Index: /trunk/BNC/newmat/bandmat.cpp
===================================================================
--- /trunk/BNC/newmat/bandmat.cpp	(revision 8900)
+++ /trunk/BNC/newmat/bandmat.cpp	(revision 8901)
Index: /trunk/BNC/newmat/cholesky.cpp
===================================================================
--- /trunk/BNC/newmat/cholesky.cpp	(revision 8900)
+++ /trunk/BNC/newmat/cholesky.cpp	(revision 8901)
Index: /trunk/BNC/newmat/controlw.h
===================================================================
--- /trunk/BNC/newmat/controlw.h	(revision 8900)
+++ /trunk/BNC/newmat/controlw.h	(revision 8901)
Index: /trunk/BNC/newmat/evalue.cpp
===================================================================
--- /trunk/BNC/newmat/evalue.cpp	(revision 8900)
+++ /trunk/BNC/newmat/evalue.cpp	(revision 8901)
Index: /trunk/BNC/newmat/fft.cpp
===================================================================
--- /trunk/BNC/newmat/fft.cpp	(revision 8900)
+++ /trunk/BNC/newmat/fft.cpp	(revision 8901)
Index: /trunk/BNC/newmat/hholder.cpp
===================================================================
--- /trunk/BNC/newmat/hholder.cpp	(revision 8900)
+++ /trunk/BNC/newmat/hholder.cpp	(revision 8901)
@@ -335,4 +335,148 @@
 }
 
+// Following previous transformation,
+// now apply the same orthogonal transformation to (MX & MU)
+// Need the X Matrix but not the U.
+// Not optimised for accessing consecutive memory
+
+void updateQRZ(const Matrix& X, Matrix& MX, Matrix& MU)
+{
+   REPORT
+   Tracer et("updateQRZ(2)");
+   int s = X.Ncols(); int n = X.Nrows();
+   if (n != MX.Nrows())
+      Throw(ProgramException("Incompatible dimensions",X,MX));
+   if (s != MU.Nrows())
+      Throw(ProgramException("Incompatible dimensions",X,MU));
+   int t = MX.Ncols();
+   if (t != MU.Ncols())
+      Throw(ProgramException("Incompatible dimensions",MX,MU));
+   
+   if (s == 0) return;
+    
+   const Real* xi0 = X.data(); Real* mx = MX.data(); Real* muj = MU.data();
+   for (int i=1; i<=s; ++i)
+   {
+		  Real sum = 0.0;
+			{
+         const Real* xi=xi0; int k=n;
+			   while(k--) { sum += square(*xi); xi+= s;}
+			}
+      Real a0 = sqrt(2.0 - sum); Real* mxj0 = mx;
+      for (int j=1; j<=t; ++j)
+      {
+         Real sum = 0.0;
+         const Real* xi=xi0; Real* mxj=mxj0; int k=n; 
+         while(--k) { sum += *xi * *mxj; xi += s; mxj += t; }
+         sum += *xi * *mxj;    // last line of loop
+         sum += a0 * *muj;
+         xi=xi0; mxj=mxj0; k=n;
+         while(--k) { *mxj -= sum * *xi; xi += s; mxj += t; }
+         *mxj -= sum * *xi;    // last line of loop
+         *muj -= sum * a0; ++mxj0; ++muj;
+      }
+			++xi0;
+   }
+}
+
+
+
+// same thing as updateQRZ(Matrix& X, UpperTriangularMatrix& U)
+// except that X is upper triangular
+// contents of X are destroyed - results are in U
+// assume we can access efficiently by columns
+// e.g. X and U will fit in cache memory
+
+void updateQRZ(UpperTriangularMatrix& X, UpperTriangularMatrix& U)
+{
+   REPORT
+   Tracer et("updateQRZ(3)");
+   int s = X.Ncols();
+   if (s != U.Ncols())
+      Throw(ProgramException("Incompatible dimensions",X,U));
+   if (s == 0) return; 
+   Real* xi0 = X.data(); Real* u = U.data();
+   for (int i=1; i<=s; ++i)
+   {
+		  Real r = *u; Real sum = 0.0;
+			{
+         Real* xi=xi0; int k=i; int l=s;
+			   while(k--) { sum += square(*xi); xi+= --l;}
+			}
+      sum = sqrt(sum + square(r));
+      if (sum == 0.0) { REPORT X.column(i) = 0.0; *u = 0.0; }
+      else
+      {
+         Real frs = fabs(r) + sum;
+         Real a0 = sqrt(frs / sum); Real alpha = a0 / frs;
+         if (r <= 0) { REPORT *u = sum; alpha = -alpha; }
+         else { REPORT *u = -sum; }
+         {
+            Real* xj0=xi0; int k=i; int l=s;
+            while(k--) { *xj0 *= alpha; --l; xj0 += l;}
+         }
+         Real* xj0=xi0; Real* uj=u;
+         for (int j=i+1; j<=s; ++j)
+         {
+            Real sum = 0.0; ++xj0; ++uj;
+            Real* xi=xi0; Real* xj=xj0; int k=i; int l=s; 
+            while(k--) { sum += *xi * *xj; --l; xi += l; xj += l; }
+            sum += a0 * *uj;
+            xi=xi0; xj=xj0; k=i; l=s;
+            while(k--) { *xj -= sum * *xi; --l; xi += l; xj += l; }
+            *uj -= sum * a0;
+         }
+      }
+			++xi0; u += s-i+1;
+   }
+}
+
+// Following previous transformation,
+// now apply the same orthogonal transformation to (MX & MU)
+// Need the X UpperTriangularMatrix but not the U.
+
+void updateQRZ(const UpperTriangularMatrix& X, Matrix& MX, Matrix& MU)
+{
+   REPORT
+   Tracer et("updateQRZ(4)");
+   int s = X.Ncols();
+   if (s != MX.Nrows())
+      Throw(ProgramException("Incompatible dimensions",X,MX));
+   if (s != MU.Nrows())
+      Throw(ProgramException("Incompatible dimensions",X,MU));
+   int t = MX.Ncols();
+   if (t != MU.Ncols())
+      Throw(ProgramException("Incompatible dimensions",MX,MU));
+   if (s == 0) return;
+    
+   const Real* xi0 = X.data(); Real* mx = MX.data(); Real* muj = MU.data();
+   for (int i=1; i<=s; ++i)
+   {
+		  Real sum = 0.0;
+			{
+         const Real* xi=xi0; int k=i; int l=s;
+			   while(k--) { sum += square(*xi); xi+= --l;}
+			}
+      Real a0 = sqrt(2.0 - sum); Real* mxj0 = mx;
+      for (int j=1; j<=t; ++j)
+      {
+         Real sum = 0.0;
+         const Real* xi=xi0; Real* mxj=mxj0; int k=i; int l=s; 
+         while(--k) { sum += *xi * *mxj; --l; xi += l; mxj += t; }
+         sum += *xi * *mxj;    // last line of loop
+         sum += a0 * *muj;
+         xi=xi0; mxj=mxj0; k=i; l=s;
+         while(--k) { *mxj -= sum * *xi; --l; xi += l; mxj += t; }
+         *mxj -= sum * *xi;    // last line of loop
+         *muj -= sum * a0; ++mxj0; ++muj;
+      }
+			++xi0;
+   }
+}
+
+
+
+
+
 // Matrix A's first n columns are orthonormal
 // so A.Columns(1,n).t() * A.Columns(1,n) is the identity matrix.
Index: /trunk/BNC/newmat/include.h
===================================================================
--- /trunk/BNC/newmat/include.h	(revision 8900)
+++ /trunk/BNC/newmat/include.h	(revision 8901)
@@ -1,3 +1,3 @@
-/// \defgroup rbd_common RBD common library 
+/// \defgroup rbd_common RBD common library
 ///@{
 
@@ -44,5 +44,7 @@
 //#define HAS_INT64                     // if unsigned _int64 is recognised
                                         // used by newran03
-                                        
+
+//#define set_unix_options              // set if running UNIX or LINUX
+
 // comment out next line if Exception causes a problem
 #define TypeDefException
@@ -63,7 +65,19 @@
 #endif
 
+// for Intel C++ for Windows
+#if defined __ICL
+   #define _STANDARD_                   // use standard library
+   #define ios_format_flags ios::fmtflags
+#endif
+
 // for Microsoft Visual C++ 7 and above (and Intel simulating these)
 #if defined _MSC_VER && _MSC_VER >= 1300
    #define _STANDARD_                   // use standard library
+#endif
+
+// for Borland Builder C++ 2006 and above
+#if defined __BCPLUSPLUS__ && __BCPLUSPLUS__ >= 0x0570
+   #define _STANDARD_                   // use standard library
+   #define ios_format_flags ios::fmtflags
 #endif
 
@@ -90,6 +104,5 @@
       #include <fstream>
    #endif
-   ////using namespace std;
-   #define USE_STD_NAMESPACE
+   using namespace std;
 #else
 
Index: /trunk/BNC/newmat/jacobi.cpp
===================================================================
--- /trunk/BNC/newmat/jacobi.cpp	(revision 8900)
+++ /trunk/BNC/newmat/jacobi.cpp	(revision 8901)
Index: /trunk/BNC/newmat/myexcept.cpp
===================================================================
--- /trunk/BNC/newmat/myexcept.cpp	(revision 8900)
+++ /trunk/BNC/newmat/myexcept.cpp	(revision 8901)
@@ -26,7 +26,4 @@
 #endif
 
-#ifdef USE_STD_NAMESPACE
-using namespace std;
-#endif
 
 //#define REG_DEREG                    // for print out uses of new/delete
Index: /trunk/BNC/newmat/myexcept.h
===================================================================
--- /trunk/BNC/newmat/myexcept.h	(revision 8900)
+++ /trunk/BNC/newmat/myexcept.h	(revision 8901)
@@ -75,4 +75,5 @@
    static void AddTrace();               // insert trace in exception record
    static Tracer* last;                  // points to Tracer list
+   static void clear() {}                // for compatibility
    friend class BaseException;
 };
@@ -93,4 +94,5 @@
    static const char* what() { return what_error; }
                                          // for getting error message
+   static void clear() {}                // for compatibility
 };
 
Index: /trunk/BNC/newmat/newfft.cpp
===================================================================
--- /trunk/BNC/newmat/newfft.cpp	(revision 8900)
+++ /trunk/BNC/newmat/newfft.cpp	(revision 8901)
Index: /trunk/BNC/newmat/newmat.h
===================================================================
--- /trunk/BNC/newmat/newmat.h	(revision 8900)
+++ /trunk/BNC/newmat/newmat.h	(revision 8901)
@@ -447,5 +447,4 @@
 class GeneralMatrix : public BaseMatrix         // declarable matrix types
 {
-   virtual GeneralMatrix* Image() const;        // copy of matrix
 protected:
    int tag_val;                                 // shows whether can reuse
@@ -465,13 +464,8 @@
    void ReverseElements();                      // internal reverse of elements
    void ReverseElements(GeneralMatrix*);        // reverse order of elements
-   void operator=(Real);                        // set matrix to constant
    Real* GetStore();                            // get store or copy
    GeneralMatrix* BorrowStore(GeneralMatrix*, MatrixType);
                                                 // temporarily access store
    void GetMatrix(const GeneralMatrix*);        // used by = and initialise
-   void Eq(const BaseMatrix&, MatrixType);      // used by =
-   void Eq(const GeneralMatrix&);               // version with no conversion
-   void Eq(const BaseMatrix&, MatrixType, bool);// used by <<
-   void Eq2(const BaseMatrix&, MatrixType);     // cut down version of Eq
    int search(const BaseMatrix*) const;
    virtual GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
@@ -484,4 +478,5 @@
              // CleanUp when the data array has already been deleted
    void PlusEqual(const GeneralMatrix& gm);
+   void SP_Equal(const GeneralMatrix& gm);
    void MinusEqual(const GeneralMatrix& gm);
    void PlusEqual(Real f);
@@ -490,4 +485,8 @@
 public:
    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
+   void Eq(const BaseMatrix&, MatrixType);      // used by =
+   void Eq(const GeneralMatrix&);               // version with no conversion
+   void Eq(const BaseMatrix&, MatrixType, bool);// used by <<
+   void Eq2(const BaseMatrix&, MatrixType);     // cut down version of Eq
    virtual MatrixType type() const = 0;         // type of a matrix
    MatrixType Type() const { return type(); }
@@ -503,4 +502,5 @@
    const Real* data() const { return store; }
    const Real* const_data() const { return store; }
+   void operator=(Real);                        // set matrix to constant
    virtual ~GeneralMatrix();                    // delete store if set
    void tDelete();                              // delete if tag_val permits
@@ -526,4 +526,5 @@
    void Inject(const GeneralMatrix& GM) { inject(GM); }
    void operator+=(const BaseMatrix&);
+   void SP_eq(const BaseMatrix&);
    void operator-=(const BaseMatrix&);
    void operator*=(const BaseMatrix&);
@@ -579,4 +580,5 @@
 //   ReturnMatrix Reverse() const;                // reverse order of elements
    void cleanup();                              // to clear store
+   virtual GeneralMatrix* Image() const;        // copy of matrix
 
    friend class Matrix;
@@ -625,5 +627,4 @@
 class Matrix : public GeneralMatrix
 {
-   GeneralMatrix* Image() const;                // copy of matrix
 public:
    Matrix() {}
@@ -668,8 +669,13 @@
    Real minimum2(int& i, int& j) const;
    void operator+=(const Matrix& M) { PlusEqual(M); }
+   void SP_eq(const Matrix& M) { SP_Equal(M); }
    void operator-=(const Matrix& M) { MinusEqual(M); }
+   void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
+   void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
+   void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    void operator+=(Real f) { GeneralMatrix::Add(f); }
    void operator-=(Real f) { GeneralMatrix::Add(-f); }
    void swap(Matrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
+   GeneralMatrix* Image() const;                // copy of matrix
    friend Real dotproduct(const Matrix& A, const Matrix& B);
    NEW_DELETE(Matrix)
@@ -679,5 +685,4 @@
 class SquareMatrix : public Matrix
 {
-   GeneralMatrix* Image() const;                // copy of matrix
 public:
    SquareMatrix() {}
@@ -701,8 +706,13 @@
    void ReSize(const GeneralMatrix& A) { resize(A); }
    void operator+=(const Matrix& M) { PlusEqual(M); }
+   void SP_eq(const Matrix& M) { SP_Equal(M); }
    void operator-=(const Matrix& M) { MinusEqual(M); }
+   void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
+   void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
+   void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    void operator+=(Real f) { GeneralMatrix::Add(f); }
    void operator-=(Real f) { GeneralMatrix::Add(-f); }
    void swap(SquareMatrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
+   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(SquareMatrix)
 };
@@ -711,5 +721,4 @@
 class nricMatrix : public Matrix
 {
-   GeneralMatrix* Image() const;                // copy of matrix
    Real** row_pointer;                          // points to rows
    void MakeRowPointer();                       // build rowpointer
@@ -743,8 +752,13 @@
    void MiniCleanUp();
    void operator+=(const Matrix& M) { PlusEqual(M); }
+   void SP_eq(const Matrix& M) { SP_Equal(M); }
    void operator-=(const Matrix& M) { MinusEqual(M); }
+   void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
+   void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
+   void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    void operator+=(Real f) { GeneralMatrix::Add(f); }
    void operator-=(Real f) { GeneralMatrix::Add(-f); }
    void swap(nricMatrix& gm);
+   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(nricMatrix)
 };
@@ -753,5 +767,4 @@
 class SymmetricMatrix : public GeneralMatrix
 {
-   GeneralMatrix* Image() const;                // copy of matrix
 public:
    SymmetricMatrix() {}
@@ -789,8 +802,13 @@
    void ReSize(const GeneralMatrix& A) { resize(A); }
    void operator+=(const SymmetricMatrix& M) { PlusEqual(M); }
+   void SP_eq(const SymmetricMatrix& M) { SP_Equal(M); }
    void operator-=(const SymmetricMatrix& M) { MinusEqual(M); }
+   void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
+   void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
+   void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    void operator+=(Real f) { GeneralMatrix::Add(f); }
    void operator-=(Real f) { GeneralMatrix::Add(-f); }
    void swap(SymmetricMatrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
+   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(SymmetricMatrix)
 };
@@ -799,5 +817,4 @@
 class UpperTriangularMatrix : public GeneralMatrix
 {
-   GeneralMatrix* Image() const;                // copy of matrix
 public:
    UpperTriangularMatrix() {}
@@ -837,9 +854,14 @@
    MatrixBandWidth bandwidth() const;
    void operator+=(const UpperTriangularMatrix& M) { PlusEqual(M); }
+   void SP_eq(const UpperTriangularMatrix& M) { SP_Equal(M); }
    void operator-=(const UpperTriangularMatrix& M) { MinusEqual(M); }
+   void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
+   void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
+   void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    void operator+=(Real f) { GeneralMatrix::operator+=(f); }
    void operator-=(Real f) { GeneralMatrix::operator-=(f); }
    void swap(UpperTriangularMatrix& gm)
       { GeneralMatrix::swap((GeneralMatrix&)gm); }
+   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(UpperTriangularMatrix)
 };
@@ -848,5 +870,4 @@
 class LowerTriangularMatrix : public GeneralMatrix
 {
-   GeneralMatrix* Image() const;                // copy of matrix
 public:
    LowerTriangularMatrix() {}
@@ -885,9 +906,14 @@
    MatrixBandWidth bandwidth() const;
    void operator+=(const LowerTriangularMatrix& M) { PlusEqual(M); }
+   void SP_eq(const LowerTriangularMatrix& M) { SP_Equal(M); }
    void operator-=(const LowerTriangularMatrix& M) { MinusEqual(M); }
+   void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
+   void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
+   void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    void operator+=(Real f) { GeneralMatrix::operator+=(f); }
    void operator-=(Real f) { GeneralMatrix::operator-=(f); }
    void swap(LowerTriangularMatrix& gm)
       { GeneralMatrix::swap((GeneralMatrix&)gm); }
+   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(LowerTriangularMatrix)
 };
@@ -896,5 +922,4 @@
 class DiagonalMatrix : public GeneralMatrix
 {
-   GeneralMatrix* Image() const;                // copy of matrix
 public:
    DiagonalMatrix() {}
@@ -942,9 +967,14 @@
 //   ReturnMatrix Reverse() const;                // reverse order of elements
    void operator+=(const DiagonalMatrix& M) { PlusEqual(M); }
+   void SP_eq(const DiagonalMatrix& M) { SP_Equal(M); }
    void operator-=(const DiagonalMatrix& M) { MinusEqual(M); }
+   void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
+   void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
+   void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    void operator+=(Real f) { GeneralMatrix::operator+=(f); }
    void operator-=(Real f) { GeneralMatrix::operator-=(f); }
    void swap(DiagonalMatrix& gm)
       { GeneralMatrix::swap((GeneralMatrix&)gm); }
+   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(DiagonalMatrix)
 };
@@ -953,5 +983,4 @@
 class RowVector : public Matrix
 {
-   GeneralMatrix* Image() const;                // copy of matrix
 public:
    RowVector() { nrows_val = 1; }
@@ -997,9 +1026,14 @@
    // friend ReturnMatrix GetMatrixRow(Matrix& A, int row);
    void operator+=(const Matrix& M) { PlusEqual(M); }
+   void SP_eq(const Matrix& M) { SP_Equal(M); }
    void operator-=(const Matrix& M) { MinusEqual(M); }
+   void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
+   void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
+   void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    void operator+=(Real f) { GeneralMatrix::Add(f); }
    void operator-=(Real f) { GeneralMatrix::Add(-f); }
    void swap(RowVector& gm)
       { GeneralMatrix::swap((GeneralMatrix&)gm); }
+   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(RowVector)
 };
@@ -1008,5 +1042,4 @@
 class ColumnVector : public Matrix
 {
-   GeneralMatrix* Image() const;                // copy of matrix
 public:
    ColumnVector() { ncols_val = 1; }
@@ -1046,9 +1079,14 @@
 //   ReturnMatrix Reverse() const;                // reverse order of elements
    void operator+=(const Matrix& M) { PlusEqual(M); }
+   void SP_eq(const Matrix& M) { SP_Equal(M); }
    void operator-=(const Matrix& M) { MinusEqual(M); }
+   void operator+=(const BaseMatrix& M) { GeneralMatrix::operator+=(M); }
+   void SP_eq(const BaseMatrix& M) { GeneralMatrix::SP_eq(M); }
+   void operator-=(const BaseMatrix& M) { GeneralMatrix::operator-=(M); }
    void operator+=(Real f) { GeneralMatrix::Add(f); }
    void operator-=(Real f) { GeneralMatrix::Add(-f); }
    void swap(ColumnVector& gm)
       { GeneralMatrix::swap((GeneralMatrix&)gm); }
+   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(ColumnVector)
 };
@@ -1064,5 +1102,4 @@
    void ludcmp();
    void get_aux(CroutMatrix&);                  // for copying indx[] etc
-   GeneralMatrix* Image() const;                // copy of matrix
 public:
    CroutMatrix(const BaseMatrix&);
@@ -1088,4 +1125,5 @@
    bool even_exchanges() const { return d; }
    void swap(CroutMatrix& gm);
+   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(CroutMatrix)
 };
@@ -1096,5 +1134,4 @@
 class BandMatrix : public GeneralMatrix
 {
-   GeneralMatrix* Image() const;                // copy of matrix
 protected:
    void CornerClear() const;                    // set unused elements to zero
@@ -1161,4 +1198,5 @@
    void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
    void swap(BandMatrix& gm);
+   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(BandMatrix)
 };
@@ -1167,5 +1205,4 @@
 class UpperBandMatrix : public BandMatrix
 {
-   GeneralMatrix* Image() const;                // copy of matrix
 public:
    UpperBandMatrix() {}
@@ -1200,4 +1237,5 @@
    void swap(UpperBandMatrix& gm)
       { BandMatrix::swap((BandMatrix&)gm); }
+   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(UpperBandMatrix)
 };
@@ -1206,5 +1244,4 @@
 class LowerBandMatrix : public BandMatrix
 {
-   GeneralMatrix* Image() const;                // copy of matrix
 public:
    LowerBandMatrix() {}
@@ -1239,4 +1276,5 @@
    void swap(LowerBandMatrix& gm)
       { BandMatrix::swap((BandMatrix&)gm); }
+   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(LowerBandMatrix)
 };
@@ -1245,5 +1283,4 @@
 class SymmetricBandMatrix : public GeneralMatrix
 {
-   GeneralMatrix* Image() const;                // copy of matrix
    void CornerClear() const;                    // set unused elements to zero
    short SimpleAddOK(const GeneralMatrix* gm);
@@ -1300,4 +1337,5 @@
    void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
    void swap(SymmetricBandMatrix& gm);
+   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(SymmetricBandMatrix)
 };
@@ -1314,5 +1352,4 @@
    void ludcmp();
    void get_aux(BandLUMatrix&);                 // for copying indx[] etc
-   GeneralMatrix* Image() const;                // copy of matrix
 public:
    BandLUMatrix(const BaseMatrix&);
@@ -1342,4 +1379,5 @@
    MatrixBandWidth bandwidth() const;
    void swap(BandLUMatrix& gm);
+   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(BandLUMatrix)
 };
@@ -1350,5 +1388,4 @@
 class IdentityMatrix : public GeneralMatrix
 {
-   GeneralMatrix* Image() const;          // copy of matrix
 public:
    IdentityMatrix() {}
@@ -1386,4 +1423,5 @@
    void swap(IdentityMatrix& gm)
       { GeneralMatrix::swap((GeneralMatrix&)gm); }
+   GeneralMatrix* Image() const;          // copy of matrix
    NEW_DELETE(IdentityMatrix)
 };
@@ -1409,4 +1447,5 @@
    void operator=(const BaseMatrix&);
    void operator+=(const BaseMatrix&);
+   void SP_eq(const BaseMatrix&);
    void operator-=(const BaseMatrix&);
    void operator*=(const BaseMatrix&);
@@ -1781,4 +1820,5 @@
    void operator=(const BaseMatrix&);
    void operator+=(const BaseMatrix&);
+   void SP_eq(const BaseMatrix&);
    void operator-=(const BaseMatrix&);
    void operator=(const GetSubMatrix& m) { operator=((const BaseMatrix&)m); }
Index: /trunk/BNC/newmat/newmat.pro
===================================================================
--- /trunk/BNC/newmat/newmat.pro	(revision 8900)
+++ /trunk/BNC/newmat/newmat.pro	(revision 8901)
@@ -10,15 +10,15 @@
 MOC_DIR     = .moc
 
-HEADERS += controlw.h include.h myexcept.h  \
-           newmatap.h newmat.h newmatio.h   \
-           newmatrc.h newmatrm.h precisio.h
+HEADERS += boolean.h controlw.h include.h myexcept.h  \
+           newmatap.h newmat.h newmatio.h newmatnl.h  \
+           newmatrc.h newmatrm.h precisio.h solution.h
 
-SOURCES += bandmat.cpp cholesky.cpp evalue.cpp  \
-           fft.cpp hholder.cpp jacobi.cpp       \
-           myexcept.cpp newfft.cpp newmat1.cpp  \
-           newmat2.cpp newmat3.cpp newmat4.cpp  \
-           newmat5.cpp newmat6.cpp newmat7.cpp  \
-           newmat8.cpp newmat9.cpp newmatex.cpp \
-           newmatrm.cpp nm_misc.cpp sort.cpp    \
-           submat.cpp svd.cpp
+SOURCES += bandmat.cpp cholesky.cpp evalue.cpp     \
+           fft.cpp hholder.cpp jacobi.cpp          \
+           myexcept.cpp newfft.cpp newmat1.cpp     \
+           newmat2.cpp newmat3.cpp newmat4.cpp     \
+           newmat5.cpp newmat6.cpp newmat7.cpp     \
+           newmat8.cpp newmat9.cpp newmatex.cpp    \
+           newmatrm.cpp newmatnl.cpp nm_misc.cpp   \
+	   sort.cpp submat.cpp svd.cpp solution.cpp
 
Index: /trunk/BNC/newmat/newmat1.cpp
===================================================================
--- /trunk/BNC/newmat/newmat1.cpp	(revision 8900)
+++ /trunk/BNC/newmat/newmat1.cpp	(revision 8901)
Index: /trunk/BNC/newmat/newmat2.cpp
===================================================================
--- /trunk/BNC/newmat/newmat2.cpp	(revision 8900)
+++ /trunk/BNC/newmat/newmat2.cpp	(revision 8901)
Index: /trunk/BNC/newmat/newmat3.cpp
===================================================================
--- /trunk/BNC/newmat/newmat3.cpp	(revision 8900)
+++ /trunk/BNC/newmat/newmat3.cpp	(revision 8901)
Index: /trunk/BNC/newmat/newmat4.cpp
===================================================================
--- /trunk/BNC/newmat/newmat4.cpp	(revision 8900)
+++ /trunk/BNC/newmat/newmat4.cpp	(revision 8901)
Index: /trunk/BNC/newmat/newmat5.cpp
===================================================================
--- /trunk/BNC/newmat/newmat5.cpp	(revision 8900)
+++ /trunk/BNC/newmat/newmat5.cpp	(revision 8901)
Index: /trunk/BNC/newmat/newmat6.cpp
===================================================================
--- /trunk/BNC/newmat/newmat6.cpp	(revision 8900)
+++ /trunk/BNC/newmat/newmat6.cpp	(revision 8901)
@@ -495,10 +495,23 @@
    Tracer tr("GeneralMatrix::operator+=");
    // MatrixConversionCheck mcc;
-   Protect();                                   // so it cannot get deleted
-						// during Evaluate
+   Protect();  // so it cannot get deleted during Evaluate
    GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
    AddedMatrix am(this,gm);
    if (gm==this) Release(2); else Release();
    Eq2(am,type());
+}
+
+// GeneralMatrix operators
+
+void GeneralMatrix::SP_eq(const BaseMatrix& X)
+{
+   REPORT
+   Tracer tr("GeneralMatrix::SP_eq");
+   // MatrixConversionCheck mcc;
+   Protect();  // so it cannot get deleted during Evaluate
+   GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
+   SPMatrix spm(this,gm);
+   if (gm==this) Release(2); else Release();
+   Eq2(spm,type());
 }
 
@@ -586,4 +599,19 @@
    if (gmx==gm) gm->Release(2); else gm->Release();
    GeneralMatrix* gmy = am.Evaluate();
+   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
+   else { REPORT }
+   gm->Protect();
+}
+
+void GenericMatrix::SP_eq(const BaseMatrix& X)
+{
+   REPORT
+   Tracer tr("GenericMatrix::SP_eq");
+   if (!gm) Throw(ProgramException("GenericMatrix is null"));
+   gm->Protect();            // so it cannot get deleted during Evaluate
+   GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
+   SPMatrix spm(gm,gmx);
+   if (gmx==gm) gm->Release(2); else gm->Release();
+   GeneralMatrix* gmy = spm.Evaluate();
    if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
    else { REPORT }
Index: /trunk/BNC/newmat/newmat7.cpp
===================================================================
--- /trunk/BNC/newmat/newmat7.cpp	(revision 8900)
+++ /trunk/BNC/newmat/newmat7.cpp	(revision 8901)
@@ -214,11 +214,19 @@
 }
 
-static void SP(GeneralMatrix* gm, GeneralMatrix* gm2)
-{
-   REPORT
-   Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
+static void SP(GeneralMatrix* gm, const GeneralMatrix* gm2)
+{
+   REPORT
+   const Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
    while (i--)
    { *s++ *= *s2++; *s++ *= *s2++; *s++ *= *s2++; *s++ *= *s2++; }
    i=gm->Storage() & 3; while (i--) *s++ *= *s2++;
+}
+
+void GeneralMatrix::SP_Equal(const GeneralMatrix& gm)
+{
+   REPORT
+   if (nrows_val != gm.nrows_val || ncols_val != gm.ncols_val)
+      Throw(IncompatibleDimensionsException(*this, gm));
+   SP(this, &gm);
 }
 
Index: /trunk/BNC/newmat/newmat8.cpp
===================================================================
--- /trunk/BNC/newmat/newmat8.cpp	(revision 8900)
+++ /trunk/BNC/newmat/newmat8.cpp	(revision 8901)
Index: /trunk/BNC/newmat/newmat9.cpp
===================================================================
--- /trunk/BNC/newmat/newmat9.cpp	(revision 8900)
+++ /trunk/BNC/newmat/newmat9.cpp	(revision 8901)
@@ -44,5 +44,6 @@
    MatrixRow mr((GeneralMatrix*)&X, LoadOnEntry);
    int w = s.width();  int nr = X.Nrows();  ios_format_flags f = s.flags();
-   s.setf(ios::fixed, ios::floatfield);
+   if (f & ios::scientific) s.setf(ios::scientific, ios::floatfield);
+   else s.setf(ios::fixed, ios::floatfield);
    for (int i=1; i<=nr; i++)
    {
Index: /trunk/BNC/newmat/newmatap.h
===================================================================
--- /trunk/BNC/newmat/newmatap.h	(revision 8900)
+++ /trunk/BNC/newmat/newmatap.h	(revision 8901)
@@ -26,5 +26,11 @@
 void QRZ(Matrix&, UpperTriangularMatrix&);
 
-void QRZ(const Matrix&, Matrix&, Matrix&);
+void QRZ(const Matrix&, Matrix&, Matrix&); 
+
+inline void QRZT(Matrix& X, Matrix& Y, LowerTriangularMatrix& L, Matrix& M)
+   { QRZT(X, L); QRZT(X, Y, M); } 
+
+inline void QRZ(Matrix& X, Matrix& Y, UpperTriangularMatrix& U, Matrix& M)
+   { QRZ(X, U); QRZ(X, Y, M); } 
 
 inline void HHDecompose(Matrix& X, LowerTriangularMatrix& L)
@@ -38,4 +44,10 @@
 void updateQRZ(Matrix& X, UpperTriangularMatrix& U);
 
+void updateQRZ(const Matrix& X, Matrix& MX, Matrix& MU);
+
+void updateQRZ(UpperTriangularMatrix& X, UpperTriangularMatrix& U);
+
+void updateQRZ(const UpperTriangularMatrix& X, Matrix& MX, Matrix& MU);
+
 inline void UpdateQRZT(Matrix& X, LowerTriangularMatrix& L)
    { updateQRZT(X, L); }
@@ -43,4 +55,14 @@
 inline void UpdateQRZ(Matrix& X, UpperTriangularMatrix& U)
    { updateQRZ(X, U); }
+
+inline void UpdateQRZ(UpperTriangularMatrix& X, UpperTriangularMatrix& U)
+   { updateQRZ(X, U); }
+
+inline void UpdateQRZ(const UpperTriangularMatrix& X, Matrix& MX, Matrix& MU)
+   { updateQRZ(X, MX, MU); }
+   
+inline void UpdateQRZ(const Matrix& X, Matrix& MX, Matrix& MU)
+   { updateQRZ(X, MX, MU); }
+
 
 // Matrix A's first n columns are orthonormal
@@ -56,12 +78,12 @@
 
 
-// produces the Cholesky decomposition of A + x.t() * x where A = chol.t() * chol
-// and x is a RowVector
+// produces the Cholesky decomposition of A + x.t() * x
+// where A = chol.t() * chol and x is a RowVector
 void update_Cholesky(UpperTriangularMatrix& chol, RowVector x);
 inline void UpdateCholesky(UpperTriangularMatrix& chol, const RowVector& x)
    { update_Cholesky(chol, x); }
 
-// produces the Cholesky decomposition of A - x.t() * x where A = chol.t() * chol
-// and x is a RowVector
+// produces the Cholesky decomposition of A - x.t() * x
+// where A = chol.t() * chol and x is a RowVector
 void downdate_Cholesky(UpperTriangularMatrix &chol, RowVector x);
 inline void DowndateCholesky(UpperTriangularMatrix &chol, const RowVector& x)
Index: /trunk/BNC/newmat/newmatex.cpp
===================================================================
--- /trunk/BNC/newmat/newmatex.cpp	(revision 8900)
+++ /trunk/BNC/newmat/newmatex.cpp	(revision 8901)
Index: /trunk/BNC/newmat/newmatio.h
===================================================================
--- /trunk/BNC/newmat/newmatio.h	(revision 8900)
+++ /trunk/BNC/newmat/newmatio.h	(revision 8901)
@@ -16,7 +16,9 @@
 #include "newmat.h"
 
-#ifdef USE_STD_NAMESPACE
-using namespace std;
+#ifdef use_namespace
+namespace NEWMAT {
 #endif
+
+
 
 // **************************** input/output *****************************/
Index: /trunk/BNC/newmat/newmatnl.cpp
===================================================================
--- /trunk/BNC/newmat/newmatnl.cpp	(revision 8901)
+++ /trunk/BNC/newmat/newmatnl.cpp	(revision 8901)
@@ -0,0 +1,266 @@
+/// \ingroup newmat
+///@{
+
+/// \file newmatnl.cpp
+/// Non-linear optimisation.
+
+// Copyright (C) 1993,4,5,6: R B Davies
+
+
+#define WANT_MATH
+#define WANT_STREAM
+
+#include "newmatap.h"
+#include "newmatnl.h"
+
+#ifdef use_namespace
+namespace NEWMAT {
+#endif
+
+
+
+void FindMaximum2::Fit(ColumnVector& Theta, int n_it)
+{
+   Tracer tr("FindMaximum2::Fit");
+   enum State {Start, Restart, Continue, Interpolate, Extrapolate,
+      Fail, Convergence};
+   State TheState = Start;
+   Real z,w,x,x2,g,l1,l2,l3,d1,d2=0,d3;
+   ColumnVector Theta1, Theta2, Theta3;
+   int np = Theta.Nrows();
+   ColumnVector H1(np), H3, HP(np), K, K1(np);
+   bool oorg, conv;
+   int counter = 0;
+   Theta1 = Theta; HP = 0.0; g = 0.0;
+
+   // This is really a set of gotos and labels, but they do not work
+   // correctly in AT&T C++ and Sun 4.01 C++.
+
+   for(;;)
+   {
+      switch (TheState)
+      {
+      case Start:
+	 tr.ReName("FindMaximum2::Fit/Start");
+	 Value(Theta1, true, l1, oorg);
+	 if (oorg) Throw(ProgramException("invalid starting value\n"));
+
+      case Restart:
+	 tr.ReName("FindMaximum2::Fit/ReStart");
+	 conv = NextPoint(H1, d1);
+	 if (conv) { TheState = Convergence; break; }
+	 if (counter++ > n_it) { TheState = Fail; break; }
+
+	 z = 1.0 / sqrt(d1);
+	 H3 = H1 * z; K = (H3 - HP) * g; HP = H3;
+	 g = 0.0;                     // de-activate to use curved projection
+	 if ( g == 0.0 ) K1 = 0.0; else K1 = K * 0.2 + K1 * 0.6;
+	 // (K - K1) * alpha + K1 * (1 - alpha)
+	 //     = K * alpha + K1 * (1 - 2 * alpha)
+	 K = K1 * d1; g = z;
+
+      case Continue:
+	 tr.ReName("FindMaximum2::Fit/Continue");
+	 Theta2 = Theta1 + H1 + K;
+	 Value(Theta2, false, l2, oorg);
+	 if (counter++ > n_it) { TheState = Fail; break; }
+	 if (oorg)
+	 {
+	    H1 *= 0.5; K *= 0.25; d1 *= 0.5; g *= 2.0;
+	    TheState =  Continue; break;
+	 }
+	 d2 = LastDerivative(H1 + K * 2.0);
+
+      case Interpolate:
+	 tr.ReName("FindMaximum2::Fit/Interpolate");
+	 z = d1 + d2 - 3.0 * (l2 - l1);
+	 w = z * z - d1 * d2;
+	 if (w < 0.0) { TheState = Extrapolate; break; }
+	 w = z + sqrt(w);
+	 if (1.5 * w + d1 < 0.0)
+	    { TheState = Extrapolate; break; }
+	 if (d2 > 0.0 && l2 > l1 && w > 0.0)
+	    { TheState = Extrapolate; break; }
+	 x = d1 / (w + d1); x2 = x * x; g /= x;
+	 Theta3 = Theta1 + H1 * x + K * x2;
+	 Value(Theta3, true, l3, oorg);
+	 if (counter++ > n_it) { TheState = Fail; break; }
+	 if (oorg)
+	 {
+	    if (x <= 1.0)
+	       { x *= 0.5; x2 = x*x; g *= 2.0; d1 *= x; H1 *= x; K *= x2; }
+	    else
+	    {
+	       x = 0.5 * (x-1.0); x2 = x*x; Theta1 = Theta2;
+	       H1 = (H1 + K * 2.0) * x;
+	       K *= x2; g = 0.0; d1 = x * d2; l1 = l2;
+	    }
+	    TheState = Continue; break;
+	 }
+
+	 if (l3 >= l1 && l3 >= l2)
+	    { Theta1 = Theta3; l1 = l3; TheState =  Restart; break; }
+
+	 d3 = LastDerivative(H1 + K * 2.0);
+	 if (l1 > l2)
+	    { H1 *= x; K *= x2; Theta2 = Theta3; d1 *= x; d2 = d3*x; }
+	 else
+	 {
+	    Theta1 = Theta2; Theta2 = Theta3;
+	    x -= 1.0; x2 = x*x; g = 0.0; H1 = (H1 + K * 2.0) * x;
+	    K *= x2; l1 = l2; l2 = l3; d1 = x*d2; d2 = x*d3;
+	    if (d1 <= 0.0) { TheState = Start; break; }
+	 }
+	 TheState =  Interpolate; break;
+
+      case Extrapolate:
+	 tr.ReName("FindMaximum2::Fit/Extrapolate");
+	 Theta1 = Theta2; g = 0.0; K *= 4.0; H1 = (H1 * 2.0 + K);
+	 d1 = 2.0 * d2; l1 = l2;
+	 TheState = Continue; break;
+
+      case Fail:
+	 Throw(ConvergenceException(Theta));
+
+      case Convergence:
+	 Theta = Theta1; return;
+      }
+   }
+}
+
+
+
+void NonLinearLeastSquares::Value
+   (const ColumnVector& Parameters, bool, Real& v, bool& oorg)
+{
+   Tracer tr("NonLinearLeastSquares::Value");
+   Y.resize(n_obs); X.resize(n_obs,n_param);
+   // put the fitted values in Y, the derivatives in X.
+   Pred.Set(Parameters);
+   if (!Pred.IsValid()) { oorg=true; return; }
+   for (int i=1; i<=n_obs; i++)
+   {
+      Y(i) = Pred(i);
+      X.Row(i) = Pred.Derivatives();
+   }
+   if (!Pred.IsValid()) { oorg=true; return; }  // check afterwards as well
+   Y = *DataPointer - Y; Real ssq = Y.SumSquare();
+   errorvar =  ssq / (n_obs - n_param);
+   cout << endl;
+   cout << setw(15) << setprecision(10) << " " << errorvar;
+   Derivs = Y.t() * X;          // get the derivative and stash it
+   oorg = false; v = -0.5 * ssq;
+}
+
+bool NonLinearLeastSquares::NextPoint(ColumnVector& Adj, Real& test)
+{
+   Tracer tr("NonLinearLeastSquares::NextPoint");
+   QRZ(X, U); QRZ(X, Y, M);     // do the QR decomposition
+   test = M.SumSquare();
+   cout << " " << setw(15) << setprecision(10)
+      << test << " " << Y.SumSquare() / (n_obs - n_param);
+   Adj = U.i() * M;
+   if (test < errorvar * criterion) return true;
+   else return false;
+}
+
+Real NonLinearLeastSquares::LastDerivative(const ColumnVector& H)
+{ return (Derivs * H).AsScalar(); }
+
+void NonLinearLeastSquares::Fit(const ColumnVector& Data,
+   ColumnVector& Parameters)
+{
+   Tracer tr("NonLinearLeastSquares::Fit");
+   n_param = Parameters.Nrows(); n_obs = Data.Nrows();
+   DataPointer = &Data;
+   FindMaximum2::Fit(Parameters, Lim);
+   cout << "\nConverged" << endl;
+}
+
+void NonLinearLeastSquares::MakeCovariance()
+{
+   if (Covariance.Nrows()==0)
+   {
+      UpperTriangularMatrix UI = U.i();
+      Covariance << UI * UI.t() * errorvar;
+      SE << Covariance;                 // get diagonals
+      for (int i = 1; i<=n_param; i++) SE(i) = sqrt(SE(i));
+   }
+}
+
+void NonLinearLeastSquares::GetStandardErrors(ColumnVector& SEX)
+   { MakeCovariance(); SEX = SE.AsColumn(); }
+
+void NonLinearLeastSquares::GetCorrelations(SymmetricMatrix& Corr)
+   { MakeCovariance(); Corr << SE.i() * Covariance * SE.i(); }
+
+void NonLinearLeastSquares::GetHatDiagonal(DiagonalMatrix& Hat) const
+{
+   Hat.resize(n_obs);
+   for (int i = 1; i<=n_obs; i++) Hat(i) = X.Row(i).SumSquare();
+}
+
+
+// the MLE_D_FI routines
+
+void MLE_D_FI::Value
+   (const ColumnVector& Parameters, bool wg, Real& v, bool& oorg)
+{
+   Tracer tr("MLE_D_FI::Value");
+   if (!LL.IsValid(Parameters,wg)) { oorg=true; return; }
+   v = LL.LogLikelihood();
+   if (!LL.IsValid()) { oorg=true; return; }     // check validity again
+   cout << endl;
+   cout << setw(20) << setprecision(10) << v;
+   oorg = false;
+   Derivs = LL.Derivatives();                    // Get derivatives
+}
+
+bool MLE_D_FI::NextPoint(ColumnVector& Adj, Real& test)
+{
+   Tracer tr("MLE_D_FI::NextPoint");
+   SymmetricMatrix FI = LL.FI();
+   LT = Cholesky(FI);
+   ColumnVector Adj1 = LT.i() * Derivs;
+   Adj = LT.t().i() * Adj1;
+   test = SumSquare(Adj1);
+   cout << "   " << setw(20) << setprecision(10) << test;
+   return (test < Criterion);
+}
+
+Real MLE_D_FI::LastDerivative(const ColumnVector& H)
+{ return (Derivs.t() * H).AsScalar(); }
+
+void MLE_D_FI::Fit(ColumnVector& Parameters)
+{
+   Tracer tr("MLE_D_FI::Fit");
+   FindMaximum2::Fit(Parameters,Lim);
+   cout << "\nConverged" << endl;
+}
+  
+void MLE_D_FI::MakeCovariance()
+{
+   if (Covariance.Nrows()==0)
+   {
+      LowerTriangularMatrix LTI = LT.i();
+      Covariance << LTI.t() * LTI;
+      SE << Covariance;                // get diagonal
+      int n = Covariance.Nrows();
+      for (int i=1; i <= n; i++) SE(i) = sqrt(SE(i));
+   }
+}
+
+void MLE_D_FI::GetStandardErrors(ColumnVector& SEX)
+{ MakeCovariance(); SEX = SE.AsColumn(); }
+   
+void MLE_D_FI::GetCorrelations(SymmetricMatrix& Corr)
+{ MakeCovariance(); Corr << SE.i() * Covariance * SE.i(); }
+
+
+
+#ifdef use_namespace
+}
+#endif
+
+
+///@}
Index: /trunk/BNC/newmat/newmatnl.h
===================================================================
--- /trunk/BNC/newmat/newmatnl.h	(revision 8901)
+++ /trunk/BNC/newmat/newmatnl.h	(revision 8901)
@@ -0,0 +1,327 @@
+/// \ingroup newmat
+///@{
+
+/// \file newmatnl.h
+/// Header file for non-linear optimisation
+
+// Copyright (C) 1993,4,5: R B Davies
+
+#ifndef NEWMATNL_LIB
+#define NEWMATNL_LIB 0
+
+#include "newmat.h"
+
+#ifdef use_namespace
+namespace NEWMAT {
+#endif
+
+
+
+/*
+
+This is a beginning of a series of classes for non-linear optimisation.
+
+At present there are two classes. FindMaximum2 is the basic optimisation
+strategy when one is doing an optimisation where one has first
+derivatives and estimates of the second derivatives. Class
+NonLinearLeastSquares is derived from FindMaximum2. This provides the
+functions that calculate function values and derivatives.
+
+A third class is now added. This is for doing maximum-likelihood when
+you have first derviatives and something like the Fisher Information
+matrix (eg the variance covariance matrix of the first derivatives or
+minus the second derivatives - this matrix is assumed to be positive
+definite).
+
+
+
+   class FindMaximum2
+
+Suppose T is the ColumnVector of parameters, F(T) the function we want
+to maximise, D(T) the ColumnVector of derivatives of F with respect to
+T, and S(T) the matrix of second derivatives.
+
+Then the basic iteration is given a value of T, update it to
+
+           T - S.i() * D
+
+where .i() denotes inverse.
+
+If F was quadratic this would give exactly the right answer (except it
+might get a minimum rather than a maximum). Since F is not usually
+quadratic, the simple procedure would be to recalculate S and D with the
+new value of T and keep iterating until the process converges. This is
+known as the method of conjugate gradients.
+
+In practice, this method may not converge. FindMaximum2 considers an
+iteration of the form
+
+           T - x * S.i() * D
+
+where x is a number. It tries x = 1 and uses the values of F and its
+slope with respect to x at x = 0 and x = 1 to fit a cubic in x. It then
+choses x to maximise the resulting function. This gives our new value of
+T. The program checks that the value of F is getting better and carries
+out a variety of strategies if it is not.
+
+The program also has a second strategy. If the successive values of T
+seem to be lying along a curve - eg we are following along a curved
+ridge, the program will try to fit this ridge and project along it. This
+does not work at present and is commented out.
+
+FindMaximum2 has three virtual functions which need to be over-ridden by
+a derived class.
+
+   void Value(const ColumnVector& T, bool wg, Real& f, bool& oorg);
+
+T is the column vector of parameters. The function returns the value of
+the function to f, but may instead set oorg to true if the parameter
+values are not valid. If wg is true it may also calculate and store the
+second derivative information.
+
+   bool NextPoint(ColumnVector& H, Real& d);
+
+Using the value of T provided in the previous call of Value, find the
+conjugate gradients adjustment to T, that is - S.i() * D. Also return
+
+           d = D.t() * S.i() * D.
+
+NextPoint should return true if it considers that the process has
+converged (d very small) and false otherwise. The previous call of Value
+will have set wg to true, so that S will be available.
+
+   Real LastDerivative(const ColumnVector& H);
+
+Return the scalar product of H and the vector of derivatives at the last
+value of T.
+
+The function Fit is the function that calls the iteration.
+
+   void Fit(ColumnVector&, int);
+
+The arguments are the trial parameter values as a ColumnVector and the
+maximum number of iterations. The program calls a DataException if the
+initial parameters are not valid and a ConvergenceException if the
+process fails to converge.
+
+
+   class NonLinearLeastSquares
+
+This class is derived from FindMaximum2 and carries out a non-linear
+least squares fit. It uses a QR decomposition to carry out the
+operations required by FindMaximum2.
+
+A prototype class R1_Col_I_D is provided. The user needs to derive a
+class from this which includes functions the predicted value of each
+observation its derivatives. An object from this class has to be
+provided to class NonLinearLeastSquares.
+
+Suppose we observe n normal random variables with the same unknown
+variance and such the i-th one has expected value given by f(i,P)
+where P is a column vector of unknown parameters and f is a known
+function. We wish to estimate P.
+
+First derive a class from R1_Col_I_D and override Real operator()(int i)
+to give the value of the function f in terms of i and the ColumnVector
+para defined in class R1_CoL_I_D. Also override ReturnMatrix
+Derivatives() to give the derivates of f at para and the value of i
+used in the preceeding call to operator(). Return the result as a
+RowVector. Construct an object from this class. Suppose in what follows
+it is called pred.
+
+Now constuct a NonLinearLeastSquaresObject accessing pred and optionally
+an iteration limit and an accuracy critierion.
+
+   NonLinearLeastSquares NLLS(pred, 1000, 0.0001);
+
+The accuracy critierion should be somewhat less than one and 0.0001 is
+about the smallest sensible value.
+
+Define a ColumnVector P containing a guess at the value of the unknown
+parameter, and a ColumnVector Y containing the unknown data. Call
+
+   NLLS.Fit(Y,P);
+
+If the process converges, P will contain the estimates of the unknown
+parameters. If it does not converge an exception will be generated.
+
+The following member functions can be called after you have done a fit.
+
+Real ResidualVariance() const;
+
+The estimate of the variance of the observations.
+
+void GetResiduals(ColumnVector& Z) const;
+
+The residuals of the individual observations.
+
+void GetStandardErrors(ColumnVector&);
+
+The standard errors of the observations.
+
+void GetCorrelations(SymmetricMatrix&);
+
+The correlations of the observations.
+
+void GetHatDiagonal(DiagonalMatrix&) const;
+
+Forms a diagonal matrix of values between 0 and 1. If the i-th value is
+larger than, say 0.2, then the i-th data value could have an undue
+influence on your estimates.
+
+
+*/
+
+class FindMaximum2
+{
+   virtual void Value(const ColumnVector&, bool, Real&, bool&) = 0;
+   virtual bool NextPoint(ColumnVector&, Real&) = 0;
+   virtual Real LastDerivative(const ColumnVector&) = 0;
+public:
+   void Fit(ColumnVector&, int);
+   virtual ~FindMaximum2() {}            // to keep gnu happy
+};
+
+class R1_Col_I_D
+{
+   // The prototype for a Real function of a ColumnVector and an
+   // integer.
+   // You need to derive your function from this one and put in your
+   // function for operator() and Derivatives() at least.
+   // You may also want to set up a constructor to enter in additional
+   // parameter values (that will not vary during the solve).
+
+protected:
+   ColumnVector para;                 // Current x value
+
+public:
+   virtual bool IsValid() { return true; }
+                                       // is the current x value OK
+   virtual Real operator()(int i) = 0; // i-th function value at current para
+   virtual void Set(const ColumnVector& X) { para = X; }
+                                       // set current para
+   bool IsValid(const ColumnVector& X)
+      { Set(X); return IsValid(); }
+                                       // set para, check OK
+   Real operator()(int i, const ColumnVector& X)
+      { Set(X); return operator()(i); }
+                                       // set para, return value
+   virtual ReturnMatrix Derivatives() = 0;
+                                       // return derivatives as RowVector
+   virtual ~R1_Col_I_D() {}            // to keep gnu happy
+};
+
+
+class NonLinearLeastSquares : public FindMaximum2
+{
+   // these replace the corresponding functions in FindMaximum2
+   void Value(const ColumnVector&, bool, Real&, bool&);
+   bool NextPoint(ColumnVector&, Real&);
+   Real LastDerivative(const ColumnVector&);
+
+   Matrix X;                         // the things we need to do the
+   ColumnVector Y;                   // QR triangularisation
+   UpperTriangularMatrix U;          // see the write-up in newmata.txt
+   ColumnVector M;
+   Real errorvar, criterion;
+   int n_obs, n_param;
+   const ColumnVector* DataPointer;
+   RowVector Derivs;
+   SymmetricMatrix Covariance;
+   DiagonalMatrix SE;
+   R1_Col_I_D& Pred;                 // Reference to predictor object
+   int Lim;                          // maximum number of iterations
+
+public:
+   NonLinearLeastSquares(R1_Col_I_D& pred, int lim=1000, Real crit=0.0001)
+      : criterion(crit), Pred(pred), Lim(lim) {}
+   void Fit(const ColumnVector&, ColumnVector&);
+   Real ResidualVariance() const { return errorvar; }
+   void GetResiduals(ColumnVector& Z) const { Z = Y; }
+   void GetStandardErrors(ColumnVector&);
+   void GetCorrelations(SymmetricMatrix&);
+   void GetHatDiagonal(DiagonalMatrix&) const;
+
+private:
+   void MakeCovariance();
+};
+
+
+// The next class is the prototype class for calculating the
+// log-likelihood.
+// I assume first derivatives are available and something like the 
+// Fisher Information or variance/covariance matrix of the first
+// derivatives or minus the matrix of second derivatives is
+// available. This matrix must be positive definite.
+
+class LL_D_FI
+{
+protected:
+   ColumnVector para;                  // current parameter values
+   bool wg;                         // true if FI matrix wanted
+
+public:
+   virtual void Set(const ColumnVector& X) { para = X; }
+                                       // set parameter values
+   virtual void WG(bool wgx) { wg = wgx; }
+                                       // set wg
+
+   virtual bool IsValid() { return true; }
+                                       // return true is para is OK
+   bool IsValid(const ColumnVector& X, bool wgx=true)
+      { Set(X); WG(wgx); return IsValid(); }
+
+   virtual Real LogLikelihood() = 0;   // return the loglikelihhod
+   Real LogLikelihood(const ColumnVector& X, bool wgx=true)
+      { Set(X); WG(wgx); return LogLikelihood(); }
+
+   virtual ReturnMatrix Derivatives() = 0;
+                                       // column vector of derivatives
+   virtual ReturnMatrix FI() = 0;      // Fisher Information matrix
+   virtual ~LL_D_FI() {}               // to keep gnu happy
+};
+
+// This is the class for doing the maximum likelihood estimation
+
+class MLE_D_FI : public FindMaximum2
+{
+   // these replace the corresponding functions in FindMaximum2
+   void Value(const ColumnVector&, bool, Real&, bool&);
+   bool NextPoint(ColumnVector&, Real&);
+   Real LastDerivative(const ColumnVector&);
+
+   // the things we need for the analysis
+   LL_D_FI& LL;                        // reference to log-likelihood
+   int Lim;                            // maximum number of iterations
+   Real Criterion;                     // convergence criterion
+   ColumnVector Derivs;                // for the derivatives
+   LowerTriangularMatrix LT;           // Cholesky decomposition of FI
+   SymmetricMatrix Covariance;
+   DiagonalMatrix SE;
+
+public:
+   MLE_D_FI(LL_D_FI& ll, int lim=1000, Real criterion=0.0001)
+      : LL(ll), Lim(lim), Criterion(criterion) {}
+   void Fit(ColumnVector& Parameters);
+   void GetStandardErrors(ColumnVector&);
+   void GetCorrelations(SymmetricMatrix&);
+
+private:
+   void MakeCovariance();
+};
+
+
+#ifdef use_namespace
+}
+#endif
+
+
+
+#endif
+
+// body file: newmatnl.cpp
+
+
+
+///@}
+
Index: /trunk/BNC/newmat/newmatrc.h
===================================================================
--- /trunk/BNC/newmat/newmatrc.h	(revision 8900)
+++ /trunk/BNC/newmat/newmatrc.h	(revision 8901)
Index: /trunk/BNC/newmat/newmatrm.cpp
===================================================================
--- /trunk/BNC/newmat/newmatrm.cpp	(revision 8900)
+++ /trunk/BNC/newmat/newmatrm.cpp	(revision 8901)
Index: /trunk/BNC/newmat/newmatrm.h
===================================================================
--- /trunk/BNC/newmat/newmatrm.h	(revision 8900)
+++ /trunk/BNC/newmat/newmatrm.h	(revision 8901)
Index: /trunk/BNC/newmat/nm_misc.cpp
===================================================================
--- /trunk/BNC/newmat/nm_misc.cpp	(revision 8900)
+++ /trunk/BNC/newmat/nm_misc.cpp	(revision 8901)
Index: /trunk/BNC/newmat/precisio.h
===================================================================
--- /trunk/BNC/newmat/precisio.h	(revision 8900)
+++ /trunk/BNC/newmat/precisio.h	(revision 8901)
@@ -25,4 +25,6 @@
 #endif
 
+using namespace std;
+	
 /// Floating point precision.
 class FloatingPointPrecision
@@ -30,20 +32,20 @@
 public:
    static int Dig()              // number of decimal digits or precision
-      { return std::numeric_limits<Real>::digits10 ; }
+      { return numeric_limits<Real>::digits10 ; }
 
    static Real Epsilon()         // smallest number such that 1+Eps!=Eps
-      { return std::numeric_limits<Real>::epsilon(); }
+      { return numeric_limits<Real>::epsilon(); }
 
    static int Mantissa()         // bits in mantisa
-      { return std::numeric_limits<Real>::digits; }
+      { return numeric_limits<Real>::digits; }
 
    static Real Maximum()         // maximum value
-      { return std::numeric_limits<Real>::max(); }
+      { return numeric_limits<Real>::max(); }
 
    static int MaximumDecimalExponent()  // maximum decimal exponent
-      { return std::numeric_limits<Real>::max_exponent10; }
+      { return numeric_limits<Real>::max_exponent10; }
 
    static int MaximumExponent()  // maximum binary exponent
-      { return std::numeric_limits<Real>::max_exponent; }
+      { return numeric_limits<Real>::max_exponent; }
 
    static Real LnMaximum()       // natural log of maximum
@@ -51,11 +53,11 @@
 
    static Real Minimum()         // minimum positive value
-      { return std::numeric_limits<Real>::min(); } 
+      { return numeric_limits<Real>::min(); } 
 
    static int MinimumDecimalExponent() // minimum decimal exponent
-      { return std::numeric_limits<Real>::min_exponent10; }
+      { return numeric_limits<Real>::min_exponent10; }
 
    static int MinimumExponent()  // minimum binary exponent
-      { return std::numeric_limits<Real>::min_exponent; }
+      { return numeric_limits<Real>::min_exponent; }
 
    static Real LnMinimum()       // natural log of minimum
@@ -63,10 +65,10 @@
 
    static int Radix()            // exponent radix
-      { return std::numeric_limits<Real>::radix; }
+      { return numeric_limits<Real>::radix; }
 
    static int Rounds()           // addition rounding (1 = does round)
    {
-	  return std::numeric_limits<Real>::round_style ==
-		 std::round_to_nearest ? 1 : 0;
+	  return numeric_limits<Real>::round_style ==
+		 round_to_nearest ? 1 : 0;
    }
 
Index: /trunk/BNC/newmat/solution.cpp
===================================================================
--- /trunk/BNC/newmat/solution.cpp	(revision 8901)
+++ /trunk/BNC/newmat/solution.cpp	(revision 8901)
@@ -0,0 +1,205 @@
+/// \ingroup newmat
+///@{
+
+/// \file solution.cpp
+/// One dimensional solve routine.
+
+// Copyright (C) 1994: R B Davies
+
+
+#define WANT_STREAM                  // include.h will get stream fns
+#define WANT_MATH                    // include.h will get math fns
+
+#include "include.h"
+#include "myexcept.h"
+
+#include "solution.h"
+
+#ifdef use_namespace
+namespace RBD_COMMON {
+#endif
+
+
+void R1_R1::Set(Real X)
+{
+   if ((!minXinf && X <= minX) || (!maxXinf && X >= maxX))
+       Throw(SolutionException("X value out of range"));
+   x = X; xSet = true;
+}
+
+R1_R1::operator Real()
+{
+   if (!xSet) Throw(SolutionException("Value of X not set"));
+   Real y = operator()();
+   return y;
+}
+
+unsigned long SolutionException::Select;
+
+SolutionException::SolutionException(const char* a_what) : BaseException()
+{
+   Select = BaseException::Select;
+   AddMessage("Error detected by solution package\n");
+   AddMessage(a_what); AddMessage("\n");
+   if (a_what) Tracer::AddTrace();
+}
+
+inline Real square(Real x) { return x*x; }
+
+void OneDimSolve::LookAt(int V)
+{
+   lim--;
+   if (!lim) Throw(SolutionException("Does not converge"));
+   Last = V;
+   Real yy = function(x[V]) - YY;
+   Finish = (fabs(yy) <= accY) || (Captured && fabs(x[L]-x[U]) <= accX );
+   y[V] = vpol*yy;
+}
+
+void OneDimSolve::HFlip() { hpol=-hpol; State(U,C,L); }
+
+void OneDimSolve::VFlip()
+   { vpol = -vpol; y[0] = -y[0]; y[1] = -y[1]; y[2] = -y[2]; }
+
+void OneDimSolve::Flip()
+{
+   hpol=-hpol; vpol=-vpol; State(U,C,L);
+   y[0] = -y[0]; y[1] = -y[1]; y[2] = -y[2];
+}
+
+void OneDimSolve::State(int I, int J, int K) { L=I; C=J; U=K; }
+
+void OneDimSolve::Linear(int I, int J, int K)
+{
+   x[J] = (x[I]*y[K] - x[K]*y[I])/(y[K] - y[I]);
+   // cout << "Linear\n";
+}
+
+void OneDimSolve::Quadratic(int I, int J, int K)
+{
+   // result to overwrite I
+   Real YJK, YIK, YIJ, XKI, XKJ;
+   YJK = y[J] - y[K]; YIK = y[I] - y[K]; YIJ = y[I] - y[J];
+   XKI = (x[K] - x[I]);
+   XKJ = (x[K]*y[J] - x[J]*y[K])/YJK;
+   if ( square(YJK/YIK)>(x[K] - x[J])/XKI ||
+      square(YIJ/YIK)>(x[J] - x[I])/XKI )
+   {
+      x[I] = XKJ;
+      // cout << "Quadratic - exceptional\n";
+   }
+   else
+   {
+      XKI = (x[K]*y[I] - x[I]*y[K])/YIK;
+      x[I] = (XKJ*y[I] - XKI*y[J])/YIJ;
+      // cout << "Quadratic - normal\n";
+   }
+}
+
+Real OneDimSolve::Solve(Real Y, Real X, Real Dev, int Lim)
+{
+   enum Loop { start, captured1, captured2, binary, finish };
+   Tracer et("OneDimSolve::Solve");
+   lim=Lim; Captured = false;
+   if ( Dev == 0.0 ) Throw(SolutionException("Dev is zero"));
+   L=0; C=1; U=2; vpol=1; hpol=1; y[C]=0.0; y[U]=0.0;
+   if (Dev<0.0) { hpol=-1; Dev = -Dev; }
+   YY=Y;                                // target value
+   x[L] = X;                            // initial trial value
+   if (!function.IsValid(X))
+      Throw(SolutionException("Starting value is invalid"));
+   Loop TheLoop = start;
+   for (;;)
+   {
+      switch (TheLoop)
+      {
+      case start:
+         LookAt(L); if (Finish) { TheLoop = finish; break; }
+         if (y[L]>0.0) VFlip();               // so Y[L] < 0
+
+         x[U] = X + Dev * hpol;
+         if (!function.maxXinf && x[U] > function.maxX)
+            x[U] = (function.maxX + X) / 2.0;
+         if (!function.minXinf && x[U] < function.minX)
+            x[U] = (function.minX + X) / 2.0;
+
+         LookAt(U); if (Finish) { TheLoop = finish; break; }
+         if (y[U] > 0.0) { TheLoop = captured1; Captured = true; break; }
+         if (y[U] == y[L])
+            Throw(SolutionException("Function is flat"));
+         if (y[U] < y[L]) HFlip();             // Change direction
+         State(L,U,C);
+         for (i=0; i<20; i++)
+         {
+            // cout << "Searching for crossing point\n";
+            // Have L C then crossing point, Y[L]<Y[C]<0
+            x[U] = x[C] + Dev * hpol;
+            if (!function.maxXinf && x[U] > function.maxX)
+            x[U] = (function.maxX + x[C]) / 2.0;
+            if (!function.minXinf && x[U] < function.minX)
+            x[U] = (function.minX + x[C]) / 2.0;
+
+            LookAt(U); if (Finish) { TheLoop = finish; break; }
+            if (y[U] > 0) { TheLoop = captured2; Captured = true; break; }
+            if (y[U] < y[C])
+                Throw(SolutionException("Function is not monotone"));
+            Dev *= 2.0;
+            State(C,U,L);
+         }
+         if (TheLoop != start ) break;
+         Throw(SolutionException("Cannot locate a crossing point"));
+
+      case captured1:
+         // cout << "Captured - 1\n";
+         // We have 2 points L and U with crossing between them
+         Linear(L,C,U);                   // linear interpolation
+                                          // - result to C
+         LookAt(C); if (Finish) { TheLoop = finish; break; }
+         if (y[C] > 0.0) Flip();            // Want y[C] < 0
+         if (y[C] < 0.5*y[L]) { State(C,L,U); TheLoop = binary; break; }
+
+      case captured2:
+         // cout << "Captured - 2\n";
+         // We have L,C before crossing, U after crossing
+         Quadratic(L,C,U);                // quad interpolation
+                                          // - result to L
+         State(C,L,U);
+         if ((x[C] - x[L])*hpol <= 0.0 || (x[C] - x[U])*hpol >= 0.0)
+            { TheLoop = captured1; break; }
+         LookAt(C); if (Finish) { TheLoop = finish; break; }
+         // cout << "Through first stage\n";
+         if (y[C] > 0.0) Flip();
+         if (y[C] > 0.5*y[L]) { TheLoop = captured2; break; }
+         else { State(C,L,U); TheLoop = captured1; break; }
+
+      case binary:
+         // We have L, U around crossing - do binary search
+         // cout << "Binary\n";
+         for (i=3; i; i--)
+         {
+            x[C] = 0.5*(x[L]+x[U]);
+            LookAt(C); if (Finish) { TheLoop = finish; break; }
+            if (y[C]>0.0) State(L,U,C); else State(C,L,U);
+         }
+         if (TheLoop != binary) break;
+         TheLoop = captured1; break;
+
+      case finish:
+	 return x[Last];
+
+      }
+   }
+}
+
+bool R1_R1::IsValid(Real X)
+{
+   Set(X);
+   return (minXinf || x > minX) && (maxXinf || x < maxX);
+}
+
+#ifdef use_namespace
+}
+#endif
+
+
+///@}
Index: /trunk/BNC/newmat/solution.h
===================================================================
--- /trunk/BNC/newmat/solution.h	(revision 8901)
+++ /trunk/BNC/newmat/solution.h	(revision 8901)
@@ -0,0 +1,99 @@
+/// \ingroup newmat
+///@{
+
+/// \file solution.h
+/// One dimensional solve routine.
+
+#include "myexcept.h"
+
+#ifdef use_namespace
+namespace RBD_COMMON {
+#endif
+
+
+// Solve the equation f(x)=y for x where f is a monotone continuous
+// function of x
+// Essentially Brent s method
+
+// You need to derive a class from R1_R1 and override "operator()"
+// with the function you want to solve.
+// Use an object from this class in OneDimSolve
+
+class R1_R1
+{
+   // the prototype for a Real function of a Real variable
+   // you need to derive your function from this one and put in your
+   // function for operator() at least. You probably also want to set up a
+   // constructor to put in additional parameter values (e.g. that will not
+   // vary during a solve)
+
+protected:
+   Real x;                             // Current x value
+   bool xSet;                          // true if a value assigned to x
+
+public:
+   Real minX, maxX;                    // range of value x
+   bool minXinf, maxXinf;              // true if these are infinite
+   R1_R1() : xSet(false), minXinf(true), maxXinf(true) {}
+   virtual Real operator()() = 0;      // function value at current x
+                                       // set current x
+   virtual void Set(Real X);           // set x, check OK
+   Real operator()(Real X) { Set(X); return operator()(); }
+                                       // set x, return value
+   virtual bool IsValid(Real X);
+   operator Real();                    // implicit conversion
+   virtual ~R1_R1() {}                 // keep gnu happy
+};
+
+class SolutionException : public BaseException
+{
+public:
+   static unsigned long Select;
+   SolutionException(const char* a_what = 0);
+};
+
+class OneDimSolve
+{
+   R1_R1& function;                     // reference to the function
+   Real accX;                           // accuracy in X direction
+   Real accY;                           // accuracy in Y direction
+   int lim;                             // maximum number of iterations
+
+public:
+   OneDimSolve(R1_R1& f, Real AccY = 0.0001, Real AccX = 0.0)
+      : function(f), accX(AccX), accY(AccY) {}
+                       // f is an R1_R1 function
+   Real Solve(Real Y, Real X, Real Dev, int Lim=100);
+                       // Solve for x in Y=f(x)
+                       // X is the initial trial value of x
+                       // X+Dev is the second trial value
+                       // program returns a value of x such that
+                       // |Y-f(x)| <= accY or |f.inv(Y)-x| <= accX
+
+private:
+   Real x[3], y[3];                         // Trial values of X and Y
+   int L,C,U,Last;                          // Locations of trial values
+   int vpol, hpol;                          // polarities
+   Real YY;                                 // target value
+   int i;
+   void LookAt(int);                        // get new value of function
+   bool Finish;                             // true if LookAt finds conv.
+   bool Captured;                           // true when target surrounded
+   void VFlip();
+   void HFlip();
+   void Flip();
+   void State(int I, int J, int K);
+   void Linear(int, int, int);
+   void Quadratic(int, int, int);
+};
+
+
+#ifdef use_namespace
+}
+#endif
+
+// body file: solution.cpp
+
+
+///@}
+
Index: /trunk/BNC/newmat/sort.cpp
===================================================================
--- /trunk/BNC/newmat/sort.cpp	(revision 8900)
+++ /trunk/BNC/newmat/sort.cpp	(revision 8901)
Index: /trunk/BNC/newmat/submat.cpp
===================================================================
--- /trunk/BNC/newmat/submat.cpp	(revision 8900)
+++ /trunk/BNC/newmat/submat.cpp	(revision 8901)
@@ -257,4 +257,7 @@
       if (row_number != gmx->Nrows() || col_number != gmx->Ncols())
          Throw(IncompatibleDimensionsException());
+      if (gm->type().is_symmetric() && 
+         ( ! gmx->type().is_symmetric() || row_skip != col_skip) )
+         Throw(ProgramException("Illegal operation on symmetric"));
       MatrixRow mrx(gmx, LoadOnEntry);
       MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
@@ -277,8 +280,8 @@
 }
 
-void GetSubMatrix::operator-=(const BaseMatrix& bmx)
-{
-   REPORT
-   Tracer tr("SubMatrix(-=)"); GeneralMatrix* gmx = 0;
+void GetSubMatrix::SP_eq(const BaseMatrix& bmx)
+{
+   REPORT
+   Tracer tr("SubMatrix(SP_eq)"); GeneralMatrix* gmx = 0;
    // MatrixConversionCheck mcc;         // Check for loss of info
    Try
@@ -287,4 +290,39 @@
       if (row_number != gmx->Nrows() || col_number != gmx->Ncols())
          Throw(IncompatibleDimensionsException());
+      if (gm->type().is_symmetric() && 
+         ( ! gmx->type().is_symmetric() || row_skip != col_skip) )
+         Throw(ProgramException("Illegal operation on symmetric"));
+      MatrixRow mrx(gmx, LoadOnEntry);
+      MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
+                                     // do need LoadOnEntry
+      MatrixRowCol sub; int i = row_number;
+      while (i--)
+      {
+         mr.SubRowCol(sub, col_skip, col_number);   // put values in sub
+         sub.Multiply(mrx); mr.Next(); mrx.Next();
+      }
+      gmx->tDelete();
+   }
+
+   CatchAll
+   {
+      if (gmx) gmx->tDelete();
+      ReThrow;
+   }
+}
+
+void GetSubMatrix::operator-=(const BaseMatrix& bmx)
+{
+   REPORT
+   Tracer tr("SubMatrix(-=)"); GeneralMatrix* gmx = 0;
+   // MatrixConversionCheck mcc;         // Check for loss of info
+   Try
+   {
+      SetUpLHS(); gmx = ((BaseMatrix&)bmx).Evaluate();
+      if (row_number != gmx->Nrows() || col_number != gmx->Ncols())
+         Throw(IncompatibleDimensionsException());
+      if (gm->type().is_symmetric() && 
+         ( ! gmx->type().is_symmetric() || row_skip != col_skip) )
+         Throw(ProgramException("Illegal operation on symmetric"));
       MatrixRow mrx(gmx, LoadOnEntry);
       MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
Index: /trunk/BNC/newmat/svd.cpp
===================================================================
--- /trunk/BNC/newmat/svd.cpp	(revision 8900)
+++ /trunk/BNC/newmat/svd.cpp	(revision 8901)
Index: /trunk/BNC/src/PPP_SSR_I/pppFilter.cpp
===================================================================
--- /trunk/BNC/src/PPP_SSR_I/pppFilter.cpp	(revision 8900)
+++ /trunk/BNC/src/PPP_SSR_I/pppFilter.cpp	(revision 8901)
@@ -332,5 +332,5 @@
   xRec(3) = z();
 
-  double rho0 = (satData->xx - xRec).norm_Frobenius();
+  double rho0 = (satData->xx - xRec).NormFrobenius();
   double dPhi = t_CST::omega * rho0 / t_CST::c;
 
@@ -341,5 +341,5 @@
   xRec += _tides->displacement(_time, xRec);
 
-  satData->rho = (satData->xx - xRec).norm_Frobenius();
+  satData->rho = (satData->xx - xRec).NormFrobenius();
 
   double tropDelay = delay_saast(satData->eleSat) +
@@ -796,12 +796,12 @@
     // --------------------------------------
     ColumnVector rho = rRec - rSat;
-    rho /= rho.norm_Frobenius();
+    rho /= rho.NormFrobenius();
 
     // GPS Satellite unit Vectors sz, sy, sx
     // -------------------------------------
-    ColumnVector sz = -rSat / rSat.norm_Frobenius();
+    ColumnVector sz = -rSat / rSat.NormFrobenius();
 
     ColumnVector xSun = t_astro::Sun(Mjd);
-    xSun /= xSun.norm_Frobenius();
+    xSun /= xSun.NormFrobenius();
 
     ColumnVector sy = crossproduct(sz, xSun);
@@ -839,5 +839,5 @@
     // ----------------
     double alpha = DotProduct(dipSat,dipRec) /
-                      (dipSat.norm_Frobenius() * dipRec.norm_Frobenius());
+                      (dipSat.NormFrobenius() * dipRec.NormFrobenius());
 
     if (alpha >  1.0) alpha =  1.0;
@@ -861,5 +861,5 @@
   Tracer tracer("t_pppFilter::cmpEle");
   ColumnVector rr = satData->xx - _xcBanc.Rows(1,3);
-  double       rho = rr.norm_Frobenius();
+  double       rho = rr.NormFrobenius();
 
   double neu[3];
Index: /trunk/BNC/src/combination/bnccomb.cpp
===================================================================
--- /trunk/BNC/src/combination/bnccomb.cpp	(revision 8900)
+++ /trunk/BNC/src/combination/bnccomb.cpp	(revision 8901)
@@ -1141,6 +1141,6 @@
       }
       else {
-        double normMax = maxDiff[prn]->_diffRao.norm_Frobenius();
-        double norm    = corr->_diffRao.norm_Frobenius();
+        double normMax = maxDiff[prn]->_diffRao.NormFrobenius();
+        double norm    = corr->_diffRao.NormFrobenius();
         if (norm > normMax) {
           maxDiff[prn] = corr;
@@ -1165,5 +1165,5 @@
       }
       else if (corr == maxDiff[prn]) {
-        double norm = corr->_diffRao.norm_Frobenius();
+        double norm = corr->_diffRao.NormFrobenius();
         if (norm > MAX_DISPLACEMENT) {
           out << _resTime.datestr().c_str()    << " "
Index: /trunk/BNC/src/rinex/reqcanalyze.cpp
===================================================================
--- /trunk/BNC/src/rinex/reqcanalyze.cpp	(revision 8900)
+++ /trunk/BNC/src/rinex/reqcanalyze.cpp	(revision 8901)
@@ -278,5 +278,5 @@
         ++nSatUsed;
         ColumnVector dx = xSat.Rows(1,3) - xyzSta;
-        double rho = dx.norm_Frobenius();
+        double rho = dx.NormFrobenius();
         AA(nSatUsed,1) = dx(1) / rho;
         AA(nSatUsed,2) = dx(2) / rho;