Index: trunk/npstat/nm/HeatEq1DNeumannBoundary.cc
===================================================================
--- trunk/npstat/nm/HeatEq1DNeumannBoundary.cc	(revision 694)
+++ trunk/npstat/nm/HeatEq1DNeumannBoundary.cc	(revision 695)
@@ -1,238 +1,241 @@
 #include <cassert>
 #include <stdexcept>
 #include <cfloat>
 #include <complex>
 #include <cmath>
 
 #include "npstat/nm/HeatEq1DNeumannBoundary.hh"
 #include "npstat/nm/ArrayND.hh"
 
 #include "npstat/stat/WeightedStatAccumulator.hh"
 
 namespace npstat {
     HeatEq1DNeumannBoundary::HeatEq1DNeumannBoundary(
         const double* thermalDiffusivity, const unsigned nx,
         const double h, const double i_dt)
         : gammas_(nx),
           D_(nx),
           DL_(nx-1U),
           DU_(nx-1U),
           DU2_(nx-2U),
           m0_(nx, nx, 1),
           m1_(nx, nx),
           IPIV_(nx),
           dt_(i_dt),
           intervalWidth_(h),
           cnt_(0),
           nx_(nx)
     {
         if (nx < 3U) throw std::invalid_argument(
             "In npstat::HeatEq1DNeumannBoundary constructor: "
             "not enough spatial discretization cells");
         if (h <= 0.0) throw std::invalid_argument(
             "In npstat::HeatEq1DNeumannBoundary constructor: "
             "spatial step must be positive");
         if (i_dt <= 0.0) throw std::invalid_argument(
             "In npstat::HeatEq1DNeumannBoundary constructor: "
             "time step must be positive");
         assert(thermalDiffusivity);
 
         const double r = i_dt/h/h;
         for (unsigned i=0; i<nx; ++i)
         {
             if (thermalDiffusivity[i] <= 0.0) throw std::invalid_argument(
                 "In npstat::HeatEq1DNeumannBoundary constructor: "
                 "thermal diffusivity must be positive");
             gammas_[i] = thermalDiffusivity[i]*r;
             D_[i] = 1.0 + gammas_[i];
         }
 
         const unsigned nxm1 = nx - 1U;
         const unsigned nxm2 = nx - 2U;
         D_[0] = 1.0 + gammas_[0]/2.0 + (gammas_[1]-gammas_[0])/8.0;
         D_[nxm1] = 1.0 + gammas_[nxm1]/2.0 - (gammas_[nxm1]-gammas_[nxm2])/8.0;
 
         for (unsigned j=1; j<nx; ++j)
         {
             const unsigned jm1 = j - 1U;
             const unsigned jp1 = j < nxm1 ? j + 1U : nxm1;
             DL_[jm1] = -gammas_[j]/2.0 + (gammas_[jp1] - gammas_[jm1])/8.0;
         }
 
         for (unsigned j=0; j<nxm1; ++j)
         {
             const unsigned jm1 = j ? j - 1U : 0U;
             const unsigned jp1 = j + 1U;
             DU_[j] = -gammas_[j]/2.0 - (gammas_[jp1] - gammas_[jm1])/8.0;
         }
 
         int INFO = 1;
         int N = nx;
         dgttrf_(&N, &DL_[0], &D_[0], &DU_[0], &DU2_[0], &IPIV_[0], &INFO);
         Private::check_lapack_status(INFO, "npstat::HeatEq1DNeumannBoundary constructor", "DGTTRF");
     }
 
     void HeatEq1DNeumannBoundary::step()
     {
         Matrix<double> *gf, *rhs;
         if (cnt_ % 2)
         {
             gf = &m1_;
             rhs = &m0_;            
         }
         else
         {
             gf = &m0_;
             rhs = &m1_;
         }
 
         const unsigned nxm1 = nx_ - 1U;
         for (unsigned cell=0; cell<nx_; ++cell)
         {
             const double* u = (*gf)[cell];
             double* r = (*rhs)[cell];
 
             for (unsigned j=0; j<nx_; ++j)
             {
                 const unsigned jm1 = j ? j - 1U : 0U;
                 const unsigned jp1 = j < nxm1 ? j + 1U : nxm1;
                 r[j] = u[j] + gammas_[j]/2.0*((u[jp1] - u[j]) + (u[jm1] - u[j])) +
                     (gammas_[jp1] - gammas_[jm1])/8.0*(u[jp1] - u[jm1]);
             }
         }
 
         char T[] = "N";
         int N = nx_;
         int INFO = 1;
         dgttrs_(T, &N, &N, &DL_[0], &D_[0], &DU_[0], &DU2_[0], &IPIV_[0],
                 const_cast<double*>(rhs->data()), &N, &INFO, 1);
         Private::check_lapack_status(INFO, "npstat::HeatEq1DNeumannBoundary::step", "DGTTRS");
 
         ++cnt_;
     }
 
     Matrix<double> HeatEq1DNeumannBoundary::GF() const
     {
         const Matrix<double>* gf;
         if (cnt_ % 2)
             gf = &m1_;
         else
             gf = &m0_;
         return gf->T();
     }
 
     double HeatEq1DNeumannBoundary::gfWidth(const unsigned cell) const
     {
         if (cell >= nx_) throw std::invalid_argument(
             "In npstat::HeatEq1DNeumannBoundary::gfWidth: "
             "cell index out of range");
         const Matrix<double>* gf;
         if (cnt_ % 2)
             gf = &m1_;
         else
             gf = &m0_;
         const double* u = (*gf)[cell];
 
         WeightedStatAccumulator acc;
         for (unsigned j=0; j<nx_; ++j)
             if (u[j] > 0.0)
                 acc.accumulate(j*intervalWidth_, u[j]);
 
         const double effco = acc.count();
-        return acc.stdev()*sqrt((effco - 1.0)/effco);
+        if (effco > 1.0)
+            return acc.stdev()*sqrt((effco - 1.0)/effco);
+        else
+            return 0.0;
     }
 
     Matrix<double> HeatEq1DNeumannBoundary::doublyStochasticGF(
         const double tol, const unsigned maxIter) const
     {
         const Matrix<double>* gf;
         if (cnt_ % 2)
             gf = &m1_;
         else
             gf = &m0_;
 
         if (maxIter)
         {
             Matrix<double> gfc(gf->T());
             gfc *= nx_;
             const unsigned long len = gfc.length();
             ArrayND<double> arr(nx_, nx_);
             arr.setData(gfc.data(), len);
             arr.makeCopulaSteps(tol, maxIter);
             gfc.setData(arr.data(), len);
             gfc /= nx_;
             return gfc;
         }
         else
         {
             Matrix<double> Jn(nx_, nx_);
             Jn.constFill(1.0/nx_);
             Matrix<double> W(nx_, nx_, 1);
             W -= Jn;
             return W*gf->T()*W + Jn;
         }
     }
 
     Matrix<double> HeatEq1DNeumannBoundary::diffScheme() const
     {
         Matrix<double> lmat(nx_, nx_, 0);
         const unsigned nxm1 = nx_ - 1U;
         const unsigned nxm2 = nx_ - 2U;
 
         lmat[0][0] = 1.0 + gammas_[0]/2.0 + (gammas_[1]-gammas_[0])/8.0;
         for (unsigned i=1; i<nxm1; ++i)
             lmat[i][i] = 1.0 + gammas_[i];
         lmat[nxm1][nxm1] = 1.0 + gammas_[nxm1]/2.0 - (gammas_[nxm1]-gammas_[nxm2])/8.0;
 
         for (unsigned j=1; j<nx_; ++j)
         {
             const unsigned jm1 = j - 1U;
             const unsigned jp1 = j < nxm1 ? j + 1U : nxm1;
             lmat[j][jm1] = -gammas_[j]/2.0 + (gammas_[jp1] - gammas_[jm1])/8.0;
         }
 
         for (unsigned j=0; j<nxm1; ++j)
         {
             const unsigned jm1 = j ? j - 1U : 0U;
             const unsigned jp1 = j + 1U;
             lmat[j][jp1] = -gammas_[j]/2.0 - (gammas_[jp1] - gammas_[jm1])/8.0;
         }
 
         Matrix<double> rmat(nx_, nx_, 0);
         rmat[0][0] = 1.0 - gammas_[0]/2.0 - (gammas_[1]-gammas_[0])/8.0;
         for (unsigned i=1; i<nxm1; ++i)
             rmat[i][i] = 1.0 - gammas_[i];
         rmat[nxm1][nxm1] = 1.0 - gammas_[nxm1]/2.0 + (gammas_[nxm1]-gammas_[nxm2])/8.0;
 
         for (unsigned j=1; j<nx_; ++j)
         {
             const unsigned jm1 = j - 1U;
             const unsigned jp1 = j < nxm1 ? j + 1U : nxm1;
             rmat[j][jm1] = gammas_[j]/2.0 - (gammas_[jp1] - gammas_[jm1])/8.0;
         }
 
         for (unsigned j=0; j<nxm1; ++j)
         {
             const unsigned jm1 = j ? j - 1U : 0U;
             const unsigned jp1 = j + 1U;
             rmat[j][jp1] = gammas_[j]/2.0 + (gammas_[jp1] - gammas_[jm1])/8.0;
         }
 
         return lmat.triDiagInv() * rmat;
     }
 
     double HeatEq1DNeumannBoundary::maxAbsEigen() const
     {
         const Matrix<double>& L = diffScheme();
         std::vector<std::complex<double> > eigen(nx_);
         L.genEigen(&eigen[0], nx_);
         double maxEigen = -DBL_MAX;
         for (unsigned i=0; i<nx_; ++i)
         {
             const double absval = std::abs(eigen[i]);
             if (absval > maxEigen)
                 maxEigen = absval;
         }
         return maxEigen;
     }
 }