Index: trunk/BNC/newmat/bandmat.cpp
===================================================================
--- trunk/BNC/newmat/bandmat.cpp	(revision 9433)
+++ trunk/BNC/newmat/bandmat.cpp	(revision 9434)
Index: trunk/BNC/newmat/cholesky.cpp
===================================================================
--- trunk/BNC/newmat/cholesky.cpp	(revision 9433)
+++ trunk/BNC/newmat/cholesky.cpp	(revision 9434)
Index: trunk/BNC/newmat/controlw.h
===================================================================
--- trunk/BNC/newmat/controlw.h	(revision 9433)
+++ trunk/BNC/newmat/controlw.h	(revision 9434)
Index: trunk/BNC/newmat/evalue.cpp
===================================================================
--- trunk/BNC/newmat/evalue.cpp	(revision 9433)
+++ trunk/BNC/newmat/evalue.cpp	(revision 9434)
Index: trunk/BNC/newmat/fft.cpp
===================================================================
--- trunk/BNC/newmat/fft.cpp	(revision 9433)
+++ trunk/BNC/newmat/fft.cpp	(revision 9434)
Index: trunk/BNC/newmat/hholder.cpp
===================================================================
--- trunk/BNC/newmat/hholder.cpp	(revision 9433)
+++ trunk/BNC/newmat/hholder.cpp	(revision 9434)
@@ -335,148 +335,4 @@
 }
 
-// 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 9433)
+++ trunk/BNC/newmat/include.h	(revision 9434)
@@ -1,3 +1,3 @@
-/// \defgroup rbd_common RBD common library
+/// \defgroup rbd_common RBD common library 
 ///@{
 
@@ -44,7 +44,5 @@
 //#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
@@ -65,19 +63,7 @@
 #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
 
@@ -104,5 +90,6 @@
       #include <fstream>
    #endif
-   using namespace std;
+   ////using namespace std;
+   #define USE_STD_NAMESPACE
 #else
 
Index: trunk/BNC/newmat/jacobi.cpp
===================================================================
--- trunk/BNC/newmat/jacobi.cpp	(revision 9433)
+++ trunk/BNC/newmat/jacobi.cpp	(revision 9434)
Index: trunk/BNC/newmat/myexcept.cpp
===================================================================
--- trunk/BNC/newmat/myexcept.cpp	(revision 9433)
+++ trunk/BNC/newmat/myexcept.cpp	(revision 9434)
@@ -26,4 +26,7 @@
 #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 9433)
+++ trunk/BNC/newmat/myexcept.h	(revision 9434)
@@ -75,5 +75,4 @@
    static void AddTrace();               // insert trace in exception record
    static Tracer* last;                  // points to Tracer list
-   static void clear() {}                // for compatibility
    friend class BaseException;
 };
@@ -94,5 +93,4 @@
    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 9433)
+++ trunk/BNC/newmat/newfft.cpp	(revision 9434)
Index: trunk/BNC/newmat/newmat.h
===================================================================
--- trunk/BNC/newmat/newmat.h	(revision 9433)
+++ trunk/BNC/newmat/newmat.h	(revision 9434)
@@ -447,4 +447,5 @@
 class GeneralMatrix : public BaseMatrix         // declarable matrix types
 {
+   virtual GeneralMatrix* Image() const;        // copy of matrix
 protected:
    int tag_val;                                 // shows whether can reuse
@@ -464,8 +465,13 @@
    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);
@@ -478,5 +484,4 @@
              // 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);
@@ -485,8 +490,4 @@
 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(); }
@@ -502,5 +503,4 @@
    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,5 +526,4 @@
    void Inject(const GeneralMatrix& GM) { inject(GM); }
    void operator+=(const BaseMatrix&);
-   void SP_eq(const BaseMatrix&);
    void operator-=(const BaseMatrix&);
    void operator*=(const BaseMatrix&);
@@ -580,5 +579,4 @@
 //   ReturnMatrix Reverse() const;                // reverse order of elements
    void cleanup();                              // to clear store
-   virtual GeneralMatrix* Image() const;        // copy of matrix
 
    friend class Matrix;
@@ -627,4 +625,5 @@
 class Matrix : public GeneralMatrix
 {
+   GeneralMatrix* Image() const;                // copy of matrix
 public:
    Matrix() {}
@@ -669,13 +668,8 @@
    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)
@@ -685,4 +679,5 @@
 class SquareMatrix : public Matrix
 {
+   GeneralMatrix* Image() const;                // copy of matrix
 public:
    SquareMatrix() {}
@@ -706,13 +701,8 @@
    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)
 };
@@ -721,4 +711,5 @@
 class nricMatrix : public Matrix
 {
+   GeneralMatrix* Image() const;                // copy of matrix
    Real** row_pointer;                          // points to rows
    void MakeRowPointer();                       // build rowpointer
@@ -752,13 +743,8 @@
    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)
 };
@@ -767,4 +753,5 @@
 class SymmetricMatrix : public GeneralMatrix
 {
+   GeneralMatrix* Image() const;                // copy of matrix
 public:
    SymmetricMatrix() {}
@@ -802,13 +789,8 @@
    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)
 };
@@ -817,4 +799,5 @@
 class UpperTriangularMatrix : public GeneralMatrix
 {
+   GeneralMatrix* Image() const;                // copy of matrix
 public:
    UpperTriangularMatrix() {}
@@ -854,14 +837,9 @@
    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)
 };
@@ -870,4 +848,5 @@
 class LowerTriangularMatrix : public GeneralMatrix
 {
+   GeneralMatrix* Image() const;                // copy of matrix
 public:
    LowerTriangularMatrix() {}
@@ -906,14 +885,9 @@
    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)
 };
@@ -922,4 +896,5 @@
 class DiagonalMatrix : public GeneralMatrix
 {
+   GeneralMatrix* Image() const;                // copy of matrix
 public:
    DiagonalMatrix() {}
@@ -967,14 +942,9 @@
 //   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)
 };
@@ -983,4 +953,5 @@
 class RowVector : public Matrix
 {
+   GeneralMatrix* Image() const;                // copy of matrix
 public:
    RowVector() { nrows_val = 1; }
@@ -1026,14 +997,9 @@
    // 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)
 };
@@ -1042,4 +1008,5 @@
 class ColumnVector : public Matrix
 {
+   GeneralMatrix* Image() const;                // copy of matrix
 public:
    ColumnVector() { ncols_val = 1; }
@@ -1079,14 +1046,9 @@
 //   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)
 };
@@ -1102,4 +1064,5 @@
    void ludcmp();
    void get_aux(CroutMatrix&);                  // for copying indx[] etc
+   GeneralMatrix* Image() const;                // copy of matrix
 public:
    CroutMatrix(const BaseMatrix&);
@@ -1125,5 +1088,4 @@
    bool even_exchanges() const { return d; }
    void swap(CroutMatrix& gm);
-   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(CroutMatrix)
 };
@@ -1134,4 +1096,5 @@
 class BandMatrix : public GeneralMatrix
 {
+   GeneralMatrix* Image() const;                // copy of matrix
 protected:
    void CornerClear() const;                    // set unused elements to zero
@@ -1198,5 +1161,4 @@
    void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
    void swap(BandMatrix& gm);
-   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(BandMatrix)
 };
@@ -1205,4 +1167,5 @@
 class UpperBandMatrix : public BandMatrix
 {
+   GeneralMatrix* Image() const;                // copy of matrix
 public:
    UpperBandMatrix() {}
@@ -1237,5 +1200,4 @@
    void swap(UpperBandMatrix& gm)
       { BandMatrix::swap((BandMatrix&)gm); }
-   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(UpperBandMatrix)
 };
@@ -1244,4 +1206,5 @@
 class LowerBandMatrix : public BandMatrix
 {
+   GeneralMatrix* Image() const;                // copy of matrix
 public:
    LowerBandMatrix() {}
@@ -1276,5 +1239,4 @@
    void swap(LowerBandMatrix& gm)
       { BandMatrix::swap((BandMatrix&)gm); }
-   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(LowerBandMatrix)
 };
@@ -1283,4 +1245,5 @@
 class SymmetricBandMatrix : public GeneralMatrix
 {
+   GeneralMatrix* Image() const;                // copy of matrix
    void CornerClear() const;                    // set unused elements to zero
    short SimpleAddOK(const GeneralMatrix* gm);
@@ -1337,5 +1300,4 @@
    void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
    void swap(SymmetricBandMatrix& gm);
-   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(SymmetricBandMatrix)
 };
@@ -1352,4 +1314,5 @@
    void ludcmp();
    void get_aux(BandLUMatrix&);                 // for copying indx[] etc
+   GeneralMatrix* Image() const;                // copy of matrix
 public:
    BandLUMatrix(const BaseMatrix&);
@@ -1379,5 +1342,4 @@
    MatrixBandWidth bandwidth() const;
    void swap(BandLUMatrix& gm);
-   GeneralMatrix* Image() const;                // copy of matrix
    NEW_DELETE(BandLUMatrix)
 };
@@ -1388,4 +1350,5 @@
 class IdentityMatrix : public GeneralMatrix
 {
+   GeneralMatrix* Image() const;          // copy of matrix
 public:
    IdentityMatrix() {}
@@ -1423,5 +1386,4 @@
    void swap(IdentityMatrix& gm)
       { GeneralMatrix::swap((GeneralMatrix&)gm); }
-   GeneralMatrix* Image() const;          // copy of matrix
    NEW_DELETE(IdentityMatrix)
 };
@@ -1447,5 +1409,4 @@
    void operator=(const BaseMatrix&);
    void operator+=(const BaseMatrix&);
-   void SP_eq(const BaseMatrix&);
    void operator-=(const BaseMatrix&);
    void operator*=(const BaseMatrix&);
@@ -1820,5 +1781,4 @@
    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 9433)
+++ trunk/BNC/newmat/newmat.pro	(revision 9434)
@@ -10,15 +10,15 @@
 MOC_DIR     = .moc
 
-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
+HEADERS += controlw.h include.h myexcept.h  \
+           newmatap.h newmat.h newmatio.h   \
+           newmatrc.h newmatrm.h precisio.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 newmatnl.cpp nm_misc.cpp   \
-	   sort.cpp submat.cpp svd.cpp solution.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 nm_misc.cpp sort.cpp    \
+           submat.cpp svd.cpp
 
Index: trunk/BNC/newmat/newmat1.cpp
===================================================================
--- trunk/BNC/newmat/newmat1.cpp	(revision 9433)
+++ trunk/BNC/newmat/newmat1.cpp	(revision 9434)
Index: trunk/BNC/newmat/newmat2.cpp
===================================================================
--- trunk/BNC/newmat/newmat2.cpp	(revision 9433)
+++ trunk/BNC/newmat/newmat2.cpp	(revision 9434)
Index: trunk/BNC/newmat/newmat3.cpp
===================================================================
--- trunk/BNC/newmat/newmat3.cpp	(revision 9433)
+++ trunk/BNC/newmat/newmat3.cpp	(revision 9434)
Index: trunk/BNC/newmat/newmat4.cpp
===================================================================
--- trunk/BNC/newmat/newmat4.cpp	(revision 9433)
+++ trunk/BNC/newmat/newmat4.cpp	(revision 9434)
Index: trunk/BNC/newmat/newmat5.cpp
===================================================================
--- trunk/BNC/newmat/newmat5.cpp	(revision 9433)
+++ trunk/BNC/newmat/newmat5.cpp	(revision 9434)
Index: trunk/BNC/newmat/newmat6.cpp
===================================================================
--- trunk/BNC/newmat/newmat6.cpp	(revision 9433)
+++ trunk/BNC/newmat/newmat6.cpp	(revision 9434)
@@ -495,23 +495,10 @@
    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());
 }
 
@@ -599,19 +586,4 @@
    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 9433)
+++ trunk/BNC/newmat/newmat7.cpp	(revision 9434)
@@ -214,19 +214,11 @@
 }
 
-static void SP(GeneralMatrix* gm, const GeneralMatrix* gm2)
-{
-   REPORT
-   const Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
+static void SP(GeneralMatrix* gm, GeneralMatrix* gm2)
+{
+   REPORT
+   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 9433)
+++ trunk/BNC/newmat/newmat8.cpp	(revision 9434)
Index: trunk/BNC/newmat/newmat9.cpp
===================================================================
--- trunk/BNC/newmat/newmat9.cpp	(revision 9433)
+++ trunk/BNC/newmat/newmat9.cpp	(revision 9434)
@@ -44,6 +44,5 @@
    MatrixRow mr((GeneralMatrix*)&X, LoadOnEntry);
    int w = s.width();  int nr = X.Nrows();  ios_format_flags f = s.flags();
-   if (f & ios::scientific) s.setf(ios::scientific, ios::floatfield);
-   else s.setf(ios::fixed, ios::floatfield);
+   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 9433)
+++ trunk/BNC/newmat/newmatap.h	(revision 9434)
@@ -26,11 +26,5 @@
 void QRZ(Matrix&, UpperTriangularMatrix&);
 
-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); } 
+void QRZ(const Matrix&, Matrix&, Matrix&);
 
 inline void HHDecompose(Matrix& X, LowerTriangularMatrix& L)
@@ -44,10 +38,4 @@
 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); }
@@ -55,14 +43,4 @@
 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
@@ -78,12 +56,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 9433)
+++ trunk/BNC/newmat/newmatex.cpp	(revision 9434)
Index: trunk/BNC/newmat/newmatio.h
===================================================================
--- trunk/BNC/newmat/newmatio.h	(revision 9433)
+++ trunk/BNC/newmat/newmatio.h	(revision 9434)
@@ -16,9 +16,7 @@
 #include "newmat.h"
 
-#ifdef use_namespace
-namespace NEWMAT {
+#ifdef USE_STD_NAMESPACE
+using namespace std;
 #endif
-
-
 
 // **************************** input/output *****************************/
Index: trunk/BNC/newmat/newmatnl.cpp
===================================================================
--- trunk/BNC/newmat/newmatnl.cpp	(revision 9433)
+++ 	(revision )
@@ -1,266 +1,0 @@
-/// \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 9433)
+++ 	(revision )
@@ -1,327 +1,0 @@
-/// \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 9433)
+++ trunk/BNC/newmat/newmatrc.h	(revision 9434)
Index: trunk/BNC/newmat/newmatrm.cpp
===================================================================
--- trunk/BNC/newmat/newmatrm.cpp	(revision 9433)
+++ trunk/BNC/newmat/newmatrm.cpp	(revision 9434)
Index: trunk/BNC/newmat/newmatrm.h
===================================================================
--- trunk/BNC/newmat/newmatrm.h	(revision 9433)
+++ trunk/BNC/newmat/newmatrm.h	(revision 9434)
Index: trunk/BNC/newmat/nm_misc.cpp
===================================================================
--- trunk/BNC/newmat/nm_misc.cpp	(revision 9433)
+++ trunk/BNC/newmat/nm_misc.cpp	(revision 9434)
Index: trunk/BNC/newmat/precisio.h
===================================================================
--- trunk/BNC/newmat/precisio.h	(revision 9433)
+++ trunk/BNC/newmat/precisio.h	(revision 9434)
@@ -25,6 +25,4 @@
 #endif
 
-using namespace std;
-	
 /// Floating point precision.
 class FloatingPointPrecision
@@ -32,20 +30,20 @@
 public:
    static int Dig()              // number of decimal digits or precision
-      { return numeric_limits<Real>::digits10 ; }
+      { return std::numeric_limits<Real>::digits10 ; }
 
    static Real Epsilon()         // smallest number such that 1+Eps!=Eps
-      { return numeric_limits<Real>::epsilon(); }
+      { return std::numeric_limits<Real>::epsilon(); }
 
    static int Mantissa()         // bits in mantisa
-      { return numeric_limits<Real>::digits; }
+      { return std::numeric_limits<Real>::digits; }
 
    static Real Maximum()         // maximum value
-      { return numeric_limits<Real>::max(); }
+      { return std::numeric_limits<Real>::max(); }
 
    static int MaximumDecimalExponent()  // maximum decimal exponent
-      { return numeric_limits<Real>::max_exponent10; }
+      { return std::numeric_limits<Real>::max_exponent10; }
 
    static int MaximumExponent()  // maximum binary exponent
-      { return numeric_limits<Real>::max_exponent; }
+      { return std::numeric_limits<Real>::max_exponent; }
 
    static Real LnMaximum()       // natural log of maximum
@@ -53,11 +51,11 @@
 
    static Real Minimum()         // minimum positive value
-      { return numeric_limits<Real>::min(); } 
+      { return std::numeric_limits<Real>::min(); } 
 
    static int MinimumDecimalExponent() // minimum decimal exponent
-      { return numeric_limits<Real>::min_exponent10; }
+      { return std::numeric_limits<Real>::min_exponent10; }
 
    static int MinimumExponent()  // minimum binary exponent
-      { return numeric_limits<Real>::min_exponent; }
+      { return std::numeric_limits<Real>::min_exponent; }
 
    static Real LnMinimum()       // natural log of minimum
@@ -65,10 +63,10 @@
 
    static int Radix()            // exponent radix
-      { return numeric_limits<Real>::radix; }
+      { return std::numeric_limits<Real>::radix; }
 
    static int Rounds()           // addition rounding (1 = does round)
    {
-	  return numeric_limits<Real>::round_style ==
-		 round_to_nearest ? 1 : 0;
+	  return std::numeric_limits<Real>::round_style ==
+		 std::round_to_nearest ? 1 : 0;
    }
 
Index: trunk/BNC/newmat/solution.cpp
===================================================================
--- trunk/BNC/newmat/solution.cpp	(revision 9433)
+++ 	(revision )
@@ -1,205 +1,0 @@
-/// \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 9433)
+++ 	(revision )
@@ -1,99 +1,0 @@
-/// \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 9433)
+++ trunk/BNC/newmat/sort.cpp	(revision 9434)
Index: trunk/BNC/newmat/submat.cpp
===================================================================
--- trunk/BNC/newmat/submat.cpp	(revision 9433)
+++ trunk/BNC/newmat/submat.cpp	(revision 9434)
@@ -257,7 +257,4 @@
       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);
@@ -280,8 +277,8 @@
 }
 
-void GetSubMatrix::SP_eq(const BaseMatrix& bmx)
-{
-   REPORT
-   Tracer tr("SubMatrix(SP_eq)"); GeneralMatrix* gmx = 0;
+void GetSubMatrix::operator-=(const BaseMatrix& bmx)
+{
+   REPORT
+   Tracer tr("SubMatrix(-=)"); GeneralMatrix* gmx = 0;
    // MatrixConversionCheck mcc;         // Check for loss of info
    Try
@@ -290,39 +287,4 @@
       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 9433)
+++ trunk/BNC/newmat/svd.cpp	(revision 9434)