blockedrect.cc

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

//
// blockedrect.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.
//

#include <math.h>

#include <util/misc/formio.h>
#include <util/keyval/keyval.h>
#include <util/state/stateio.h>
#include <math/scmat/blocked.h>
#include <math/scmat/cmatrix.h>
#include <math/scmat/elemop.h>

#include <math/scmat/local.h>
#include <math/scmat/repl.h>

using namespace std;
using namespace sc;

/////////////////////////////////////////////////////////////////////////////
// BlockedSCMatrix member functions

static ClassDesc BlockedSCMatrix_cd(
  typeid(BlockedSCMatrix),"BlockedSCMatrix",1,"public SCMatrix",
  0, 0, 0);

void
BlockedSCMatrix::resize(SCDimension *a, SCDimension *b)
{
  if (mats_) {
    delete[] mats_;
    mats_=0;
  }
  
  d1 = a;
  d2 = b;
  
  if (!a || !b || !a->blocks()->nblock() || !b->blocks()->nblock())
    return;

  if (a->blocks()->nblock() > 1 && b->blocks()->nblock() == 1) {
    nblocks_ = d1->blocks()->nblock();
    mats_ = new RefSCMatrix[d1->blocks()->nblock()];
    for (int i=0; i < d1->blocks()->nblock(); i++)
      if (d1->blocks()->size(i) && d2->blocks()->size(0))
        mats_[i] = subkit->matrix(d1->blocks()->subdim(i),
                                  d2->blocks()->subdim(0));

  } else if (a->blocks()->nblock() == 1 && b->blocks()->nblock() > 1) {
    nblocks_ = d2->blocks()->nblock();
    mats_ = new RefSCMatrix[d2->blocks()->nblock()];
    for (int i=0; i < d2->blocks()->nblock(); i++)
      if (d2->blocks()->size(i) && d1->blocks()->size(0))
        mats_[i] = subkit->matrix(d1->blocks()->subdim(0),
                                  d2->blocks()->subdim(i));

  } else if (a->blocks()->nblock() == b->blocks()->nblock()) {
    nblocks_ = d2->blocks()->nblock();
    mats_ = new RefSCMatrix[d1->blocks()->nblock()];
    for (int i=0; i < d1->blocks()->nblock(); i++)
      if (d2->blocks()->size(i) && d1->blocks()->size(i))
        mats_[i] = subkit->matrix(d1->blocks()->subdim(i),
                                  d2->blocks()->subdim(i));

  } else {
    ExEnv::errn() << indent << "BlockedSCMatrix::resize: wrong number of blocks\n";
    abort();
  }

}

BlockedSCMatrix::BlockedSCMatrix(const RefSCDimension&a,
                                 const RefSCDimension&b,
                                 BlockedSCMatrixKit*k):
  SCMatrix(a,b,k),
  subkit(k->subkit()),
  mats_(0), nblocks_(0)
{
  resize(a,b);
}

BlockedSCMatrix::~BlockedSCMatrix()
{
  if (mats_) {
    delete[] mats_;
    mats_=0;
  }
  nblocks_=0;
}

void
BlockedSCMatrix::assign_val(double v)
{
  for (int i=0; i < nblocks_; i++)
    if (mats_[i].nonnull())
      mats_[i]->assign(v);
}

double
BlockedSCMatrix::get_element(int i,int j) const
{
  int block_i, block_j;
  int elem_i, elem_j;

  d1->blocks()->elem_to_block(i,block_i,elem_i);
  d2->blocks()->elem_to_block(j,block_j,elem_j);

  if (d1->blocks()->nblock() == 1 && d2->blocks()->nblock() > 1) {
    return mats_[block_j]->get_element(elem_i,elem_j);

  } else if (d1->blocks()->nblock() > 1 && d2->blocks()->nblock() == 1) {
    return mats_[block_i]->get_element(elem_i,elem_j);
  } else if (d1->blocks()->nblock() == d2->blocks()->nblock()
             && block_i == block_j) {
    return mats_[block_i]->get_element(elem_i,elem_j);
  } else {
    return 0;
  }
}

void
BlockedSCMatrix::set_element(int i,int j,double a)
{
  int block_i, block_j;
  int elem_i, elem_j;

  d1->blocks()->elem_to_block(i,block_i,elem_i);
  d2->blocks()->elem_to_block(j,block_j,elem_j);
  
  if (d1->blocks()->nblock() == 1 && d2->blocks()->nblock() > 1) {
    mats_[block_j]->set_element(elem_i,elem_j,a);

  } else if (d1->blocks()->nblock() > 1 && d2->blocks()->nblock()) {
    mats_[block_i]->set_element(elem_i,elem_j,a);
  } else if (d1->blocks()->nblock() == d2->blocks()->nblock()
             && block_i == block_j) {
    mats_[block_i]->set_element(elem_i,elem_j,a);
  }
}

void
BlockedSCMatrix::accumulate_element(int i,int j,double a)
{
  int block_i, block_j;
  int elem_i, elem_j;

  d1->blocks()->elem_to_block(i,block_i,elem_i);
  d2->blocks()->elem_to_block(j,block_j,elem_j);
  
  if (d1->blocks()->nblock() == 1 && d2->blocks()->nblock() > 1) {
    mats_[block_j]->accumulate_element(elem_i,elem_j,a);

  } else if (d1->blocks()->nblock() > 1 && d2->blocks()->nblock() == 1 ||
             d1->blocks()->nblock() == d2->blocks()->nblock()) {
    mats_[block_i]->accumulate_element(elem_i,elem_j,a);
  }
}

SCMatrix *
BlockedSCMatrix::get_subblock(int br, int er, int bc, int ec)
{
  ExEnv::errn() << indent << "BlockedSCMatrix::get_subblock: cannot get subblock\n";
  abort();
  return 0;
}

void
BlockedSCMatrix::assign_subblock(SCMatrix*sb, int br, int er, int bc, int ec,
                               int source_br, int source_bc)
{
  ExEnv::errn() << indent
       << "BlockedSCMatrix::assign_subblock: cannot assign subblock\n";
  abort();
}

void
BlockedSCMatrix::accumulate_subblock(SCMatrix*sb,
                                     int br, int er, int bc, int ec,
                                     int source_br, int source_bc)
{
  ExEnv::errn() << indent << "BlockedSCMatrix::accumulate_subblock:"
       << " cannot accumulate subblock\n";
  abort();
}

SCVector *
BlockedSCMatrix::get_row(int i)
{
  ExEnv::errn() << indent << "BlockedSCMatrix::get_row: cannot get row\n";
  abort();

  return 0;
}

void
BlockedSCMatrix::assign_row(SCVector *v, int i)
{
  ExEnv::errn() << indent << "BlockedSCMatrix::assign_row: cannot assign row\n";
  abort();
}

void
BlockedSCMatrix::accumulate_row(SCVector *v, int i)
{
  ExEnv::errn() << indent << "BlockedSCMatrix::accumulate_row: cannot accumulate row\n";
  abort();
}

SCVector *
BlockedSCMatrix::get_column(int i)
{
  ExEnv::errn() << indent << "BlockedSCMatrix::get_column: cannot get column\n";
  abort();

  return 0;
}

void
BlockedSCMatrix::assign_column(SCVector *v, int i)
{
  ExEnv::errn() << indent << "BlockedSCMatrix::assign_column: cannot assign column\n";
  abort();
}

void
BlockedSCMatrix::accumulate_column(SCVector *v, int i)
{
  ExEnv::errn() << indent
       << "BlockedSCMatrix::accumulate_column: cannot accumulate column\n";
  abort();
}

// does the outer product a x b.  this must have rowdim() == a->dim() and
// coldim() == b->dim()
void
BlockedSCMatrix::accumulate_outer_product(SCVector*a,SCVector*b)
{
  const char* name = "BlockedSCMatrix::accumulate_outer_product";
  // make sure that the arguments are of the correct type
  BlockedSCVector* la = require_dynamic_cast<BlockedSCVector*>(a,name);
  BlockedSCVector* lb = require_dynamic_cast<BlockedSCVector*>(b,name);

  // make sure that the dimensions match
  if (!rowdim()->equiv(la->dim()) || !coldim()->equiv(lb->dim())) {
    ExEnv::errn() << indent
         << "BlockedSCMatrix::accumulate_outer_product(SCVector*,SCVector*): "
         << "dimensions don't match\n";
    abort();
  }

  for (int i=0; i < d1->blocks()->nblock(); i++)
    if (mats_[i].nonnull())
      mats_[i]->accumulate_outer_product(la->vecs_[i], lb->vecs_[i]);
}

void
BlockedSCMatrix::accumulate_product_rr(SCMatrix*a,SCMatrix*b)
{
  int i, zero = 0;

  const char* name = "BlockedSCMatrix::accumulate_product";
  // make sure that the arguments are of the correct type
  BlockedSCMatrix* la = require_dynamic_cast<BlockedSCMatrix*>(a,name);
  BlockedSCMatrix* lb = require_dynamic_cast<BlockedSCMatrix*>(b,name);

  // make sure that the dimensions match
  if (!rowdim()->equiv(la->rowdim()) || !coldim()->equiv(lb->coldim()) ||
      !la->coldim()->equiv(lb->rowdim())) {
    ExEnv::errn() << indent
         << "BlockedSCMatrix::accumulate_product_rr(SCMatrix*a,SCMatrix*b): "
         << "dimensions don't match\n";
    abort();
  }

  // find out the number of blocks we need to process.
  int mxnb = (nblocks_ > la->nblocks_) ? nblocks_ : la->nblocks_;
  
  int nrba = la->d1->blocks()->nblock();
  int ncba = la->d2->blocks()->nblock();
  int nrbb = lb->d1->blocks()->nblock();
  int ncbb = lb->d2->blocks()->nblock();
  
  int &mi = (nrba==1 && ncba > 1 && nrbb > 1 && ncbb==1) ? zero : i;
  int &ai = (nrba==1 && ncba==1) ? zero : i;
  int &bi = (nrbb==1 && ncbb==1) ? zero : i;

  for (i=0; i < mxnb; i++) {
    if (mats_[mi].null() || la->mats_[ai].null() || lb->mats_[bi].null())
      continue;
    mats_[mi]->accumulate_product(la->mats_[ai], lb->mats_[bi]);
  }
}

void
BlockedSCMatrix::accumulate_product_rs(SCMatrix*a,SymmSCMatrix*b)
{
  int i, zero=0;
  
  const char* name = "BlockedSCMatrix::accumulate_product";
  // make sure that the arguments are of the correct type
  BlockedSCMatrix* la = require_dynamic_cast<BlockedSCMatrix*>(a,name);
  BlockedSymmSCMatrix* lb = require_dynamic_cast<BlockedSymmSCMatrix*>(b,name);

  // make sure that the dimensions match
  if (!rowdim()->equiv(la->rowdim()) || !coldim()->equiv(lb->dim()) ||
      !la->coldim()->equiv(lb->dim())) {
    ExEnv::errn() << indent
         << "BlockedSCMatrix::accumulate_product_rs(SCMatrix*a,SymmSCMatrix*b): "
         << "dimensions don't match\n";
    abort();
  }

  int &bi = (lb->d->blocks()->nblock()==1) ? zero : i;
  
  for (i=0; i < nblocks_; i++) {
    if (mats_[i].null() || la->mats_[i].null() || lb->mats_[bi].null())
      continue;
    mats_[i]->accumulate_product(la->mats_[i], lb->mats_[bi]);
  }
}


void
BlockedSCMatrix::accumulate_product_rd(SCMatrix*a,DiagSCMatrix*b)
{
  int i, zero=0;
  
  const char* name = "BlockedSCMatrix::accumulate_product";
  // make sure that the arguments are of the correct type
  BlockedSCMatrix* la = require_dynamic_cast<BlockedSCMatrix*>(a,name);
  BlockedDiagSCMatrix* lb = require_dynamic_cast<BlockedDiagSCMatrix*>(b,name);

  // make sure that the dimensions match
  if (!rowdim()->equiv(la->rowdim()) || !coldim()->equiv(lb->dim()) ||
      !la->coldim()->equiv(lb->dim())) {
    ExEnv::errn() << indent
         << "BlockedSCMatrix::accumulate_product_rd(SCMatrix*a,DiagSCMatrix*b): "
         << "dimensions don't match\n";
    abort();
  }

  int &bi = (lb->d->blocks()->nblock()==1) ? zero : i;
  
  for (i=0; i < nblocks_; i++) {
    if (mats_[i].null() || la->mats_[i].null() || lb->mats_[bi].null())
      continue;
    mats_[i]->accumulate_product(la->mats_[i], lb->mats_[bi]);
  }
}

void
BlockedSCMatrix::accumulate(const SCMatrix*a)
{
  // make sure that the arguments is of the correct type
  const BlockedSCMatrix* la
    = require_dynamic_cast<const BlockedSCMatrix*>(a,"BlockedSCMatrix::accumulate");

  // make sure that the dimensions match
  if (!rowdim()->equiv(la->rowdim()) || !coldim()->equiv(la->coldim())) {
    ExEnv::errn() << indent << "BlockedSCMatrix::accumulate(SCMatrix*a): "
         << "dimensions don't match\n";
    abort();
  }

  for (int i=0; i < nblocks_; i++)
    if (mats_[i].nonnull())
      mats_[i]->accumulate(la->mats_[i]);
}

void
BlockedSCMatrix::accumulate(const SymmSCMatrix*a)
{
  // make sure that the arguments is of the correct type
  const BlockedSymmSCMatrix* la
    = require_dynamic_cast<const BlockedSymmSCMatrix*>(a,"BlockedSCMatrix::accumulate");

  // make sure that the dimensions match
  if (!rowdim()->equiv(la->dim()) || !coldim()->equiv(la->dim())) {
    ExEnv::errn() << indent << "BlockedSCMatrix::accumulate(SymmSCMatrix*a): "
         << "dimensions don't match\n";
    abort();
  }

  for (int i=0; i < nblocks_; i++)
    if (mats_[i].nonnull())
      mats_[i]->accumulate(la->mats_[i]);
}