BlockStructure.cpp

Go to the documentation of this file.

/* 
 * @BEGIN LICENSE
 *
 * Copyright (C) 2014-2015  Ward Poelmans
 *
 * This file is part of v2DM-DOCI.
 * 
 * v2DM-DOCI is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Foobar is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Foobar.  If not, see <http://www.gnu.org/licenses/>.
 *
 * @END LICENSE
 */

#include <assert.h>

#include "Matrix.h"
#include "Vector.h"
#include "BlockStructure.h"

using namespace doci2DM;

   template<class BlockType>
BlockStructure<BlockType>::BlockStructure(int nr)
{
   blocks.resize(nr);

   degen.resize(nr,0);
}

   template<class BlockType>
BlockStructure<BlockType>::BlockStructure(const BlockStructure<BlockType> &blockmat_copy)
{
   blocks.resize(blockmat_copy.blocks.size());

#pragma omp parallel for
   for(int i=0;i<blocks.size();i++)
      blocks[i].reset(new BlockType(blockmat_copy[i]));

   degen = blockmat_copy.degen;
}

   template<class BlockType>
BlockStructure<BlockType>::BlockStructure(BlockStructure<BlockType> &&blockmat_copy)
{
   blocks = std::move(blockmat_copy.blocks);

   degen = std::move(blockmat_copy.degen);
}

   template<class BlockType>
void BlockStructure<BlockType>::setDim(int block,int dim,int degeneracy)
{
   assert(block < blocks.size());

   this->degen[block] = degeneracy;

   blocks[block].reset(new BlockType(dim));
}

   template<class BlockType>
BlockType &BlockStructure<BlockType>::operator[](int block)
{
   return *blocks[block];
}

template<class BlockType>
const BlockType &BlockStructure<BlockType>::operator[](int block) const
{
   return *blocks[block];
}


   template<class BlockType>
BlockStructure<BlockType> &BlockStructure<BlockType>::operator=(const BlockStructure<BlockType> &blockmat_copy)
{
   blocks.resize(blockmat_copy.blocks.size());

#pragma omp parallel for
   for(int i=0;i<blocks.size();i++)
      blocks[i].reset(new BlockType(blockmat_copy[i]));

   return *this;
}

   template<class BlockType>
BlockStructure<BlockType> &BlockStructure<BlockType>::operator=(double a)
{
#pragma omp parallel for
   for(int i = 0;i < blocks.size();++i)
      *blocks[i] = a;

   return *this;
}

   template<class BlockType>
BlockStructure<BlockType> &BlockStructure<BlockType>::operator+=(const BlockStructure<BlockType> &blockmat_pl)
{
#pragma omp parallel for
   for(int i = 0;i < blocks.size();++i)
      *blocks[i] += blockmat_pl[i];

   return *this;
}

   template<class BlockType>
BlockStructure<BlockType> &BlockStructure<BlockType>::operator-=(const BlockStructure<BlockType> &blockmat_pl)
{
#pragma omp parallel for
   for(int i = 0;i < blocks.size();++i)
      *blocks[i] -= blockmat_pl[i];

   return *this;
}

   template<class BlockType>
BlockStructure<BlockType> &BlockStructure<BlockType>::daxpy(double alpha,const BlockStructure<BlockType> &blockmat_pl)
{
#pragma omp parallel for
   for(int i = 0;i < blocks.size();++i)
      blocks[i]->daxpy(alpha,blockmat_pl[i]);

   return *this;
}

   template<class BlockType>
BlockStructure<BlockType> &BlockStructure<BlockType>::operator/=(double c)
{
#pragma omp parallel for
   for(int i = 0;i < blocks.size();++i)
      *blocks[i] /= c;

   return *this;
}

   template<class BlockType>
BlockStructure<BlockType> &BlockStructure<BlockType>::operator*=(double c)
{
#pragma omp parallel for
   for(int i = 0;i < blocks.size();++i)
      *blocks[i] *= c;

   return *this;
}

namespace doci2DM {
   template<>
double &BlockStructure<Matrix>::operator()(int block,int i,int j)
{
   return (*blocks[block])(i,j);
}

template<>
double BlockStructure<Matrix>::operator()(int block,int i,int j) const
{
   return (*blocks[block])(i,j);
}

   template<>
double& BlockStructure<Vector>::operator()(int block,int index)
{
   return (*blocks[block])[index];
}

template<>
double BlockStructure<Vector>::operator()(int block,int index) const
{
   return (*blocks[block])[index];
}

}

template<class BlockType>
int BlockStructure<BlockType>::gnr() const
{
   return blocks.size();
}

template<class BlockType>
int BlockStructure<BlockType>::gdim(int i) const
{
   return blocks[i]->gn();
}

template<class BlockType>
int BlockStructure<BlockType>::gdeg(int i) const
{
   return degen[i];
}

template<class BlockType>
double BlockStructure<BlockType>::trace() const
{
   double ward = 0;

#pragma omp parallel for reduction(+:ward)
   for(int i = 0;i < blocks.size();++i)
      ward += degen[i]*blocks[i]->trace();

   return ward;
}

template<class BlockType>
double BlockStructure<BlockType>::ddot(const BlockStructure<BlockType> &blocks_in) const
{
   double ward = 0.0;

#pragma omp parallel for reduction(+:ward)
   for(int i = 0;i < blocks.size();i++)
      ward += degen[i]*blocks[i]->ddot(blocks_in[i]);

   return ward;
}

   template<class BlockType>
void BlockStructure<BlockType>::invert()
{
#pragma omp parallel for
   for(int i = 0;i < blocks.size();++i)
      blocks[i]->invert();
}

   template<class BlockType>
void BlockStructure<BlockType>::dscal(double alpha)
{
#pragma omp parallel for
   for(int i = 0;i < blocks.size();++i)
      blocks[i]->dscal(alpha);
}

   template<class BlockType>
void BlockStructure<BlockType>::fill_Random()
{
   for(int i = 0;i < blocks.size();++i)
      blocks[i]->fill_Random();
}

   template<class BlockType>
void BlockStructure<BlockType>::fill_Random(int seed)
{
   for(int i = 0;i < blocks.size();++i)
      blocks[i]->fill_Random(seed);
}

   template<class BlockType>
void BlockStructure<BlockType>::sqrt(int option)
{
#pragma omp parallel for
   for(int i = 0;i < blocks.size();++i)
      blocks[i]->sqrt(option);
}

   template<class BlockType>
void BlockStructure<BlockType>::L_map(const BlockStructure<BlockType> &map,const BlockStructure<BlockType> &object)
{
#pragma omp parallel for
   for(int i = 0;i < blocks.size();++i)
      blocks[i]->L_map(map[i],object[i]);
}

   template<class BlockType>
BlockStructure<BlockType> &BlockStructure<BlockType>::mprod(const BlockStructure<BlockType> &A, const BlockStructure<BlockType> &B)
{
#pragma omp parallel for
   for(int i = 0;i < blocks.size();++i)
      blocks[i]->mprod(A[i],B[i]);

   return *this;
}

   template<class BlockType>
void BlockStructure<BlockType>::symmetrize()
{
#pragma omp parallel for
   for(int i = 0;i < blocks.size();++i)
      blocks[i]->symmetrize();
}

namespace doci2DM {
template<>
void BlockStructure<Vector>::sort()
{
   for(unsigned int i = 0;i < blocks.size();++i)
      blocks[i]->sort();
}
}

   template<class BlockType>
void BlockStructure<BlockType>::sep_pm(BlockStructure<BlockType> &pos, BlockStructure<BlockType> &neg)
{
   for(unsigned int i=0;i<blocks.size();++i)
      blocks[i]->sep_pm(*pos.blocks[i], *neg.blocks[i]);
}


namespace doci2DM
{
   template<class MyBlockType>
   std::ostream &operator<< (std::ostream &output,const doci2DM::BlockStructure<MyBlockType> &blocks_p)
   {
      for(int i = 0;i < blocks_p.blocks.size();++i)
      {
         output << i << "\t" << blocks_p.blocks[i]->gn() << "\t" << blocks_p.degen[i] << std::endl;
         output << std::endl;

         output << *blocks_p.blocks[i] << std::endl;
      }

      return output;
   }

   // explicit instantiation must occur in namespace
   template class BlockStructure<doci2DM::Matrix>;
   template class BlockStructure<doci2DM::Vector>;

   // also instantiation for friend functions
   template std::ostream &operator<< (std::ostream &output,const doci2DM::BlockStructure<doci2DM::Matrix> &blocks_p);
   template std::ostream &operator<< (std::ostream &output,const doci2DM::BlockStructure<doci2DM::Vector> &blocks_p);
}


/* vim: set ts=3 sw=3 expandtab :*/