blockedrect.cc

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

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

  // make sure that the dimensions match
  if (!rowdim()->equiv(la->dim()) || !coldim()->equiv(la->dim())) {
    ExEnv::errn() << indent << "BlockedSCMatrix::accumulate(DiagSCMatrix*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 SCVector*a)
{
  // make sure that the arguments is of the correct type
  const BlockedSCVector* la
    = require_dynamic_cast<const BlockedSCVector*>(a,"BlockedSCVector::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 << "BlockedSCMatrix::accumulate(SCVector*a): "
         << "dimensions don't match\n";
    abort();
  }

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

void
BlockedSCMatrix::transpose_this()
{
  for (int i=0; i < nblocks_; i++)
    if (mats_[i].nonnull())
      mats_[i]->transpose_this();
  
  RefSCDimension tmp = d1;
  d1 = d2;
  d2 = tmp;
}

// hack, hack, hack!  One day we'll get svd working everywhere.
double
BlockedSCMatrix::invert_this()
{
  int i;
  double res=1;

  // if this matrix is block diagonal, then give a normal inversion a shot
  if (d1->blocks()->nblock() == d2->blocks()->nblock()) {
    for (i=0; i < nblocks_; i++)
      if (mats_[i].nonnull()) res *= mats_[i]->invert_this();
    return res;
  }

  // ok, let's make sure that the matrix is at least square
  if (d1->n() != d2->n()) {
    ExEnv::errn() << indent
         << "BlockedSCMatrix::invert_this: SVD not implemented yet\n";
    abort();
  }

  if (d1->blocks()->nblock() == 1) {
    RefSCMatrix tdim = subkit->matrix(d1->blocks()->subdim(0),
                                      d1->blocks()->subdim(0));
    tdim->convert(this);
    res = tdim->invert_this();
    transpose_this();

    // d1 and d2 were swapped by now
    for (i=0; i < d1->blocks()->nblock(); i++)
      if (mats_[i].nonnull())
        mats_[i]->convert(tdim.get_subblock(d1->blocks()->start(i),
                                            d1->blocks()->fence(i)-1,
                                            0, d2->n()-1));
    
    return res;

  } else if (d2->blocks()->nblock() == 1) {
    RefSCMatrix tdim = subkit->matrix(d2->blocks()->subdim(0),
                                      d2->blocks()->subdim(0));

    tdim->convert(this);
    res = tdim->invert_this();
    transpose_this();

    // d1 and d2 were swapped by now
    for (i=0; i < d2->blocks()->nblock(); i++)
      if (mats_[i].nonnull())
        mats_[i]->convert(tdim.get_subblock(0, d1->n()-1,
                                            d2->blocks()->start(i),
                                            d2->blocks()->fence(i)-1));
    
    return res;

  } else {
    ExEnv::errn() << indent
         << "BlockedSCMatrix::invert_this: SVD not implemented yet\n";
    abort();
  }

  return 0.0;
}

void
BlockedSCMatrix::gen_invert_this()
{
  ExEnv::errn() << indent
       << "BlockedSCMatrix::gen_invert_this: SVD not implemented yet\n";
  abort();
}

double
BlockedSCMatrix::determ_this()
{
  double res=1;
  
  for (int i=0; i < nblocks_; i++)
    if (mats_[i].nonnull())
      res *= mats_[i]->determ_this();

  return res;
}

double
BlockedSCMatrix::trace()
{
  double ret=0;
  for (int i=0; i < nblocks_; i++)
    if (mats_[i].nonnull())
      ret += mats_[i]->trace();
  
  return ret;
}

void
BlockedSCMatrix::svd_this(SCMatrix *U, DiagSCMatrix *sigma, SCMatrix *V)
{
  BlockedSCMatrix* lU =
    require_dynamic_cast<BlockedSCMatrix*>(U,"BlockedSCMatrix::svd_this");
  BlockedSCMatrix* lV =
    require_dynamic_cast<BlockedSCMatrix*>(V,"BlockedSCMatrix::svd_this");
  BlockedDiagSCMatrix* lsigma =
    require_dynamic_cast<BlockedDiagSCMatrix*>(sigma,"BlockedSCMatrix::svd_this");

  for (int i=0; i < nblocks_; i++)
    if (mats_[i].nonnull())
      mats_[i]->svd_this(lU->mats_[i], lsigma->mats_[i], lV->mats_[i]);
}

double
BlockedSCMatrix::solve_this(SCVector*v)
{
  double res=1;
  
  BlockedSCVector* lv =
    require_dynamic_cast<BlockedSCVector*>(v,"BlockedSCMatrix::solve_this");
  
  // make sure that the dimensions match
  if (!rowdim()->equiv(lv->dim())) {
    ExEnv::errn() << indent << "BlockedSCMatrix::solve_this(SCVector*v): "
         << "dimensions don't match\n";
    abort();
  }

  for (int i=0; i < nblocks_; i++)
    if (mats_[i].nonnull())
      res *= mats_[i]->solve_this(lv->vecs_[i]);

  return res;
}

void
BlockedSCMatrix::schmidt_orthog(SymmSCMatrix *S, int nc)
{
  BlockedSymmSCMatrix* lS =
    require_dynamic_cast<BlockedSymmSCMatrix*>(S,"BlockedSCMatrix::schmidt_orthog");
  
  // make sure that the dimensions match
  if (!rowdim()->equiv(lS->dim())) {
    ExEnv::errn() << indent << "BlockedSCMatrix::schmidt_orthog(): "
         << "dimensions don't match\n";
    abort();
  }

  for (int i=0; i < nblocks_; i++)
    if (mats_[i].nonnull())
      mats_[i]->schmidt_orthog(lS->mats_[i].pointer(),
                               lS->dim()->blocks()->subdim(i).n());
}

int
BlockedSCMatrix::schmidt_orthog_tol(SymmSCMatrix *S, double tol, double *res)
{
  ExEnv::err0()
               << "ERROR: BlockedSCMatrix::schmidt_orthog_tol doesn't exist"
               << endl;
  abort();
  return 0;
}

void
BlockedSCMatrix::element_op(const Ref<SCElementOp>& op)
{
  BlockedSCElementOp *bop = dynamic_cast<BlockedSCElementOp*>(op.pointer());

  op->defer_collect(1);
  for (int i=0; i < nblocks_; i++) {
    if (bop)
      bop->working_on(i);
    if (mats_[i].nonnull())
      mats_[i]->element_op(op);
  }
  op->defer_collect(0);
  if (op->has_collect()) op->collect(messagegrp());
}

void
BlockedSCMatrix::element_op(const Ref<SCElementOp2>& op,
                          SCMatrix* m)
{
  BlockedSCMatrix *lm
    = require_dynamic_cast<BlockedSCMatrix*>(m,"BlockedSCMatrix::element_op");

  if (!rowdim()->equiv(lm->rowdim()) || !coldim()->equiv(lm->coldim())) {
    ExEnv::errn() << indent << "BlockedSCMatrix: bad element_op\n";
    abort();
  }

  BlockedSCElementOp2 *bop = dynamic_cast<BlockedSCElementOp2*>(op.pointer());

  op->defer_collect(1);
  for (int i=0; i < nblocks_; i++) {
    if (bop)
      bop->working_on(i);
    if (mats_[i].nonnull())
      mats_[i]->element_op(op,lm->mats_[i].pointer());
  }
  op->defer_collect(0);
  if (op->has_collect()) op->collect(messagegrp());
}

void
BlockedSCMatrix::element_op(const Ref<SCElementOp3>& op,
                          SCMatrix* m,SCMatrix* n)
{
  BlockedSCMatrix *lm
    = require_dynamic_cast<BlockedSCMatrix*>(m,"BlockedSCMatrix::element_op");
  BlockedSCMatrix *ln
    = require_dynamic_cast<BlockedSCMatrix*>(n,"BlockedSCMatrix::element_op");

  if (!rowdim()->equiv(lm->rowdim()) || !coldim()->equiv(lm->coldim()) ||
      !rowdim()->equiv(ln->rowdim()) || !coldim()->equiv(ln->coldim())) {
    ExEnv::errn() << indent << "BlockedSCMatrix: bad element_op\n";
    abort();
  }

  BlockedSCElementOp3 *bop = dynamic_cast<BlockedSCElementOp3*>(op.pointer());

  op->defer_collect(1);
  for (int i=0; i < nblocks_; i++) {
    if (bop)
      bop->working_on(i);
    if (mats_[i].nonnull())
      mats_[i]->element_op(op,lm->mats_[i].pointer(),
                              ln->mats_[i].pointer());
  }
  op->defer_collect(0);
  if (op->has_collect()) op->collect(messagegrp());
}

void
BlockedSCMatrix::vprint(const char *title, ostream& os, int prec) const
{
  int len = (title) ? strlen(title) : 0;
  char *newtitle = new char[len + 80];

  for (int i=0; i < nblocks_; i++) {
    if (mats_[i].null())
      continue;

    sprintf(newtitle,"%s:  block %d",title,i+1);
    mats_[i]->print(newtitle, os, prec);
  }

  delete[] newtitle;
}

RefSCDimension
BlockedSCMatrix::rowdim(int i) const
{
  return d1->blocks()->subdim(i);
}

RefSCDimension
BlockedSCMatrix::coldim(int i) const
{
  return d2->blocks()->subdim(i);
}

int
BlockedSCMatrix::nblocks() const
{
  return nblocks_;
}

RefSCMatrix
BlockedSCMatrix::block(int i)
{
  if (mats_)
      return mats_[i];
  else
      return (SCMatrix*)0;
}

Ref<SCMatrixSubblockIter>
BlockedSCMatrix::local_blocks(SCMatrixSubblockIter::Access access)
{
  Ref<SCMatrixCompositeSubblockIter> iter
      = new SCMatrixCompositeSubblockIter(access, nblocks());
  for (int i=0; i<nblocks(); i++) {
      if (block(i).null())
          iter->set_iter(i, new SCMatrixNullSubblockIter(access));
      else
          iter->set_iter(i, block(i)->local_blocks(access));
    }
  Ref<SCMatrixSubblockIter> ret = iter.pointer();
  return ret;
}

Ref<SCMatrixSubblockIter>
BlockedSCMatrix::all_blocks(SCMatrixSubblockIter::Access access)
{
  Ref<SCMatrixCompositeSubblockIter> iter
      = new SCMatrixCompositeSubblockIter(access, nblocks());
  for (int i=0; i<nblocks(); i++) {
      if (block(i).null())
          iter->set_iter(i, new SCMatrixNullSubblockIter(access));
      else
          iter->set_iter(i, block(i)->all_blocks(access));
    }
  Ref<SCMatrixSubblockIter> ret = iter.pointer();
  return ret;
}

void
BlockedSCMatrix::save(StateOut&s)
{
  int nr = nrow();
  int nc = ncol();
  s.put(nr);
  s.put(nc);
  int has_subblocks = 1;
  s.put(has_subblocks);
  s.put(nblocks());
  for (int i=0; i<nblocks(); i++) {
      block(i).save(s);
    }
}

void
BlockedSCMatrix::restore(StateIn& s)
{
  int nrt, nr = nrow();
  int nct, nc = ncol();
  s.get(nrt);
  s.get(nct);
  if (nrt != nr || nct != nc) {
      ExEnv::errn() << indent << "BlockedSCMatrix::restore(): bad dimensions" << endl;
      abort();
    }
  int has_subblocks;
  s.get(has_subblocks);
  if (has_subblocks) {
      int nblock;
      s.get(nblock);
      if (nblock != nblocks()) {
          ExEnv::errn() << indent
               << "BlockedSCMatrix::restore(): nblock differs\n" << endl;
          abort();
        }
      for (int i=0; i<nblocks(); i++) {
          block(i).restore(s);
        }
    }
  else {
      ExEnv::errn() << indent
           << "BlockedSCMatrix::restore(): no subblocks--cannot restore"
           << endl;
      abort();
    }
}

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

// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End: