Index: trunk/BNC/src/orbComp/sp3Comp.cpp
===================================================================
--- trunk/BNC/src/orbComp/sp3Comp.cpp	(revision 11076)
+++ trunk/BNC/src/orbComp/sp3Comp.cpp	(revision 11083)
@@ -178,42 +178,62 @@
   }
 
-  int nPar = epochs.size() + clkSats.size();
-  SymmetricMatrix NN(nPar); NN = 0.0;
-  ColumnVector    bb(nPar); bb = 0.0;
-
-  // Create Matrix A'A and vector b
-  // ------------------------------
-  for (unsigned ie = 0; ie < epochs.size(); ie++) {
+  // Least-squares adjustment of an epoch-specific offset (common to all
+  // satellites) and a satellite-specific offset (common to all epochs), the
+  // sum of the satellite offsets constrained to zero (weight 1). The normal
+  // equations of the epoch offsets are diagonal, so these offsets are
+  // eliminated (Schur complement) and only the satellite offsets are solved
+  // for; the epoch offsets follow from them. The solution is identical to
+  // that of the full system, but the effort grows only linearly with the
+  // number of epochs (before: cubic).
+  // -----------------------------------------------------------------------
+  int nEpo = epochs.size();
+  int nSat = clkSats.size();
+  int nPar = nEpo + nSat;
+
+  vector<vector<int> >    satIdx(nEpo);   // satellites (index) of each epoch
+  ColumnVector    nEp(nEpo); nEp = 0.0;   // number of satellites per epoch
+  ColumnVector    bEp(nEpo); bEp = 0.0;   // sum of the differences per epoch
+  SymmetricMatrix NS(nSat);  NS  = 1.0;   // regularization: sum of satellite offsets
+  ColumnVector    bS(nSat);  bS  = 0.0;
+  for (int ie = 0; ie < nEpo; ie++) {
     const map<t_prn, double>& dc = epochs[ie]->_dc;
-    Matrix       AA(dc.size(), nPar); AA = 0.0;
-    ColumnVector ll(dc.size());       ll = 0.0;
-    map<t_prn, double>::const_iterator it;
-    int ii = -1;
-    for (it = dc.begin(); it != dc.end(); it++) {
-      const t_prn& prn = it->first;
-      if (satIndex(clkSats, prn) != -1) {
-        ++ii;
-        int index = epochs.size() + satIndex(clkSats, prn);
-        AA[ii][ie]    = 1.0; // epoch-specfic offset (common for all satellites)
-        AA[ii][index] = 1.0; // satellite-specific offset (common for all epochs)
-        ll[ii]        = it->second;
-      }
-    }
-    SymmetricMatrix dN; dN << AA.t() * AA;
-    NN += dN;
-    bb += AA.t() * ll;
-  }
-
-  // Regularize NN
-  // -------------
-  RowVector HH(nPar);
-  HH.columns(1, epochs.size())      = 0.0;
-  HH.columns(epochs.size()+1, nPar) = 1.0;
-  SymmetricMatrix dN; dN << HH.t() * HH;
-  NN += dN;
-
-  // Estimate Parameters
-  // -------------------
-  ColumnVector xx = NN.i() * bb;
+    for (map<t_prn, double>::const_iterator it = dc.begin(); it != dc.end(); it++) {
+      int is = satIndex(clkSats, it->first);
+      if (is != -1) {
+        satIdx[ie].push_back(is);
+        nEp[ie]     += 1.0;
+        bEp[ie]     += it->second;
+        NS(is+1, is+1) += 1.0;
+        bS[is]      += it->second;
+      }
+    }
+  }
+
+  // Eliminate the epoch offsets
+  // ---------------------------
+  for (int ie = 0; ie < nEpo; ie++) {
+    const vector<int>& idx = satIdx[ie];
+    for (unsigned ia = 0; ia < idx.size(); ia++) {
+      bS[idx[ia]] -= bEp[ie] / nEp[ie];
+      for (unsigned ib = 0; ib <= ia; ib++) {
+        NS(idx[ia]+1, idx[ib]+1) -= 1.0 / nEp[ie];
+      }
+    }
+  }
+
+  // Estimate Parameters (epoch offsets first, then satellite offsets)
+  // ---------------------------------------------------------------
+  ColumnVector yy = NS.i() * bS;
+  ColumnVector xx(nPar);
+  for (int ie = 0; ie < nEpo; ie++) {
+    double sum = bEp[ie];
+    for (unsigned ia = 0; ia < satIdx[ie].size(); ia++) {
+      sum -= yy[satIdx[ie][ia]];
+    }
+    xx[ie] = sum / nEp[ie];
+  }
+  for (int is = 0; is < nSat; is++) {
+    xx[nEpo + is] = yy[is];
+  }
 
   // Compute clock residuals