DSymmMatrixParam.h

#ifndef UTIL_D_SYMM_MATRIX_PARAM_H
#define UTIL_D_SYMM_MATRIX_PARAM_H

/*
* Util Package - C++ Utilities for Scientific Computation
*
* Copyright 2010 - 2017, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/

#include <util/param/Parameter.h>
#include <util/containers/DMatrix.h>
#include <util/format/Int.h>
#ifdef UTIL_MPI
#include <util/mpi/MpiSendRecv.h>
#endif
#include <util/global.h>

#include <iomanip> 

namespace Util
{

   template <class Type>
   class DSymmMatrixParam : public Parameter
   {
      
   public:
  
      DSymmMatrixParam(const char *label, DMatrix<Type>& matrix, int n, bool isRequired = true);
 
      void writeParam(std::ostream &out);

   protected:
      
      virtual void readValue(std::istream& in);

      virtual void loadValue(Serializable::IArchive& ar);

      virtual void saveValue(Serializable::OArchive& ar);

      #ifdef UTIL_MPI

      virtual void bcastValue();
      #endif

   private:
   
      DMatrix<Type>* matrixPtr_;
   
      int n_; 
   
   };

   /*
   * DMatrix constructor.
   */
   template <class Type>
   DSymmMatrixParam<Type>::DSymmMatrixParam(const char* label, DMatrix<Type>& matrix, int n, bool isRequired)
    : Parameter(label, isRequired),
      matrixPtr_(&matrix),
      n_(n)
   {}

   /*
   * Read a DMatrix from isteam.
   */
   template <class Type>
   void DSymmMatrixParam<Type>::readValue(std::istream &in)
   {  
      // Preconditions
      if (!(matrixPtr_->isAllocated())) {
         UTIL_THROW("Cannot read unallocated DMatrix");
      }
      if (n_ != matrixPtr_->capacity1()) {
         UTIL_THROW("Error: Logical size n_ != DMatrix<Type>::capacity1()");
      }
      if (n_ != matrixPtr_->capacity2()) {
         UTIL_THROW("Error: Logical size n_ != DMatrix<Type>::capacity2()");
      }

      // Create matrix of int/bool flags for error checking
      DMatrix<int> flags;
      flags.allocate(n_, n_);
      int i, j;
      for (i = 0; i < n_; ++i) {
         for (j = 0; j < n_; ++j) {
            flags(i,j) = 0;
         }
      }

      double value;
      for (int k = 0; k < n_*(n_ + 1)/2; ++k) {
         in >> i >> j >> value;
         UTIL_CHECK(flags(i,j) == 0);
         (*matrixPtr_)(i, j) = value;
         flags(i, j) = 1;
         if (i != j) {
            UTIL_CHECK(flags(j,i) == 0);
            (*matrixPtr_)(j, i) = value;
            flags(j, i) = 1;
         }
      }

   }

   /*
   * Load a DMatrix from input archive.
   */
   template <class Type>
   void DSymmMatrixParam<Type>::loadValue(Serializable::IArchive& ar)
   {  
      if (!(matrixPtr_->isAllocated())) {
         matrixPtr_->allocate(n_, n_);
      } else {
         if (n_ != matrixPtr_->capacity1()) {
            UTIL_THROW("Error: Logical size n_ != DMatrix<Type>::capacity1()");
         }
         if (n_ != matrixPtr_->capacity2()) {
            UTIL_THROW("Error: Logical size n_ != DMatrix<Type>::capacity2()");
         }
      }
      ar >> *matrixPtr_;
   }

   /*
   * Save a DMatrix to an output archive.
   */
   template <class Type>
   void DSymmMatrixParam<Type>::saveValue(Serializable::OArchive& ar)
   {
      if (n_ != matrixPtr_->capacity1()) {
         UTIL_THROW("Error: Logical size n_ != DMatrix<Type>::capacity1()");
      }
      if (n_ != matrixPtr_->capacity2()) {
         UTIL_THROW("Error: Logical size n_ != DMatrix<Type>::capacity2()");
      }
      ar << *matrixPtr_; 
   }

   #ifdef UTIL_MPI
   /*
   * Broadcast a DMatrix.
   */
   template <class Type>
   void DSymmMatrixParam<Type>::bcastValue()
   {  
      if (!(matrixPtr_->isAllocated())) {
         matrixPtr_->allocate(n_, n_);
      } else {
         if (n_ != matrixPtr_->capacity1()) {
            UTIL_THROW("Error: Logical size n_ > DMatrix<Type>::capacity1()");
         }
         if (n_ != matrixPtr_->capacity2()) {
            UTIL_THROW("Error: Logical size n_ > DMatrix<Type>::capacity2()");
         }
      }
      bcast<Type>(ioCommunicator(), *matrixPtr_, n_, n_, 0); 
   }
   #endif

   /*
   * Write a DSymmMatrixParam.
   */
   template <class Type>
   void DSymmMatrixParam<Type>::writeParam(std::ostream &out)
   {
      if (isActive()) {
         // Preconditions
         if (!(matrixPtr_->isAllocated())) {
            UTIL_THROW("Cannot read unallocated DMatrix");
         }
         if (n_ > matrixPtr_->capacity1()) {
            UTIL_THROW("Error: Logical size n_ > DMatrix<Type>::capacity1()");
         }
         if (n_ > matrixPtr_->capacity2()) {
            UTIL_THROW("Error: Logical size n_ > DMatrix<Type>::capacity2()");
         }
   
         Label space("");
         int i, j;
         for (i = 0; i < n_; ++i) {
            for (j = 0; j <= i; ++j) {
               if (i == 0 && j == 0) {
                  out << indent() << label_;
               } else {
                  out << indent() << space;
               }
               out << Int(i, 4) << " " << Int(j, 4) << " "
                   << std::right << std::scientific 
                   << std::setprecision(Parameter::Precision) 
                   << std::setw(Parameter::Width)
                   << (*matrixPtr_)(i, j)
                   << std::endl;
            }
         }
      }
   }

} 
#endif