elemop.cc

大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CHARMM并行计算。非常具有应用价值。

//
// elemop.cc
//
// Copyright (C) 1996 Limit Point Systems, Inc.
//
// Author: Curtis Janssen <[e-mail]>
// Maintainer: LPS
//
// This file is part of the SC Toolkit.
//
// The SC Toolkit is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
//
// The SC Toolkit 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 Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
//
// The U.S. Government is granted a limited license as per AL 91-7.
//

#ifdef __GNUC__
#pragma implementation
#endif

#include <stdexcept>

#include <stdlib.h>
#include <math.h>

#include <util/misc/formio.h>
#include <util/state/stateio.h>
#include <math/scmat/block.h>
#include <math/scmat/blkiter.h>
#include <math/scmat/elemop.h>
#include <math/scmat/abstract.h>

using namespace sc;

/////////////////////////////////////////////////////////////////////////////
// SCElementOp member functions

static ClassDesc SCElementOp_cd(
  typeid(SCElementOp),"SCElementOp",1,"public SavableState",
  0, 0, 0);

SCElementOp::SCElementOp()
{
}

SCElementOp::~SCElementOp()
{
}

int
SCElementOp::has_collect()
{
  return 0;
}

void
SCElementOp::defer_collect(int)
{
}

int
SCElementOp::has_side_effects()
{
  return 0;
}

void
SCElementOp::collect(const Ref<MessageGrp>&)
{
}

bool
SCElementOp::threadsafe()
{
  return false;
}

bool
SCElementOp::cloneable()
{
  return false;
}

Ref<SCElementOp>
SCElementOp::clone()
{
  throw std::runtime_error("SCElementOp::clone: not implemented");
}

void
SCElementOp::collect(const Ref<SCElementOp> &)
{
  throw std::runtime_error("SCElementOp::collect(const Ref<SCElementOp> &): "
                           "not implemented");
}

void
SCElementOp::process_base(SCMatrixBlock* a)
{
  if (dynamic_cast<SCMatrixRectBlock*>(a))
    process_spec_rect(dynamic_cast<SCMatrixRectBlock*>(a));
  else if (dynamic_cast<SCMatrixLTriBlock*>(a))
    process_spec_ltri(dynamic_cast<SCMatrixLTriBlock*>(a));
  else if (dynamic_cast<SCMatrixDiagBlock*>(a))
    process_spec_diag(dynamic_cast<SCMatrixDiagBlock*>(a));
  else if (dynamic_cast<SCVectorSimpleBlock*>(a))
    process_spec_vsimp(dynamic_cast<SCVectorSimpleBlock*>(a));
  else if (dynamic_cast<SCMatrixRectSubBlock*>(a))
    process_spec_rectsub(dynamic_cast<SCMatrixRectSubBlock*>(a));
  else if (dynamic_cast<SCMatrixLTriSubBlock*>(a))
    process_spec_ltrisub(dynamic_cast<SCMatrixLTriSubBlock*>(a));
  else if (dynamic_cast<SCMatrixDiagSubBlock*>(a))
    process_spec_diagsub(dynamic_cast<SCMatrixDiagSubBlock*>(a));
  else if (dynamic_cast<SCVectorSimpleSubBlock*>(a))
    process_spec_vsimpsub(dynamic_cast<SCVectorSimpleSubBlock*>(a));
  else
    a->process(this);
}

// If specializations of SCElementOp do not handle a particle
// block type, then these functions will be called and will
// set up an appropiate block iterator which specializations
// of SCElementOp must handle since it is pure virtual.

void
SCElementOp::process_spec_rect(SCMatrixRectBlock* a)
{
  SCMatrixBlockIter*i = new SCMatrixRectBlockIter(a);
  SCMatrixBlockIter&r=*i;
  process(r);
  // this causes a SCMatrixRectBlock::operator int() to be
  // called with this = 0x0 using gcc 2.5.6
  // process(*i,b);
  delete i;
}
void
SCElementOp::process_spec_ltri(SCMatrixLTriBlock* a)
{
  SCMatrixBlockIter*i = new SCMatrixLTriBlockIter(a);
  process(*i);
  delete i;
}
void
SCElementOp::process_spec_diag(SCMatrixDiagBlock* a)
{
  SCMatrixBlockIter*i = new SCMatrixDiagBlockIter(a);
  process(*i);
  delete i;
}
void
SCElementOp::process_spec_vsimp(SCVectorSimpleBlock* a)
{
  SCMatrixBlockIter*i = new SCVectorSimpleBlockIter(a);
  process(*i);
  delete i;
}
void
SCElementOp::process_spec_rectsub(SCMatrixRectSubBlock* a)
{
  SCMatrixBlockIter*i = new SCMatrixRectSubBlockIter(a);
  SCMatrixBlockIter&r=*i;
  process(r);
  // this causes a SCMatrixRectBlock::operator int() to be
  // called with this = 0x0 using gcc 2.5.6
  // process(*i,b);
  delete i;
}
void
SCElementOp::process_spec_ltrisub(SCMatrixLTriSubBlock* a)
{
  SCMatrixBlockIter*i = new SCMatrixLTriSubBlockIter(a);
  process(*i);
  delete i;
}
void
SCElementOp::process_spec_diagsub(SCMatrixDiagSubBlock* a)
{
  SCMatrixBlockIter*i = new SCMatrixDiagSubBlockIter(a);
  process(*i);
  delete i;
}
void
SCElementOp::process_spec_vsimpsub(SCVectorSimpleSubBlock* a)
{
  SCMatrixBlockIter*i = new SCVectorSimpleSubBlockIter(a);
  process(*i);
  delete i;
}

/////////////////////////////////////////////////////////////////////////////
// SCElementOp2 member functions

static ClassDesc SCElementOp2_cd(
  typeid(SCElementOp2),"SCElementOp2",1,"public SavableState",
  0, 0, 0);

SCElementOp2::SCElementOp2()
{
}

SCElementOp2::~SCElementOp2()
{
}

int
SCElementOp2::has_collect()
{
  return 0;
}

void
SCElementOp2::defer_collect(int)
{
}

int
SCElementOp2::has_side_effects()
{
  return 0;
}

int
SCElementOp2::has_side_effects_in_arg()
{
  return 0;
}

void
SCElementOp2::collect(const Ref<MessageGrp>&)
{
}

void
SCElementOp2::process_base(SCMatrixBlock* a, SCMatrixBlock* b)
{
  a->process(this, b);
}

// If specializations of SCElementOp2 do not handle a particle
// block type, then these functions will be called and will
// set up an appropiate block iterator which specializations
// of SCElementOp2 must handle since it is pure virtual.

void
SCElementOp2::process_spec_rect(SCMatrixRectBlock* a,SCMatrixRectBlock* b)
{
  SCMatrixBlockIter*i = new SCMatrixRectBlockIter(a);
  SCMatrixBlockIter*j = new SCMatrixRectBlockIter(b);
  process(*i,*j);
  // this causes a SCMatrixRectBlock::operator int() to be
  // called with this = 0x0 using gcc 2.5.6
  // process(*i,b);
  delete i;
  delete j;
}
void
SCElementOp2::process_spec_ltri(SCMatrixLTriBlock* a,SCMatrixLTriBlock* b)
{
  SCMatrixBlockIter*i = new SCMatrixLTriBlockIter(a);
  SCMatrixBlockIter*j = new SCMatrixLTriBlockIter(b);
  process(*i,*j);
  delete i;
  delete j;
}
void
SCElementOp2::process_spec_diag(SCMatrixDiagBlock* a,SCMatrixDiagBlock* b)
{
  SCMatrixBlockIter*i = new SCMatrixDiagBlockIter(a);
  SCMatrixBlockIter*j = new SCMatrixDiagBlockIter(b);
  process(*i,*j);
  delete i;
  delete j;
}
void
SCElementOp2::process_spec_vsimp(SCVectorSimpleBlock* a,SCVectorSimpleBlock* b)
{
  SCMatrixBlockIter*i = new SCVectorSimpleBlockIter(a);
  SCMatrixBlockIter*j = new SCVectorSimpleBlockIter(b);
  process(*i,*j);
  delete i;
  delete j;
}

/////////////////////////////////////////////////////////////////////////////
// SCElementOp3 member functions

static ClassDesc SCElementOp3_cd(
  typeid(SCElementOp3),"SCElementOp3",1,"public SavableState",
  0, 0, 0);

SCElementOp3::SCElementOp3()
{
}

SCElementOp3::~SCElementOp3()
{
}

int
SCElementOp3::has_collect()
{
  return 0;
}

void
SCElementOp3::defer_collect(int)
{
}

int
SCElementOp3::has_side_effects()
{
  return 0;
}

int
SCElementOp3::has_side_effects_in_arg1()
{
  return 0;
}

int
SCElementOp3::has_side_effects_in_arg2()
{
  return 0;
}

void
SCElementOp3::collect(const Ref<MessageGrp>&)
{
}

void
SCElementOp3::process_base(SCMatrixBlock* a,
                           SCMatrixBlock* b,
                           SCMatrixBlock* c)
{
  a->process(this, b, c);
}

// If specializations of SCElementOp3 do not handle a particle
// block type, then these functions will be called and will
// set up an appropiate block iterator which specializations
// of SCElementOp3 must handle since it is pure virtual.

void
SCElementOp3::process_spec_rect(SCMatrixRectBlock* a,
                                SCMatrixRectBlock* b,
                                SCMatrixRectBlock* c)
{
  SCMatrixBlockIter*i = new SCMatrixRectBlockIter(a);
  SCMatrixBlockIter*j = new SCMatrixRectBlockIter(b);
  SCMatrixBlockIter*k = new SCMatrixRectBlockIter(c);
  process(*i,*j,*k);
  delete i;
  delete j;
  delete k;
}
void
SCElementOp3::process_spec_ltri(SCMatrixLTriBlock* a,
                                SCMatrixLTriBlock* b,
                                SCMatrixLTriBlock* c)
{
  SCMatrixBlockIter*i = new SCMatrixLTriBlockIter(a);
  SCMatrixBlockIter*j = new SCMatrixLTriBlockIter(b);
  SCMatrixBlockIter*k = new SCMatrixLTriBlockIter(c);
  process(*i,*j,*k);
  delete i;
  delete j;
  delete k;
}
void
SCElementOp3::process_spec_diag(SCMatrixDiagBlock* a,
                                SCMatrixDiagBlock* b,
                                SCMatrixDiagBlock* c)
{
  SCMatrixBlockIter*i = new SCMatrixDiagBlockIter(a);
  SCMatrixBlockIter*j = new SCMatrixDiagBlockIter(b);
  SCMatrixBlockIter*k = new SCMatrixDiagBlockIter(c);
  process(*i,*j,*k);
  delete i;
  delete j;
  delete k;
}
void
SCElementOp3::process_spec_vsimp(SCVectorSimpleBlock* a,
                                 SCVectorSimpleBlock* b,
                                 SCVectorSimpleBlock* c)
{
  SCMatrixBlockIter*i = new SCVectorSimpleBlockIter(a);
  SCMatrixBlockIter*j = new SCVectorSimpleBlockIter(b);
  SCMatrixBlockIter*k = new SCVectorSimpleBlockIter(c);
  process(*i,*j,*k);
  delete i;
  delete j;
  delete k;
}

/////////////////////////////////////////////////////////////////////////
// SCElementScale members

static ClassDesc SCElementScale_cd(
  typeid(SCElementScale),"SCElementScale",1,"public SCElementOp",
  0, 0, create<SCElementScale>);
SCElementScale::SCElementScale(double a):scale(a) {}
SCElementScale::SCElementScale(StateIn&s):
  SCElementOp(s)
{
  s.get(scale);
}
void
SCElementScale::save_data_state(StateOut&s)
{
  s.put(scale);
}
SCElementScale::~SCElementScale() {}
void
SCElementScale::process(SCMatrixBlockIter&i)
{
  for (i.reset(); i; ++i) {
      i.set(scale*i.get());
    }
}

int
SCElementScale::has_side_effects()
{
  return 1;
}

/////////////////////////////////////////////////////////////////////////
// SCElementScalarProduct members

static ClassDesc SCElementScalarProduct_cd(
  typeid(SCElementScalarProduct),"SCElementScalarProduct",1,"public SCElementOp2",
  0, 0, create<SCElementScalarProduct>);

SCElementScalarProduct::SCElementScalarProduct():
  deferred_(0), product(0.0)
{
}

SCElementScalarProduct::SCElementScalarProduct(StateIn&s):
  SCElementOp2(s)
{
  s.get(product);
  s.get(deferred_);
}

void
SCElementScalarProduct::save_data_state(StateOut&s)
{
  s.put(product);
  s.put(deferred_);
}

SCElementScalarProduct::~SCElementScalarProduct()
{
}

void
SCElementScalarProduct::process(SCMatrixBlockIter&i,
                                SCMatrixBlockIter&j)
{
  for (i.reset(),j.reset(); i; ++i,++j) {
      product += i.get()*j.get();
    }
}

int
SCElementScalarProduct::has_collect()
{
  return 1;
}

void
SCElementScalarProduct::defer_collect(int h)
{
  deferred_=h;
}

void
SCElementScalarProduct::collect(const Ref<MessageGrp>&grp)
{
  if (!deferred_)
    grp->sum(product);
}

double
SCElementScalarProduct::result()
{
  return product;
}

/////////////////////////////////////////////////////////////////////////
// SCDestructiveElementProduct members

static ClassDesc SCDestructiveElementProduct_cd(
  typeid(SCDestructiveElementProduct),"SCDestructiveElementProduct",1,"public SCElementOp2",
  0, 0, create<SCDestructiveElementProduct>);
SCDestructiveElementProduct::SCDestructiveElementProduct() {}
SCDestructiveElementProduct::SCDestructiveElementProduct(StateIn&s):
  SCElementOp2(s)
{
}
void
SCDestructiveElementProduct::save_data_state(StateOut&s)
{
}
SCDestructiveElementProduct::~SCDestructiveElementProduct() {}
void
SCDestructiveElementProduct::process(SCMatrixBlockIter&i,
                                     SCMatrixBlockIter&j)
{
  for (i.reset(),j.reset(); i; ++i,++j) {
      i.set(i.get()*j.get());
    }
}

int
SCDestructiveElementProduct::has_side_effects()
{
  return 1;