Add counting of subnormal/zero numbers in sparse factors

examples/sparse/testMMdouble.cpp

@@ -60,19 +60,26 @@ test(int argc, char* argv[], CSRMatrix<scalar_t,integer_t>& A) {
   }
   spss.solve(b.data(), x.data());
 
-  integer_t neg, zero, pos;
-  auto err = spss.inertia(neg, zero, pos);
-  std::cout << "# INERTIA neg,zero,pos = "
-            << neg << ", " << zero << ", " << pos
-            <<  " (" << err << ")" << std::endl;
+  // std::size_t subs = 0, zeros = 0;
+  // auto err0 = spss.subnormals(subs, zeros);
+  // std::cout << "# SUBNORMALS = " << subs
+  //           << "   ZEROS = " << zeros
+  //           << " (" << err0 << ")" << std::endl;
+
+  // integer_t neg, zero, pos;
+  // auto err = spss.inertia(neg, zero, pos);
+  // std::cout << "# INERTIA neg,zero,pos = "
+  //           << neg << ", " << zero << ", " << pos
+  //           <<  " (" << err << ")" << std::endl;
 
   std::cout << "# COMPONENTWISE SCALED RESIDUAL = "
             << A.max_scaled_residual(x.data(), b.data()) << std::endl;
 
   strumpack::blas::axpy(N, scalar_t(-1.), x_exact.data(), 1, x.data(), 1);
   auto nrm_error = strumpack::blas::nrm2(N, x.data(), 1);
   auto nrm_x_exact = strumpack::blas::nrm2(N, x_exact.data(), 1);
-  std::cout << "# RELATIVE ERROR = " << (nrm_error/nrm_x_exact) << std::endl;
+  std::cout << "# RELATIVE ERROR (|x-xtrue|_2/|xtrue|_2) = "
+            << (nrm_error/nrm_x_exact) << std::endl;
 }
 
 int main(int argc, char* argv[]) {

examples/sparse/testMMdoubleMPIDist.cpp

@@ -75,12 +75,19 @@ test(int argc, char* argv[], CSRMatrix<scalar,integer>& A) {
   }
   spss.solve(b.data(), x.data());
 
-  integer neg, zero, pos;
-  auto err = spss.inertia(neg, zero, pos);
-  if (!rank)
-    std::cout << "# INERTIA neg,zero,pos = "
-              << neg << ", " << zero << ", " << pos
-              <<  " (" << err << ")" << std::endl;
+  // std::size_t subs = 0, zeros = 0;
+  // auto err0 = spss.subnormals(subs, zeros);
+  // if (!rank)
+  //   std::cout << "# SUBNORMALS = " << subs
+  //             << "   ZEROS = " << zeros
+  //             << " (" << err0 << ")" << std::endl;
+
+  // integer neg, zero, pos;
+  // auto err = spss.inertia(neg, zero, pos);
+  // if (!rank)
+  //   std::cout << "# INERTIA neg,zero,pos = "
+  //             << neg << ", " << zero << ", " << pos
+  //             <<  " (" << err << ")" << std::endl;
 
   auto scaled_res = Adist.max_scaled_residual(x.data(), b.data());
   if (!rank)
@@ -123,13 +130,16 @@ int main(int argc, char* argv[]) {
 
   bool is_complex = false;
   if (rank == 0) {
-    if (A.read_matrix_market(f)) {
-      is_complex = true;
-      if (Acomplex.read_matrix_market(f)) {
-        std::cerr << "Could not read matrix from file." << std::endl;
-        return 1;
+    if (A.read_binary(f)) {
+      if (A.read_matrix_market(f)) {
+        is_complex = true;
+        if (Acomplex.read_matrix_market(f)) {
+          std::cerr << "Could not read matrix from file." << std::endl;
+          return 1;
+        }
       }
     }
+
   }
   MPI_Bcast(&is_complex, sizeof(bool), MPI_BYTE, 0, MPI_COMM_WORLD);
   if (!is_complex)

src/BLR/BLRMatrix.cpp

@@ -37,6 +37,24 @@
 namespace strumpack {
   namespace BLR {
 
+    template<typename scalar_t> std::size_t
+    BLRMatrix<scalar_t>::subnormals() const {
+      std::size_t sn = 0;
+      for (std::size_t j=0; j<colblocks(); j++)
+        for (std::size_t i=0; i<rowblocks(); i++)
+          sn += tile(i, j).subnormals();
+      return sn;
+    }
+
+    template<typename scalar_t> std::size_t
+    BLRMatrix<scalar_t>::zeros() const {
+      std::size_t nz = 0;
+      for (std::size_t j=0; j<colblocks(); j++)
+        for (std::size_t i=0; i<rowblocks(); i++)
+          nz += tile(i, j).zeros();
+      return nz;
+    }
+
     template<typename scalar_t> BLRMatrix<scalar_t>::BLRMatrix
     (DenseM_t& A, const std::vector<std::size_t>& rowtiles,
      const std::vector<std::size_t>& coltiles, const Opts_t& opts)

src/BLR/BLRMatrix.hpp

@@ -97,6 +97,12 @@ namespace strumpack {
       std::size_t nonzeros() const override;
       std::size_t rank() const override;
 
+      /**
+       * Return number of subnormal elements in the matrix.
+       */
+      std::size_t subnormals() const;
+      std::size_t zeros() const;
+
       DenseM_t dense() const;
       void dense(DenseM_t& A) const;
 

src/BLR/BLRTile.hpp

@@ -60,6 +60,10 @@ namespace strumpack {
       virtual std::size_t memory() const = 0;
       virtual std::size_t nonzeros() const = 0;
       virtual std::size_t maximum_rank() const = 0;
+
+      virtual std::size_t subnormals() const = 0;
+      virtual std::size_t zeros() const = 0;
+
       virtual bool is_low_rank() const = 0;
       virtual void dense(DenseM_t& A) const = 0;
       virtual DenseM_t dense() const = 0;

src/BLR/DenseTile.hpp

@@ -59,6 +59,9 @@ namespace strumpack {
       std::size_t maximum_rank() const override { return 0; }
       bool is_low_rank() const override { return false; };
 
+      std::size_t subnormals() const override { return D_.subnormals(); }
+      std::size_t zeros() const override { return D_.zeros(); }
+
       void dense(DenseM_t& A) const override { A = D_; }
       DenseM_t dense() const override { return D_; }
 

src/BLR/LRTile.hpp

@@ -90,6 +90,9 @@ namespace strumpack {
       std::size_t nonzeros() const override { return U_.nonzeros() + V_.nonzeros(); }
       std::size_t maximum_rank() const override { return U_.cols(); }
 
+      std::size_t subnormals() const override { return U_.subnormals() + V_.subnormals(); }
+      std::size_t zeros() const override { return U_.zeros() + V_.zeros(); }
+
       void dense(DenseM_t& A) const override;
       DenseM_t dense() const override;
 

src/SparseSolverBase.cpp

@@ -152,6 +152,17 @@ namespace strumpack {
     return tree()->inertia(neg, zero, pos);
   }
 
+  template<typename scalar_t,typename integer_t> ReturnCode
+  SparseSolverBase<scalar_t,integer_t>::subnormals
+  (std::size_t& ns, std::size_t& nz) {
+    ns = nz = 0;
+    if (!this->factored_) {
+      ReturnCode ierr = this->factor();
+      if (ierr != ReturnCode::SUCCESS) return ierr;
+    }
+    return tree()->subnormals(ns, nz);
+  }
+
   template<typename scalar_t,typename integer_t> void
   SparseSolverBase<scalar_t,integer_t>::draw
   (const std::string& name) const {

src/SparseSolverBase.hpp

@@ -367,6 +367,8 @@ namespace strumpack {
      */
     ReturnCode inertia(integer_t& neg, integer_t& zero, integer_t& pos);
 
+    ReturnCode subnormals(std::size_t& ns, std::size_t& nz);
+
     /**
      * Create a gnuplot script to draw/plot the sparse factors. Only
      * do this for small matrices! It is very slow!

src/dense/DenseMatrix.cpp

@@ -829,8 +829,9 @@ namespace strumpack {
     DenseMatrix tmp(*this);
     auto minmn = std::min(rows(), cols());
     std::vector<scalar_t> S(minmn);
-    int info = blas::gesvd('N', 'N', rows(), cols(), tmp.data(), tmp.ld(),
-                           S.data(), NULL, 1, NULL, 1);
+    int info = blas::gesvd
+      ('N', 'N', rows(), cols(), tmp.data(), tmp.ld(),
+       S.data(), NULL, 1, NULL, 1);
     if (info)
       std::cout << "ERROR in gesvd: info = " << info << std::endl;
     return S;
@@ -844,6 +845,28 @@ namespace strumpack {
       (char(job), char(ul), rows(), data(), ld(), lambda.data());
   }
 
+  template<typename scalar_t> std::size_t
+  DenseMatrix<scalar_t>::subnormals() const {
+    if (!data_) return 0;
+    std::size_t ns = 0;
+    for (std::size_t c=0; c<cols(); c++)
+      for (std::size_t r=0; r<rows(); r++)
+        if (!std::isnormal(std::real(operator()(r, c))) &&
+            !std::isnormal(std::imag(operator()(r, c))))
+          ns++;
+    return ns;
+  }
+
+  template<typename scalar_t> std::size_t
+  DenseMatrix<scalar_t>::zeros() const {
+    if (!data_) return 0;
+    std::size_t nz = 0;
+    for (std::size_t c=0; c<cols(); c++)
+      for (std::size_t r=0; r<rows(); r++)
+        if (operator()(r, c) == scalar_t(0.)) nz++;
+    return nz;
+  }
+
   template<typename scalar_t> void
   DenseMatrix<scalar_t>::write(const std::string& fname) const {
     std::ofstream f(fname, std::ios::out | std::ios::trunc);

src/dense/DenseMatrix.hpp

@@ -839,8 +839,8 @@ namespace strumpack {
      * to be allocated.
      * \param depth current OpenMP task recursion depth
      */
-    void LQ
-    (DenseMatrix<scalar_t>& L, DenseMatrix<scalar_t>& Q, int depth) const;
+    void LQ(DenseMatrix<scalar_t>& L,
+            DenseMatrix<scalar_t>& Q, int depth) const;
 
     /**
      * Builds an orthonormal basis for the columns in this matrix,
@@ -976,6 +976,13 @@ namespace strumpack {
      */
     static DenseMatrix<scalar_t> read(const std::string& fname);
 
+
+    /**
+     * Return number of subnormal elements in the matrix.
+     */
+    std::size_t subnormals() const;
+    std::size_t zeros() const;
+
   private:
     void ID_column_GEQP3
     (DenseMatrix<scalar_t>& X, std::vector<int>& piv,

src/dense/DistributedMatrix.cpp

@@ -1081,6 +1081,30 @@ namespace strumpack {
     trsm(Side::L, UpLo::U, Trans::N, Diag::N, scalar_t(1.), R1, X);
   }
 
+
+  template<typename scalar_t> std::size_t
+  DistributedMatrix<scalar_t>::subnormals() const {
+    if (!data_) return 0;
+    std::size_t ns = 0;
+    for (int c=0; c<lcols_; c++)
+      for (int r=0; r<lrows_; r++)
+        if (!std::isnormal(std::real(operator()(r, c))) &&
+            !std::isnormal(std::imag(operator()(r, c))))
+          ns++;
+    return ns;
+  }
+  template<typename scalar_t> std::size_t
+  DistributedMatrix<scalar_t>::zeros() const {
+    if (!data_) return 0;
+    std::size_t nz = 0;
+    for (int c=0; c<lcols_; c++)
+      for (int r=0; r<lrows_; r++)
+        if (operator()(r, c) == scalar_t(0.))
+          nz++;
+    return nz;
+  }
+
+
   template<typename scalar_t> void gemm
   (Trans ta, Trans tb, scalar_t alpha, const DistributedMatrix<scalar_t>& A,
    const DistributedMatrix<scalar_t>& B,

src/dense/DistributedMatrix.hpp

@@ -675,6 +675,10 @@ namespace strumpack {
                 std::vector<std::size_t>& ind, real_t rel_tol, real_t abs_tol,
                 int max_rank, const BLACSGrid* grid_T);
 
+
+    std::size_t subnormals() const;
+    std::size_t zeros() const;
+
     /**
      * Default row blocksize used for 2D block cyclic
      * dustribution. This is set during CMake configuration.

src/sparse/EliminationTree.cpp

@@ -202,6 +202,12 @@ namespace strumpack {
     return root_->inertia(neg, zero, pos);
   }
 
+  template<typename scalar_t,typename integer_t> ReturnCode
+  EliminationTree<scalar_t,integer_t>::subnormals
+  (std::size_t& ns, std::size_t& nz) const {
+    return root_->subnormals(ns, nz);
+  }
+
   template<typename scalar_t,typename integer_t> void
   EliminationTree<scalar_t,integer_t>::draw
   (const SpMat_t& A, const std::string& name) const {

src/sparse/EliminationTree.hpp

@@ -77,6 +77,8 @@ namespace strumpack {
     virtual ReturnCode inertia(integer_t& neg,
                                integer_t& zero,
                                integer_t& pos) const;
+    virtual ReturnCode subnormals(std::size_t& ns,
+                                  std::size_t& nz) const;
 
     void print_rank_statistics(std::ostream &out) const;
 

src/sparse/EliminationTreeMPI.cpp

@@ -88,6 +88,16 @@ namespace strumpack {
     return info;
   }
 
+  template<typename scalar_t,typename integer_t> ReturnCode
+  EliminationTreeMPI<scalar_t,integer_t>::subnormals
+  (std::size_t& ns, std::size_t& nz) const {
+    auto info = EliminationTree<scalar_t,integer_t>::subnormals(ns, nz);
+    ns = comm_.all_reduce(ns, MPI_SUM);
+    nz = comm_.all_reduce(nz, MPI_SUM);
+    return info;
+  }
+
+
   // explicit template specializations
   template class EliminationTreeMPI<float,int>;
   template class EliminationTreeMPI<double,int>;

src/sparse/EliminationTreeMPI.hpp

@@ -64,6 +64,8 @@ namespace strumpack {
     ReturnCode inertia(integer_t& neg,
                        integer_t& zero,
                        integer_t& pos) const override;
+    virtual ReturnCode subnormals(std::size_t& ns,
+                                  std::size_t& nz) const;
 
   protected:
     const MPIComm& comm_;

src/sparse/fronts/FrontalMatrix.cpp

@@ -397,6 +397,17 @@ namespace strumpack {
     return node_inertia(neg, zero, pos);
   }
 
+  template<typename scalar_t,typename integer_t> ReturnCode
+  FrontalMatrix<scalar_t,integer_t>::subnormals(std::size_t& ns,
+                                                std::size_t& nz) const {
+    ReturnCode el = ReturnCode::SUCCESS, er = ReturnCode::SUCCESS;
+    if (lchild_) el = lchild_->subnormals(ns, nz);
+    if (rchild_) er = rchild_->subnormals(ns, nz);
+    if (el != ReturnCode::SUCCESS) return el;
+    if (er != ReturnCode::SUCCESS) return er;
+    return node_subnormals(ns, nz);
+  }
+
 #if defined(STRUMPACK_USE_MPI)
   template<typename scalar_t,typename integer_t> void
   FrontalMatrix<scalar_t,integer_t>::multifrontal_solve

src/sparse/fronts/FrontalMatrix.hpp

@@ -132,6 +132,8 @@ namespace strumpack {
     ReturnCode inertia(integer_t& neg,
                        integer_t& zero,
                        integer_t& pos) const;
+    ReturnCode subnormals(std::size_t& ns, std::size_t& nz) const;
+
 
     virtual void
     extend_add_to_dense(DenseM_t& paF11, DenseM_t& paF12,
@@ -373,6 +375,11 @@ namespace strumpack {
       return ReturnCode::INACCURATE_INERTIA;
     }
 
+    virtual ReturnCode node_subnormals(std::size_t& ns,
+                                       std::size_t& nz) const {
+      return ReturnCode::INACCURATE_INERTIA;
+    }
+
   private:
     FrontalMatrix(const FrontalMatrix&) = delete;
     FrontalMatrix& operator=(FrontalMatrix const&) = delete;

src/sparse/fronts/FrontalMatrixBLR.cpp

@@ -553,6 +553,21 @@ namespace strumpack {
     return F11blr_.nonzeros() + F12blr_.nonzeros() + F21blr_.nonzeros();
   }
 
+  template<typename scalar_t,typename integer_t> ReturnCode
+  FrontalMatrixBLR<scalar_t,integer_t>::node_subnormals
+  (std::size_t& ns, std::size_t& nz) const {
+    auto dns = F11blr_.subnormals() + F12blr_.subnormals() + F21blr_.subnormals();
+    auto dnz = F11blr_.zeros() + F12blr_.zeros() + F21blr_.zeros();
+    if (dns || dnz)
+      std::cout << "BLR front ds= " << this->dim_sep()
+                << " du= " << this->dim_upd()
+                << " subnormals= " << dns
+                << " zeros= " << dnz << std::endl;
+    ns += dns;
+    nz += dnz;
+    return ReturnCode::SUCCESS;
+  }
+
   template<typename scalar_t,typename integer_t> void
   FrontalMatrixBLR<scalar_t,integer_t>::partition
   (const Opts_t& opts, const SpMat_t& A,