.

include/sctl.hpp

@@ -1,3 +1,5 @@
+// Scientific Computing Template Library
+
 #ifndef _SCTL_HPP_
 #define _SCTL_HPP_
 
@@ -10,17 +12,10 @@
 // Have MPI
 //#define SCTL_HAVE_MPI
 
-// Disable assert checks.
-//#ifndef NDEBUG
-//#define NDEBUG
-//#endif
-
 // Parameters for memory manager
 #define SCTL_MEM_ALIGN 64
-#define SCTL_GLOBAL_MEM_BUFF 1024LL * 1LL  // in MB
-#ifndef NDEBUG
-#define SCTL_MEMDEBUG // Enable memory checks.
-#endif
+#define SCTL_GLOBAL_MEM_BUFF 1024LL * 0LL  // in MB
+//#define SCTL_MEMDEBUG // Enable memory checks.
 
 // Profiling parameters
 #define SCTL_PROFILE 5 // Granularity level

include/sctl/comm.hpp

@@ -24,6 +24,10 @@ class Comm {
 
   Comm();
 
+#ifdef SCTL_HAVE_MPI
+  Comm(const MPI_Comm mpi_comm) { Init(mpi_comm); }
+#endif
+
   Comm(const Comm& c);
 
   static Comm Self();
@@ -62,7 +66,7 @@ class Comm {
 
   template <class Type> void Scan(ConstIterator<Type> sbuf, Iterator<Type> rbuf, int count, CommOp op) const;
 
-  template <class Type> void PartitionW(Vector<Type>& nodeList, const Vector<Long>* wts_ = NULL) const;
+  template <class Type> void PartitionW(Vector<Type>& nodeList, const Vector<Long>* wts_ = nullptr) const;
 
   template <class Type> void PartitionN(Vector<Type>& v, Long N) const;
 
@@ -70,7 +74,7 @@ class Comm {
 
   template <class Type> void HyperQuickSort(const Vector<Type>& arr_, Vector<Type>& SortedElem) const;
 
-  template <class Type> void SortScatterIndex(const Vector<Type>& key, Vector<Long>& scatter_index, const Type* split_key_ = NULL) const;
+  template <class Type> void SortScatterIndex(const Vector<Type>& key, Vector<Long>& scatter_index, const Type* split_key_ = nullptr) const;
 
   template <class Type> void ScatterForward(Vector<Type>& data_, const Vector<Long>& scatter_index) const;
 
@@ -84,8 +88,6 @@ class Comm {
   };
 
 #ifdef SCTL_HAVE_MPI
-  Comm(const MPI_Comm mpi_comm) { Init(mpi_comm); }
-
   void Init(const MPI_Comm mpi_comm) {
     MPI_Comm_dup(mpi_comm, &mpi_comm_);
     MPI_Comm_rank(mpi_comm_, &mpi_rank_);

include/sctl/comm.txx

@@ -90,7 +90,7 @@ inline void Comm::Barrier() const {
 template <class SType> void* Comm::Isend(ConstIterator<SType> sbuf, Long scount, Integer dest, Integer tag) const {
   static_assert(std::is_trivially_copyable<SType>::value, "Data is not trivially copyable!");
 #ifdef SCTL_HAVE_MPI
-  if (!scount) return NULL;
+  if (!scount) return nullptr;
   Vector<MPI_Request>& request = *NewReq();
   request.ReInit(1);
 
@@ -108,14 +108,14 @@ template <class SType> void* Comm::Isend(ConstIterator<SType> sbuf, Long scount,
     send_req.insert(std::pair<Integer, ConstIterator<char>>(tag, (ConstIterator<char>)sbuf));
   else
     memcopy(it->second, (ConstIterator<char>)sbuf, scount * sizeof(SType));
-  return NULL;
+  return nullptr;
 #endif
 }
 
 template <class RType> void* Comm::Irecv(Iterator<RType> rbuf, Long rcount, Integer source, Integer tag) const {
   static_assert(std::is_trivially_copyable<RType>::value, "Data is not trivially copyable!");
 #ifdef SCTL_HAVE_MPI
-  if (!rcount) return NULL;
+  if (!rcount) return nullptr;
   Vector<MPI_Request>& request = *NewReq();
   request.ReInit(1);
 
@@ -129,13 +129,13 @@ template <class RType> void* Comm::Irecv(Iterator<RType> rbuf, Long rcount, Inte
     recv_req.insert(std::pair<Integer, Iterator<char>>(tag, (Iterator<char>)rbuf));
   else
     memcopy((Iterator<char>)rbuf, it->second, rcount * sizeof(RType));
-  return NULL;
+  return nullptr;
 #endif
 }
 
 inline void Comm::Wait(void* req_ptr) const {
 #ifdef SCTL_HAVE_MPI
-  if (req_ptr == NULL) return;
+  if (req_ptr == nullptr) return;
   Vector<MPI_Request>& request = *(Vector<MPI_Request>*)req_ptr;
   // std::vector<MPI_Status> status(request.Dim());
   if (request.Dim()) MPI_Waitall(request.Dim(), &request[0], MPI_STATUSES_IGNORE);  //&status[0]);
@@ -155,7 +155,7 @@ template <class SType, class RType> void Comm::Allgather(ConstIterator<SType> sb
     rbuf[0];
     rbuf[rcount * Size() - 1];
   }
-  MPI_Allgather((scount ? &sbuf[0] : NULL), scount, CommDatatype<SType>::value(), (rcount ? &rbuf[0] : NULL), rcount, CommDatatype<RType>::value(), mpi_comm_);
+  MPI_Allgather((scount ? &sbuf[0] : nullptr), scount, CommDatatype<SType>::value(), (rcount ? &rbuf[0] : nullptr), rcount, CommDatatype<RType>::value(), mpi_comm_);
 #else
   memcopy((Iterator<char>)rbuf, (ConstIterator<char>)sbuf, scount * sizeof(SType));
 #endif
@@ -181,7 +181,7 @@ template <class SType, class RType> void Comm::Allgatherv(ConstIterator<SType> s
     rbuf[0];
     rbuf[rcount_sum - 1];
   }
-  MPI_Allgatherv((scount ? &sbuf[0] : NULL), scount, CommDatatype<SType>::value(), (rcount_sum ? &rbuf[0] : NULL), &rcounts_.begin()[0], &rdispls_.begin()[0], CommDatatype<RType>::value(), mpi_comm_);
+  MPI_Allgatherv((scount ? &sbuf[0] : nullptr), scount, CommDatatype<SType>::value(), (rcount_sum ? &rbuf[0] : nullptr), &rcounts_.begin()[0], &rdispls_.begin()[0], CommDatatype<RType>::value(), mpi_comm_);
 #else
   memcopy((Iterator<char>)(rbuf + rdispls[0]), (ConstIterator<char>)sbuf, scount * sizeof(SType));
 #endif
@@ -199,7 +199,7 @@ template <class SType, class RType> void Comm::Alltoall(ConstIterator<SType> sbu
     rbuf[0];
     rbuf[rcount * Size() - 1];
   }
-  MPI_Alltoall((scount ? &sbuf[0] : NULL), scount, CommDatatype<SType>::value(), (rcount ? &rbuf[0] : NULL), rcount, CommDatatype<RType>::value(), mpi_comm_);
+  MPI_Alltoall((scount ? &sbuf[0] : nullptr), scount, CommDatatype<SType>::value(), (rcount ? &rbuf[0] : nullptr), rcount, CommDatatype<RType>::value(), mpi_comm_);
 #else
   memcopy((Iterator<char>)rbuf, (ConstIterator<char>)sbuf, scount * sizeof(SType));
 #endif
@@ -214,7 +214,7 @@ template <class SType, class RType> void* Comm::Ialltoallv_sparse(ConstIterator<
     if (rcounts[i]) request_count++;
     if (scounts[i]) request_count++;
   }
-  if (!request_count) return NULL;
+  if (!request_count) return nullptr;
   Vector<MPI_Request>& request = *NewReq();
   request.ReInit(request_count);
   Integer request_iter = 0;
@@ -238,7 +238,7 @@ template <class SType, class RType> void* Comm::Ialltoallv_sparse(ConstIterator<
   return &request;
 #else
   memcopy((Iterator<char>)(rbuf + rdispls[0]), (ConstIterator<char>)(sbuf + sdispls[0]), scounts[0] * sizeof(SType));
-  return NULL;
+  return nullptr;
 #endif
 }
 
@@ -277,7 +277,7 @@ template <class Type> void Comm::Alltoallv(ConstIterator<Type> sbuf, ConstIterat
       rtotal += rcounts[i];
     }
 
-    MPI_Alltoallv((stotal ? &sbuf[0] : NULL), &scnt[0], &sdsp[0], CommDatatype<Type>::value(), (rtotal ? &rbuf[0] : NULL), &rcnt[0], &rdsp[0], CommDatatype<Type>::value(), mpi_comm_);
+    MPI_Alltoallv((stotal ? &sbuf[0] : nullptr), &scnt[0], &sdsp[0], CommDatatype<Type>::value(), (rtotal ? &rbuf[0] : nullptr), &rcnt[0], &rdsp[0], CommDatatype<Type>::value(), mpi_comm_);
     return;
     //#endif
   }
@@ -354,7 +354,7 @@ template <class Type> void Comm::PartitionW(Vector<Type>& nodeList, const Vector
 
   Vector<Long> wts;
   Long localWt = 0;
-  if (wts_ == NULL) {  // Construct arrays of wts.
+  if (wts_ == nullptr) {  // Construct arrays of wts.
     wts.ReInit(nlSize);
 #pragma omp parallel for schedule(static)
     for (Long i = 0; i < nlSize; i++) {
@@ -557,7 +557,7 @@ template <class Type> void Comm::SortScatterIndex(const Vector<Type>& key, Vecto
   Vector<Pair_t> psorted;
   HyperQuickSort(parray, psorted);
 
-  if (npes > 1 && split_key_ != NULL) {  // Partition data
+  if (npes > 1 && split_key_ != nullptr) {  // Partition data
     Vector<Type> split_key(npes);
     Allgather(Ptr2ConstItr<Type>(split_key_, 1), 1, split_key.begin(), 1);
 
@@ -1027,7 +1027,7 @@ template <class Type> void Comm::HyperQuickSort(const Vector<Type>& arr_, Vector
             glb_splt_cnts_[i] = glb_splt_cnts[i];
             glb_splt_disp_[i] = glb_splt_disp[i];
           }
-          MPI_Allgatherv((splt_count ? &splitters[0] : NULL), splt_count, CommDatatype<Type>::value(), &glb_splitters[0], &glb_splt_cnts_[0], &glb_splt_disp_[0], CommDatatype<Type>::value(), comm);
+          MPI_Allgatherv((splt_count ? &splitters[0] : nullptr), splt_count, CommDatatype<Type>::value(), &glb_splitters[0], &glb_splt_cnts_[0], &glb_splt_disp_[0], CommDatatype<Type>::value(), comm);
         }
       }
 
@@ -1100,8 +1100,8 @@ template <class Type> void Comm::HyperQuickSort(const Vector<Type>& arr_, Vector
         rbuff.ReInit(rsize);
         rbuff_ext.ReInit(ext_rsize);
         MPI_Status status;
-        MPI_Sendrecv((ssize ? &sbuff[0] : NULL), ssize, CommDatatype<Type>::value(), partner, 0, (rsize ? &rbuff[0] : NULL), rsize, CommDatatype<Type>::value(), partner, 0, comm, &status);
-        if (extra_partner) MPI_Sendrecv(NULL, 0, CommDatatype<Type>::value(), split_id, 0, (ext_rsize ? &rbuff_ext[0] : NULL), ext_rsize, CommDatatype<Type>::value(), split_id, 0, comm, &status);
+        MPI_Sendrecv((ssize ? &sbuff[0] : nullptr), ssize, CommDatatype<Type>::value(), partner, 0, (rsize ? &rbuff[0] : nullptr), rsize, CommDatatype<Type>::value(), partner, 0, comm, &status);
+        if (extra_partner) MPI_Sendrecv(nullptr, 0, CommDatatype<Type>::value(), split_id, 0, (ext_rsize ? &rbuff_ext[0] : nullptr), ext_rsize, CommDatatype<Type>::value(), split_id, 0, comm, &status);
       }
 
       Long nbuff_size = lsize + rsize + ext_rsize;

include/sctl/matrix.hpp

@@ -13,11 +13,18 @@ template <class ValueType> class Vector;
 template <class ValueType> class Permutation;
 
 template <class ValueType> class Matrix {
-
  public:
+  typedef ValueType value_type;
+  typedef ValueType& reference;
+  typedef const ValueType& const_reference;
+  typedef Iterator<ValueType> iterator;
+  typedef ConstIterator<ValueType> const_iterator;
+  typedef Long difference_type;
+  typedef Long size_type;
+
   Matrix();
 
-  Matrix(Long dim1, Long dim2, Iterator<ValueType> data_ = NULL, bool own_data_ = true);
+  Matrix(Long dim1, Long dim2, Iterator<ValueType> data_ = nullptr, bool own_data_ = true);
 
   Matrix(const Matrix<ValueType>& M);
 
@@ -25,7 +32,7 @@ template <class ValueType> class Matrix {
 
   void Swap(Matrix<ValueType>& M);
 
-  void ReInit(Long dim1, Long dim2, Iterator<ValueType> data_ = NULL, bool own_data_ = true);
+  void ReInit(Long dim1, Long dim2, Iterator<ValueType> data_ = nullptr, bool own_data_ = true);
 
   void Write(const char* fname) const;
 

include/sctl/matrix.txx

@@ -30,7 +30,7 @@ template <class ValueType> Matrix<ValueType>::Matrix() {
   dim[0] = 0;
   dim[1] = 0;
   own_data = true;
-  data_ptr = NULL;
+  data_ptr = nullptr;
 }
 
 template <class ValueType> Matrix<ValueType>::Matrix(Long dim1, Long dim2, Iterator<ValueType> data_, bool own_data_) {
@@ -40,11 +40,11 @@ template <class ValueType> Matrix<ValueType>::Matrix(Long dim1, Long dim2, Itera
   if (own_data) {
     if (dim[0] * dim[1] > 0) {
       data_ptr = aligned_new<ValueType>(dim[0] * dim[1]);
-      if (data_ != NULL) {
+      if (data_ != nullptr) {
         memcopy(data_ptr, data_, dim[0] * dim[1]);
       }
     } else
-      data_ptr = NULL;
+      data_ptr = nullptr;
   } else
     data_ptr = data_;
 }
@@ -57,16 +57,16 @@ template <class ValueType> Matrix<ValueType>::Matrix(const Matrix<ValueType>& M)
     data_ptr = aligned_new<ValueType>(dim[0] * dim[1]);
     memcopy(data_ptr, M.data_ptr, dim[0] * dim[1]);
   } else
-    data_ptr = NULL;
+    data_ptr = nullptr;
 }
 
 template <class ValueType> Matrix<ValueType>::~Matrix() {
   if (own_data) {
-    if (data_ptr != NULL) {
+    if (data_ptr != nullptr) {
       aligned_delete(data_ptr);
     }
   }
-  data_ptr = NULL;
+  data_ptr = nullptr;
   dim[0] = 0;
   dim[1] = 0;
 }
@@ -93,7 +93,7 @@ template <class ValueType> void Matrix<ValueType>::ReInit(Long dim1, Long dim2,
   if (own_data_ && own_data && dim[0] * dim[1] >= dim1 * dim2) {
     dim[0] = dim1;
     dim[1] = dim2;
-    if (data_ != NULL) {
+    if (data_ != nullptr) {
       memcopy(data_ptr, data_, dim[0] * dim[1]);
     }
   } else {
@@ -104,7 +104,7 @@ template <class ValueType> void Matrix<ValueType>::ReInit(Long dim1, Long dim2,
 
 template <class ValueType> void Matrix<ValueType>::Write(const char* fname) const {
   FILE* f1 = fopen(fname, "wb+");
-  if (f1 == NULL) {
+  if (f1 == nullptr) {
     std::cout << "Unable to open file for writing:" << fname << '\n';
     return;
   }
@@ -118,7 +118,7 @@ template <class ValueType> void Matrix<ValueType>::Write(const char* fname) cons
 
 template <class ValueType> void Matrix<ValueType>::Read(const char* fname) {
   FILE* f1 = fopen(fname, "r");
-  if (f1 == NULL) {
+  if (f1 == nullptr) {
     std::cout << "Unable to open file for reading:" << fname << '\n';
     return;
   }
@@ -181,7 +181,7 @@ template <class ValueType> Matrix<ValueType> Matrix<ValueType>::operator+(const
   SCTL_ASSERT(M2.Dim(0) == M1.Dim(0) && M2.Dim(1) == M1.Dim(1));
   Profile::Add_FLOP(dim[0] * dim[1]);
 
-  Matrix<ValueType> M_r(M1.Dim(0), M1.Dim(1), NULL);
+  Matrix<ValueType> M_r(M1.Dim(0), M1.Dim(1));
   for (Long i = 0; i < M1.Dim(0) * M1.Dim(1); i++) M_r[0][i] = M1[0][i] + M2[0][i];
   return M_r;
 }
@@ -191,7 +191,7 @@ template <class ValueType> Matrix<ValueType> Matrix<ValueType>::operator-(const
   SCTL_ASSERT(M2.Dim(0) == M1.Dim(0) && M2.Dim(1) == M1.Dim(1));
   Profile::Add_FLOP(dim[0] * dim[1]);
 
-  Matrix<ValueType> M_r(M1.Dim(0), M1.Dim(1), NULL);
+  Matrix<ValueType> M_r(M1.Dim(0), M1.Dim(1));
   for (Long i = 0; i < M1.Dim(0) * M1.Dim(1); i++) M_r[0][i] = M1[0][i] - M2[0][i];
   return M_r;
 }
@@ -200,7 +200,7 @@ template <class ValueType> Matrix<ValueType> Matrix<ValueType>::operator*(const
   SCTL_ASSERT(dim[1] == M.dim[0]);
   Profile::Add_FLOP(2 * (((Long)dim[0]) * dim[1]) * M.dim[1]);
 
-  Matrix<ValueType> M_r(dim[0], M.dim[1], NULL);
+  Matrix<ValueType> M_r(dim[0], M.dim[1]);
   if (M.Dim(0) * M.Dim(1) == 0 || this->Dim(0) * this->Dim(1) == 0) return M_r;
   mat::gemm<ValueType>('N', 'N', M.dim[1], dim[0], dim[1], 1.0, M.data_ptr, M.dim[1], data_ptr, dim[1], 0.0, M_r.data_ptr, M_r.dim[1]);
   return M_r;
@@ -319,28 +319,28 @@ template <class ValueType> Matrix<ValueType>& Matrix<ValueType>::operator/=(Valu
 
 template <class ValueType> Matrix<ValueType> Matrix<ValueType>::operator+(ValueType s) const {
   Long N = dim[0] * dim[1];
-  Matrix<ValueType> M_r(dim(0), dim(1));
+  Matrix<ValueType> M_r(dim[0], dim[1]);
   for (Long i = 0; i < N; i++) M_r.data_ptr[i] = data_ptr[i] + s;
   return M_r;
 }
 
 template <class ValueType> Matrix<ValueType> Matrix<ValueType>::operator-(ValueType s) const {
   Long N = dim[0] * dim[1];
-  Matrix<ValueType> M_r(dim(0), dim(1));
+  Matrix<ValueType> M_r(dim[0], dim[1]);
   for (Long i = 0; i < N; i++) M_r.data_ptr[i] = data_ptr[i] - s;
   return M_r;
 }
 
 template <class ValueType> Matrix<ValueType> Matrix<ValueType>::operator*(ValueType s) const {
   Long N = dim[0] * dim[1];
-  Matrix<ValueType> M_r(dim(0), dim(1));
+  Matrix<ValueType> M_r(dim[0], dim[1]);
   for (Long i = 0; i < N; i++) M_r.data_ptr[i] = data_ptr[i] * s;
   return M_r;
 }
 
 template <class ValueType> Matrix<ValueType> Matrix<ValueType>::operator/(ValueType s) const {
   Long N = dim[0] * dim[1];
-  Matrix<ValueType> M_r(dim(0), dim(1));
+  Matrix<ValueType> M_r(dim[0], dim[1]);
   for (Long i = 0; i < N; i++) M_r.data_ptr[i] = data_ptr[i] / s;
   return M_r;
 }
@@ -582,7 +582,7 @@ template <class ValueType> Permutation<ValueType> Permutation<ValueType>::RandPe
 
 template <class ValueType> Matrix<ValueType> Permutation<ValueType>::GetMatrix() const {
   Long size = perm.Dim();
-  Matrix<ValueType> M_r(size, size, NULL);
+  Matrix<ValueType> M_r(size, size);
   for (Long i = 0; i < size; i++)
     for (Long j = 0; j < size; j++) M_r[i][j] = (perm[j] == i ? scal[j] : 0.0);
   return M_r;
@@ -645,7 +645,7 @@ template <class ValueType> Matrix<ValueType> Permutation<ValueType>::operator*(c
   Long d0 = M.Dim(0);
   Long d1 = M.Dim(1);
 
-  Matrix<ValueType> M_r(d0, d1, NULL);
+  Matrix<ValueType> M_r(d0, d1);
   for (Long i = 0; i < d0; i++) {
     const ValueType s = scal[i];
     ConstIterator<ValueType> M_ = M[i];
@@ -661,7 +661,7 @@ template <class ValueType> Matrix<ValueType> operator*(const Matrix<ValueType>&
   Long d0 = M.Dim(0);
   Long d1 = M.Dim(1);
 
-  Matrix<ValueType> M_r(d0, d1, NULL);
+  Matrix<ValueType> M_r(d0, d1);
   for (Long i = 0; i < d0; i++) {
     ConstIterator<Long> perm_ = P.perm.begin();
     ConstIterator<ValueType> scal_ = P.scal.begin();

include/sctl/mem_mgr.hpp

@@ -28,11 +28,11 @@ template <class ValueType> class ConstIterator {
   void IteratorAssertChecks(Long j = 0) const;
 
  public:
-  typedef std::random_access_iterator_tag iterator_category;
-  typedef const ValueType& reference;
   typedef Long difference_type;
   typedef ValueType value_type;
   typedef const ValueType* pointer;
+  typedef const ValueType& reference;
+  typedef std::random_access_iterator_tag iterator_category;
 
  protected:
   char* base;
@@ -42,12 +42,12 @@ template <class ValueType> class ConstIterator {
   static const Long ValueSize = sizeof(ValueType);
 
  public:
-  ConstIterator(void* base_ = NULL) {
+  ConstIterator(void* base_ = nullptr) {
     base = (char*)base_;
     len = 0;
     offset = 0;
     alloc_ctr = 0;
-    mem_head = NULL;
+    mem_head = nullptr;
   }
 
   // template <size_t LENGTH> ConstIterator(ValueType (&base_)[LENGTH]) {  // DEPRECATED
@@ -161,14 +161,14 @@ template <class ValueType> class ConstIterator {
 template <class ValueType> class Iterator : public ConstIterator<ValueType> {
 
  public:
-  typedef std::random_access_iterator_tag iterator_category;
-  typedef ValueType& reference;
   typedef Long difference_type;
   typedef ValueType value_type;
   typedef ValueType* pointer;
+  typedef ValueType& reference;
+  typedef std::random_access_iterator_tag iterator_category;
 
  public:
-  Iterator(void* base_ = NULL) : ConstIterator<ValueType>(base_) {}
+  Iterator(void* base_ = nullptr) : ConstIterator<ValueType>(base_) {}
 
   template <size_t LENGTH> Iterator(ValueType (&base_)[LENGTH]) : ConstIterator<ValueType>(base_) {}
 
@@ -274,8 +274,6 @@ template <class T> class TypeTraits {
 
  public:
   static uintptr_t ID();
-
-  static bool IsPOD();
 };
 
 /**

include/sctl/mem_mgr.txx

@@ -20,7 +20,7 @@ template <class ValueType> inline ConstIterator<ValueType>::ConstIterator(const
     alloc_ctr = mh.alloc_ctr;
     mem_head = &mh;
   } else
-    mem_head = NULL;
+    mem_head = nullptr;
 }
 
 template <class ValueType> inline void ConstIterator<ValueType>::IteratorAssertChecks(Long j) const {
@@ -30,7 +30,7 @@ template <class ValueType> inline void ConstIterator<ValueType>::IteratorAssertC
   const auto& mem_head = this->mem_head;
   const auto& alloc_ctr = this->alloc_ctr;
 
-  if (*this == NULL) SCTL_WARN("dereferencing a NULL pointer is undefined.");
+  if (*this == nullptr) SCTL_WARN("dereferencing a nullptr is undefined.");
   SCTL_ASSERT_MSG(offset >= 0 && offset + (Long)sizeof(ValueType) <= len, "access to pointer [B" << (offset < 0 ? "" : "+") << offset << ",B" << (offset + (Long)sizeof(ValueType) < 0 ? "" : "+") << offset + (Long)sizeof(ValueType) << ") is outside of the range [B,B+" << len << ").");
   if (mem_head) {
     MemoryManager::MemHead& mh = *(MemoryManager::MemHead*)(mem_head);
@@ -94,23 +94,6 @@ template <class ValueType, Long DIM> inline StaticArray<ValueType, DIM>& StaticA
 
 template <class T> inline uintptr_t TypeTraits<T>::ID() { return (uintptr_t) & ID; }
 
-template <class T> inline bool TypeTraits<T>::IsPOD() { return false; }
-
-#define SCTLDefinePOD(type)                          \
-  template <> bool inline TypeTraits<type>::IsPOD() { \
-    return true;                                      \
-  };
-SCTLDefinePOD(char);
-SCTLDefinePOD(float);
-SCTLDefinePOD(double);
-SCTLDefinePOD(int);
-SCTLDefinePOD(long long);
-SCTLDefinePOD(unsigned long);
-SCTLDefinePOD(char*);
-SCTLDefinePOD(float*);
-SCTLDefinePOD(double*);
-#undef SCTLDefinePOD
-
 inline MemoryManager::MemoryManager(Long N) {
   buff_size = N;
   {  // Allocate buff
@@ -167,7 +150,7 @@ inline MemoryManager::~MemoryManager() {
 }
 
 inline MemoryManager::MemHead& MemoryManager::GetMemHead(char* I) {
-  SCTL_ASSERT_MSG(I != NULL, "NULL pointer exception.");
+  SCTL_ASSERT_MSG(I != nullptr, "nullptr exception.");
   static uintptr_t alignment = SCTL_MEM_ALIGN - 1;
   static uintptr_t header_size = (uintptr_t)(sizeof(MemHead) + alignment) & ~(uintptr_t)alignment;
   return *(MemHead*)(((char*)I) - header_size);
@@ -187,13 +170,13 @@ inline void MemoryManager::CheckMemHead(const MemHead& mem_head) {  // Verify he
 }
 
 inline Iterator<char> MemoryManager::malloc(const Long n_elem, const Long type_size, const uintptr_t type_id) const {
-  if (!n_elem) return NULL;
+  if (!n_elem) return nullptr;
   static uintptr_t alignment = SCTL_MEM_ALIGN - 1;
   static uintptr_t header_size = (uintptr_t)(sizeof(MemHead) + alignment) & ~(uintptr_t)alignment;
 
   Long size = n_elem * type_size + header_size;
   size = (uintptr_t)(size + alignment) & ~(uintptr_t)alignment;
-  char* base = NULL;
+  char* base = nullptr;
 
   omp_set_lock(&omp_lock);
   static Long alloc_ctr = 0;
@@ -292,7 +275,7 @@ inline Iterator<char> MemoryManager::malloc(const Long n_elem, const Long type_s
 }
 
 inline void MemoryManager::free(Iterator<char> p) const {
-  if (p == NULL) return;
+  if (p == nullptr) return;
   static uintptr_t alignment = SCTL_MEM_ALIGN - 1;
   static uintptr_t header_size = (uintptr_t)(sizeof(MemHead) + alignment) & ~(uintptr_t)alignment;
 
@@ -411,7 +394,7 @@ inline void MemoryManager::test() {
 
     for (Integer j = 0; j < 3; j++) {
       tmp = (Iterator<double>)memgr.malloc(M * sizeof(double));
-      SCTL_ASSERT(tmp != NULL);
+      SCTL_ASSERT(tmp != nullptr);
       tt = omp_get_wtime();
 #pragma omp parallel for
       for (Long i = 0; i < M; i += 64) tmp[i] = (double)i;
@@ -428,7 +411,7 @@ inline void MemoryManager::test() {
     std::cout << "Without memory manager: ";
     for (Integer j = 0; j < 3; j++) {
       tmp = (double*)::malloc(M * sizeof(double));
-      SCTL_ASSERT(tmp != NULL);
+      SCTL_ASSERT(tmp != nullptr);
       tt = omp_get_wtime();
 #pragma omp parallel for
       for (Long i = 0; i < M; i += 64) tmp[i] = (double)i;
@@ -479,19 +462,19 @@ inline void MemoryManager::delete_node(Long indx) const {
   n.size = 0;
   n.prev = 0;
   n.next = 0;
-  n.mem_ptr = NULL;
+  n.mem_ptr = nullptr;
   node_stack.push(indx);
 }
 
 template <class ValueType> inline Iterator<ValueType> aligned_new(Long n_elem, const MemoryManager* mem_mgr) {
-  if (!n_elem) return NULL;
+  if (!n_elem) return nullptr;
 
   static MemoryManager def_mem_mgr(0);
   if (!mem_mgr) mem_mgr = &def_mem_mgr;
   Iterator<ValueType> A = (Iterator<ValueType>)mem_mgr->malloc(n_elem, sizeof(ValueType));
-  SCTL_ASSERT_MSG(A != NULL, "memory allocation failed.");
+  SCTL_ASSERT_MSG(A != nullptr, "memory allocation failed.");
 
-  if (!TypeTraits<ValueType>::IsPOD()) {  // Call constructors
+  if (!std::is_trivial<ValueType>::value) {  // Call constructors
                                           // printf("%s\n", __PRETTY_FUNCTION__);
 #pragma omp parallel for schedule(static)
     for (Long i = 0; i < n_elem; i++) {
@@ -514,9 +497,9 @@ template <class ValueType> inline Iterator<ValueType> aligned_new(Long n_elem, c
 }
 
 template <class ValueType> inline void aligned_delete(Iterator<ValueType> A, const MemoryManager* mem_mgr) {
-  if (A == NULL) return;
+  if (A == nullptr) return;
 
-  if (!TypeTraits<ValueType>::IsPOD()) {  // Call destructors
+  if (!std::is_trivial<ValueType>::value) {  // Call destructors
     // printf("%s\n", __PRETTY_FUNCTION__);
     MemoryManager::MemHead& mem_head = MemoryManager::GetMemHead((char*)&A[0]);
 #ifdef SCTL_MEMDEBUG
@@ -547,13 +530,12 @@ template <class ValueType> inline void aligned_delete(Iterator<ValueType> A, con
 }
 
 template <class ValueType> inline Iterator<ValueType> memcopy(Iterator<ValueType> destination, ConstIterator<ValueType> source, Long num) {
-  static_assert(std::is_trivially_copyable<ValueType>::value, "Data is not trivially copyable!");
   if (destination != source && num) {
 #ifdef SCTL_MEMDEBUG
     destination[num - 1];
     source[num - 1];
 #endif
-    if (TypeTraits<ValueType>::IsPOD()) {
+    if (std::is_trivially_copyable<ValueType>::value) {
       memcpy(&destination[0], &source[0], num * sizeof(ValueType));
     } else {
       for (Long i = 0; i < num; i++) destination[i] = source[i];

include/sctl/parallel_solver.hpp

@@ -38,7 +38,7 @@ template <class Real> int ParallelSolverMatVec(Mat M_, Vec x_, Vec Mx_) {
   VecGetLocalSize(Mx_, &N_);
   SCTL_ASSERT(N == N_);
 
-  void* data = NULL;
+  void* data = nullptr;
   MatShellGetContext(M_, &data);
   auto& M = dynamic_cast<const typename ParallelSolver<Real>::ParallelOp&>(*(typename ParallelSolver<Real>::ParallelOp*)data);
 
@@ -101,7 +101,7 @@ template <class Real> inline void ParallelSolver<Real>::operator()(Vector<Real>*
   KSPSetType(ksp, KSPGMRES);
   KSPSetNormType(ksp, KSP_NORM_UNPRECONDITIONED);
   KSPSetTolerances(ksp, tol, PETSC_DEFAULT, PETSC_DEFAULT, max_iter);
-  if (verbose_) KSPMonitorSet(ksp, KSPMonitorDefault, NULL, NULL);
+  if (verbose_) KSPMonitorSet(ksp, KSPMonitorDefault, nullptr, nullptr);
   KSPGMRESSetRestart(ksp, max_iter);
   ierr = KSPSetFromOptions(ksp);
   CHKERRABORT(comm, ierr);
@@ -152,7 +152,7 @@ template <class Real> static Real inner_prod(const Vector<Real>& x, const Vector
   for (Long i = 0; i < N; i++) x_dot_y += x[i] * y[i];
 
   Real x_dot_y_glb = 0;
-  comm.Allreduce(SCTL_PTR2CONSTITR(Real, &x_dot_y, 1), SCTL_PTR2ITR(Real, &x_dot_y_glb, 1), 1, Comm::CommOp::SUM);
+  comm.Allreduce(Ptr2ConstItr<Real>(&x_dot_y, 1), Ptr2Itr<Real>(&x_dot_y_glb, 1), 1, Comm::CommOp::SUM);
 
   return x_dot_y_glb;
 }

include/sctl/profile.hpp

@@ -23,11 +23,11 @@ class Profile {
 
   static bool Enable(bool state);
 
-  static void Tic(const char* name_, const Comm* comm_ = NULL, bool sync_ = false, Integer level = 0);
+  static void Tic(const char* name_, const Comm* comm_ = nullptr, bool sync_ = false, Integer level = 0);
 
   static void Toc();
 
-  static void print(const Comm* comm_ = NULL);
+  static void print(const Comm* comm_ = nullptr);
 
   static void reset();
 

include/sctl/profile.txx

@@ -48,10 +48,10 @@ inline void Profile::Tic(const char* name_, const Comm* comm_, bool sync_, Integ
   // sync_=true;
   if (!prof.enable_state) return;
   if (verbose <= SCTL_PROFILE && prof.verb_level.size() == prof.enable_depth) {
-    if (comm_ != NULL && sync_) comm_->Barrier();
+    if (comm_ != nullptr && sync_) comm_->Barrier();
 #ifdef SCTL_VERBOSE
     Integer rank = 0;
-    if (comm_ != NULL) rank = comm_->Rank();
+    if (comm_ != nullptr) rank = comm_->Rank();
     if (!rank) {
       for (size_t i = 0; i < prof.name.size(); i++) std::cout << "    ";
       std::cout << "\033[1;31m" << name_ << "\033[0m {\n";
@@ -98,7 +98,7 @@ inline void Profile::Toc() {
     prof.max_m_log.push_back(prof.max_mem.back());
 
 #ifndef NDEBUG
-    if (comm_ != NULL && sync_) comm_->Barrier();
+    if (comm_ != nullptr && sync_) comm_->Barrier();
 #endif
     prof.name.pop();
     prof.comm.pop();
@@ -107,7 +107,7 @@ inline void Profile::Toc() {
 
 #ifdef SCTL_VERBOSE
     Integer rank = 0;
-    if (comm_ != NULL) rank = comm_->Rank();
+    if (comm_ != nullptr) rank = comm_->Rank();
     if (!rank) {
       for (size_t i = 0; i < prof.name.size(); i++) std::cout << "    ";
       std::cout << "}\n";
@@ -125,7 +125,7 @@ inline void Profile::print(const Comm* comm_) {
   SCTL_ASSERT_MSG(prof.name.empty(), "Missing balancing Toc()");
 
   Comm c_self = Comm::Self();
-  if (comm_ == NULL) comm_ = &c_self;
+  if (comm_ == nullptr) comm_ = &c_self;
   comm_->Barrier();
 
   Integer np, rank;

include/sctl/stacktrace.h

@@ -63,7 +63,7 @@ inline void abortHandler(int signum, siginfo_t* si, void* unused) {
 #pragma omp critical(STACK_TRACE)
   if (first_time) {
     first_time = false;
-    const char* name = NULL;
+    const char* name = nullptr;
     switch (signum) {
       case SIGABRT:
         name = "SIGABRT";
@@ -98,12 +98,12 @@ inline int SetSigHandler() {
   sa.sa_sigaction = abortHandler;
   sigemptyset(&sa.sa_mask);
 
-  sigaction(SIGABRT, &sa, NULL);
-  sigaction(SIGSEGV, &sa, NULL);
-  sigaction(SIGBUS, &sa, NULL);
-  sigaction(SIGILL, &sa, NULL);
-  sigaction(SIGFPE, &sa, NULL);
-  sigaction(SIGPIPE, &sa, NULL);
+  sigaction(SIGABRT, &sa, nullptr);
+  sigaction(SIGSEGV, &sa, nullptr);
+  sigaction(SIGBUS, &sa, nullptr);
+  sigaction(SIGILL, &sa, nullptr);
+  sigaction(SIGFPE, &sa, nullptr);
+  sigaction(SIGPIPE, &sa, nullptr);
 
   return 0;
 }

include/sctl/vector.hpp

@@ -10,12 +10,18 @@
 namespace SCTL_NAMESPACE {
 
 template <class ValueType> class Vector {
-
  public:
+  typedef ValueType value_type;
+  typedef ValueType& reference;
+  typedef const ValueType& const_reference;
+  typedef Iterator<ValueType> iterator;
+  typedef ConstIterator<ValueType> const_iterator;
+  typedef Long difference_type;
+  typedef Long size_type;
 
   Vector();
 
-  Vector(Long dim_, Iterator<ValueType> data_ = Iterator<ValueType>(NULL), bool own_data_ = true);
+  Vector(Long dim_, Iterator<ValueType> data_ = nullptr, bool own_data_ = true);
 
   Vector(const Vector& V);
 
@@ -25,7 +31,7 @@ template <class ValueType> class Vector {
 
   void Swap(Vector<ValueType>& v1);
 
-  void ReInit(Long dim_, Iterator<ValueType> data_ = NULL, bool own_data_ = true);
+  void ReInit(Long dim_, Iterator<ValueType> data_ = nullptr, bool own_data_ = true);
 
   void Write(const char* fname) const;
 

include/sctl/vector.txx

@@ -11,7 +11,7 @@ template <class ValueType> Vector<ValueType>::Vector() {
   dim = 0;
   capacity = 0;
   own_data = true;
-  data_ptr = NULL;
+  data_ptr = nullptr;
 }
 
 template <class ValueType> Vector<ValueType>::Vector(Long dim_, Iterator<ValueType> data_, bool own_data_) {
@@ -21,11 +21,11 @@ template <class ValueType> Vector<ValueType>::Vector(Long dim_, Iterator<ValueTy
   if (own_data) {
     if (dim > 0) {
       data_ptr = aligned_new<ValueType>(capacity);
-      if (data_ != NULL) {
+      if (data_ != nullptr) {
         memcopy(data_ptr, data_, dim);
       }
     } else
-      data_ptr = NULL;
+      data_ptr = nullptr;
   } else
     data_ptr = data_;
 }
@@ -38,7 +38,7 @@ template <class ValueType> Vector<ValueType>::Vector(const Vector<ValueType>& V)
     data_ptr = aligned_new<ValueType>(capacity);
     memcopy(data_ptr, V.data_ptr, dim);
   } else
-    data_ptr = NULL;
+    data_ptr = nullptr;
 }
 
 template <class ValueType> Vector<ValueType>::Vector(const std::vector<ValueType>& V) {
@@ -49,16 +49,16 @@ template <class ValueType> Vector<ValueType>::Vector(const std::vector<ValueType
     data_ptr = aligned_new<ValueType>(capacity);
     memcopy(data_ptr, Ptr2ConstItr<ValueType>(&V[0], V.size()), dim);
   } else
-    data_ptr = NULL;
+    data_ptr = nullptr;
 }
 
 template <class ValueType> Vector<ValueType>::~Vector() {
   if (own_data) {
-    if (data_ptr != NULL) {
+    if (data_ptr != nullptr) {
       aligned_delete(data_ptr);
     }
   }
-  data_ptr = NULL;
+  data_ptr = nullptr;
   capacity = 0;
   dim = 0;
 }
@@ -83,7 +83,7 @@ template <class ValueType> void Vector<ValueType>::Swap(Vector<ValueType>& v1) {
 template <class ValueType> void Vector<ValueType>::ReInit(Long dim_, Iterator<ValueType> data_, bool own_data_) {
   if (own_data_ && own_data && dim_ <= capacity) {
     dim = dim_;
-    if (data_ != NULL) {
+    if (data_ != nullptr) {
       memcopy(data_ptr, data_, dim);
     }
   } else {
@@ -94,7 +94,7 @@ template <class ValueType> void Vector<ValueType>::ReInit(Long dim_, Iterator<Va
 
 template <class ValueType> void Vector<ValueType>::Write(const char* fname) const {
   FILE* f1 = fopen(fname, "wb+");
-  if (f1 == NULL) {
+  if (f1 == nullptr) {
     std::cout << "Unable to open file for writing:" << fname << '\n';
     return;
   }
@@ -108,7 +108,7 @@ template <class ValueType> void Vector<ValueType>::Write(const char* fname) cons
 
 template <class ValueType> void Vector<ValueType>::Read(const char* fname) {
   FILE* f1 = fopen(fname, "r");
-  if (f1 == NULL) {
+  if (f1 == nullptr) {
     std::cout << "Unable to open file for reading:" << fname << '\n';
     return;
   }