distrect.cc

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

//
// distrect.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 <iostream>
#include <iomanip>
#include <stdlib.h>
#include <math.h>

#include <util/misc/formio.h>
#include <util/keyval/keyval.h>
#include <math/scmat/dist.h>
#include <math/scmat/cmatrix.h>
#include <math/scmat/elemop.h>

using namespace std;
using namespace sc;

#define DEBUG 0

/////////////////////////////////////////////////////////////////////////////

static void
fail(const char *m)
{
  ExEnv::errn() << indent << "distrect.cc: error: " << m << endl;
  abort();
}

/////////////////////////////////////////////////////////////////////////////
// DistSCMatrix member functions

static ClassDesc DistSCMatrix_cd(
  typeid(DistSCMatrix),"DistSCMatrix",1,"public SCMatrix",
  0, 0, 0);

DistSCMatrix::DistSCMatrix(const RefSCDimension&a,const RefSCDimension&b,
                           DistSCMatrixKit*k):
  SCMatrix(a,b,k)
{
  init_blocklist();
}

int
DistSCMatrix::block_to_node(int i, int j) const
{
  return (i*d2->blocks()->nblock() + j)%messagegrp()->n();
}

Ref<SCMatrixBlock>
DistSCMatrix::block_to_block(int i, int j) const
{
  int offset = i*d2->blocks()->nblock() + j;
  int nproc = messagegrp()->n();

  if ((offset%nproc) != messagegrp()->me()) return 0;

  SCMatrixBlockListIter I;
  for (I=blocklist->begin(); I!=blocklist->end(); I++) {
      if (I.block()->blocki == i && I.block()->blockj == j)
          return I.block();
    }

  ExEnv::errn() << indent << "DistSCMatrix::block_to_block: internal error" << endl;
  abort();
  return 0;
}

double *
DistSCMatrix::find_element(int i, int j) const
{
  int bi, oi;
  d1->blocks()->elem_to_block(i, bi, oi);

  int bj, oj;
  d2->blocks()->elem_to_block(j, bj, oj);

  Ref<SCMatrixRectBlock> blk
      = dynamic_cast<SCMatrixRectBlock*>(block_to_block(bi, bj).pointer());
  if (blk.nonnull()) {
      return &blk->data[oi*(blk->jend-blk->jstart)+oj];
    }
  else {
      return 0;
    }
}

int
DistSCMatrix::element_to_node(int i, int j) const
{
  int bi, oi;
  d1->blocks()->elem_to_block(i, bi, oi);

  int bj, oj;
  d2->blocks()->elem_to_block(j, bj, oj);

  return block_to_node(bi,bj);
}

void
DistSCMatrix::init_blocklist()
{
  int i, j, index;
  int me = messagegrp()->me();
  blocklist = new SCMatrixBlockList;
  for (i=0, index=0; i<d1->blocks()->nblock(); i++) {
      for (j=0; j<d2->blocks()->nblock(); j++, index++) {
          if (block_to_node(i,j) == me) {
              Ref<SCMatrixBlock> b
                  = new SCMatrixRectBlock(d1->blocks()->start(i),
                                          d1->blocks()->fence(i),
                                          d2->blocks()->start(j),
                                          d2->blocks()->fence(j));
              b->blocki = i;
              b->blockj = j;
              blocklist->append(b);
            }
        }
    }
}

DistSCMatrix::~DistSCMatrix()
{
}

double
DistSCMatrix::get_element(int i,int j) const
{
  double res;
  double *e = find_element(i,j);
  if (e) {
      res = *e;
      messagegrp()->bcast(res, messagegrp()->me());
    }
  else {
      messagegrp()->bcast(res, element_to_node(i, j));
    }
  return res;
}

void
DistSCMatrix::set_element(int i,int j,double a)
{
  double *e = find_element(i,j);
  if (e) {
      *e = a;
    }
}

void
DistSCMatrix::accumulate_element(int i,int j,double a)
{
  double *e = find_element(i,j);
  if (e) {
      *e += a;
    }
}

SCMatrix *
DistSCMatrix::get_subblock(int br, int er, int bc, int ec)
{
  error("get_subblock");
  return 0;
}

void
DistSCMatrix::assign_subblock(SCMatrix*sb, int br, int er, int bc, int ec,
                              int source_br, int source_bc)
{
  error("assign_subblock");
}

void
DistSCMatrix::accumulate_subblock(SCMatrix*sb, int br, int er, int bc, int ec,
                                  int source_br, int source_bc)
{
  error("accumulate_subblock");
}

SCVector *
DistSCMatrix::get_row(int i)
{
  error("get_row");
  return 0;
}

void
DistSCMatrix::assign_row(SCVector *v, int i)
{
  error("assign_row");
}

void
DistSCMatrix::accumulate_row(SCVector *v, int i)
{
  error("accumulate_row");
}

SCVector *
DistSCMatrix::get_column(int i)
{
  double *col = new double[nrow()];
  
  Ref<SCMatrixSubblockIter> iter = local_blocks(SCMatrixSubblockIter::Read);
  for (iter->begin(); iter->ready(); iter->next()) {
    SCMatrixRectBlock *b = dynamic_cast<SCMatrixRectBlock*>(iter->block());
    if (b->jstart > i || b->jend <= i)
      continue;

    int joff = i-b->jstart;
    int jlen = b->jend-b->jstart;
    int ist = 0;
    
    for (int ii=b->istart; ii < b->iend; ii++,ist++)
      col[ii] = b->data[ist*jlen+joff];
  }
  
  SCVector * rcol = kit_->vector(rowdim());
  rcol->assign(col);

  delete[] col;

  return rcol;
}

void
DistSCMatrix::assign_column(SCVector *v, int i)
{
  error("assign_column");
}

void
DistSCMatrix::accumulate_column(SCVector *v, int i)
{
  error("accumulate_column");
}

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

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

  SCMatrixBlockListIter i1, i2;
  for (i1=la->blocklist->begin(),i2=blocklist->begin();
       i1!=la->blocklist->end() && i2!=blocklist->end();
       i1++,i2++) {
      int n = i1.block()->ndat();
      if (n != i2.block()->ndat()) {
          ExEnv::errn() << indent
               << "DistSCMatrixListSubblockIter::accumulate block mismatch: "
               << "internal error" << endl;
          abort();
        }
      double *dat1 = i1.block()->dat();
      double *dat2 = i2.block()->dat();
      for (int i=0; i<n; i++) {
          dat2[i] += dat1[i];
        }
    }
}

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

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

  Ref<SCMatrixSubblockIter> I = ((SymmSCMatrix*)a)->all_blocks(SCMatrixSubblockIter::Read);
  for (I->begin(); I->ready(); I->next()) {
      Ref<SCMatrixBlock> block = I->block();
      if (DEBUG)
          ExEnv::outn() << messagegrp()->me() << ": "
               << block->class_name()
               << "(" << block->blocki << ", " << block->blockj << ")"
               << endl;
      // see if i've got this block
      Ref<SCMatrixBlock> localblock
          = block_to_block(block->blocki,block->blockj);
      if (localblock.nonnull()) {
          // the diagonal blocks require special handling
          if (block->blocki == block->blockj) {
              int n = rowblocks()->size(block->blocki);
              double *dat1 = block->dat();
              double *dat2 = localblock->dat();
              for (int i=0; i<n; i++) {
                  for (int j=0; j<i; j++) {
                      double tmp = *dat1;
                      dat2[i*n+j] += tmp;
                      dat2[j*n+i] += tmp;
                      dat1++;
                    }
                  dat2[i*n+i] = *dat1++;
                }
            }
          else {
              int n = block->ndat();
              double *dat1 = block->dat();
              double *dat2 = localblock->dat();
              for (int i=0; i<n; i++) {
                  dat2[i] += dat1[i];
                }
            }
        }
      // now for the transpose
      if (block->blocki != block->blockj) {
          localblock = block_to_block(block->blockj,block->blocki);
          if (localblock.nonnull()) {
              int nr = rowblocks()->size(block->blocki);
              int nc = rowblocks()->size(block->blockj);
              double *dat1 = block->dat();
              double *dat2 = localblock->dat();
              for (int i=0; i<nr; i++) {
                  for (int j=0; j<nc; j++) {
                      dat2[j*nr+i] += *dat1++;
                    }
                }
            }
        }
    }
}

void
DistSCMatrix::accumulate(const DiagSCMatrix*a)
{
  // make sure that the argument is of the correct type
  const DistDiagSCMatrix* la
    = require_dynamic_cast<const DistDiagSCMatrix*>(a,"DistSCMatrix::accumulate");

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

  Ref<SCMatrixSubblockIter> I = ((DiagSCMatrix*)a)->all_blocks(SCMatrixSubblockIter::Read);
  for (I->begin(); I->ready(); I->next()) {
      Ref<SCMatrixBlock> block = I->block();
      // see if i've got this block
      Ref<SCMatrixBlock> localblock
          = block_to_block(block->blocki,block->blockj);
      if (localblock.nonnull()) {
          int n = rowblocks()->size(block->blocki);
          double *dat1 = block->dat();
          double *dat2 = localblock->dat();
          for (int i=0; i<n; i++) {
              dat2[i*n+i] += *dat1++;
            }
        }
    }
}

void
DistSCMatrix::accumulate(const SCVector*a)
{
  // make sure that the argument is of the correct type
  const DistSCVector* la
    = require_dynamic_cast<const DistSCVector*>(a,"DistSCMatrix::accumulate");

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

  SCMatrixBlockListIter I, J;
  for (I = blocklist->begin(), J = la->blocklist->begin();
       I != blocklist->end() && J != la->blocklist->end();
       I++,J++) {
      int n = I.block()->ndat();
      if (n != J.block()->ndat()) {
          ExEnv::errn() << indent << "DistSCMatrix::accumulate(SCVector*a): "
               << "block lists do not match" << endl;
          abort();
        }
      double *dati = I.block()->dat();
      double *datj = J.block()->dat();
      for (int i=0; i<n; i++) dati[i] += datj[i];
    }
}

void
DistSCMatrix::accumulate_product_rr(SCMatrix*pa,SCMatrix*pb)
{
  const char* name = "DistSCMatrix::accumulate_product_rr";
  // make sure that the arguments are of the correct type
  DistSCMatrix* a = require_dynamic_cast<DistSCMatrix*>(pa,name);
  DistSCMatrix* b = require_dynamic_cast<DistSCMatrix*>(pb,name);

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

  // i need this in row form and a in row form
  create_vecform(Row);
  vecform_zero();
  a->create_vecform(Row);
  a->vecform_op(CopyToVec);

  Ref<SCMatrixSubblockIter> I = b->all_blocks(SCMatrixSubblockIter::Read);
  for (I->begin(); I->ready(); I->next()) {
      Ref<SCMatrixRectBlock> blk
          = dynamic_cast<SCMatrixRectBlock*>(I->block());
      int kk,k,jj,j,i,nj;
      nj = blk->jend - blk->jstart;
      double *data = blk->data;
      for (i=0; i<nvec; i++) {
          double *veci = vec[i];
          double *aveci = a->vec[i];
          for (j=blk->jstart,jj=0; j<blk->jend; j++,jj++) {
              double tmp = 0.0;
              for (k=blk->istart,kk=0; k<blk->iend; k++,kk++) {
                  tmp += aveci[k] * data[kk*nj+jj];
                }
              veci[j] += tmp;
            }
        }
    }

  vecform_op(AccumFromVec);