abstract.cc

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

//
// abstract.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 <math.h>

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

using namespace std;
using namespace sc;

/////////////////////////////////////////////////////////////////////////
// These member are used by the abstract SCMatrix classes.
/////////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////////
// SCMatrixKit members

static ClassDesc SCMatrixKit_cd(
  typeid(SCMatrixKit),"SCMatrixKit",1,"public DescribedClass",
  0, 0, 0);

SCMatrixKit::SCMatrixKit()
{
  grp_ = MessageGrp::get_default_messagegrp();
}

SCMatrixKit::SCMatrixKit(const Ref<KeyVal>& keyval)
{
  grp_ << keyval->describedclassvalue("messagegrp");
  if (grp_.null()) grp_ = MessageGrp::get_default_messagegrp();
}

SCMatrixKit::~SCMatrixKit()
{
}

SCMatrix*
SCMatrixKit::restore_matrix(StateIn& s,
                            const RefSCDimension& d1,
                            const RefSCDimension& d2)
{
  SCMatrix *r = matrix(d1,d2);
  r->restore(s);
  return r;
}

SymmSCMatrix*
SCMatrixKit::restore_symmmatrix(StateIn& s, const RefSCDimension& d)
{
  SymmSCMatrix *r = symmmatrix(d);
  r->restore(s);
  return r;
}

DiagSCMatrix*
SCMatrixKit::restore_diagmatrix(StateIn& s, const RefSCDimension& d)
{
  DiagSCMatrix *r = diagmatrix(d);
  r->restore(s);
  return r;
}

SCVector*
SCMatrixKit::restore_vector(StateIn& s, const RefSCDimension& d)
{
  SCVector *r = vector(d);
  r->restore(s);
  return r;
}

Ref<MessageGrp>
SCMatrixKit::messagegrp() const
{
  return grp_;
}

/////////////////////////////////////////////////////////////////////////
// SCMatrix members

static ClassDesc SCMatrix_cd(
  typeid(SCMatrix),"SCMatrix",1,"public DescribedClass",
  0, 0, 0);

SCMatrix::SCMatrix(const RefSCDimension&a1, const RefSCDimension&a2,
                   SCMatrixKit*kit):
  d1(a1),
  d2(a2),
  kit_(kit)
{
}

SCMatrix::~SCMatrix()
{
}

void
SCMatrix::save(StateOut&s)
{
  int nr = nrow();
  int nc = ncol();
  s.put(nr);
  s.put(nc);
  int has_subblocks = 0;
  s.put(has_subblocks);
  for (int i=0; i<nr; i++) {
      for (int j=0; j<nc; j++) {
          s.put(get_element(i,j));
        }
    }
}

void
SCMatrix::restore(StateIn& s)
{
  int nrt, nr = nrow();
  int nct, nc = ncol();
  s.get(nrt);
  s.get(nct);
  if (nrt != nr || nct != nc) {
      ExEnv::errn() << "SCMatrix::restore(): bad dimensions" << endl;
      abort();
    }
  int has_subblocks;
  s.get(has_subblocks);
  if (!has_subblocks) {
      for (int i=0; i<nr; i++) {
          for (int j=0; j<nc; j++) {
              double tmp;
              s.get(tmp);
              set_element(i,j, tmp);
            }
        }
    }
  else {
      ExEnv::errn() << "SCMatrix::restore(): matrix has subblocks--cannot restore"
           << endl;
      abort();
    }
}

double
SCMatrix::maxabs() const
{
  Ref<SCElementMaxAbs> op = new SCElementMaxAbs();
  Ref<SCElementOp> abop = op.pointer();
  ((SCMatrix *)this)->element_op(abop);
  return op->result();
}

void
SCMatrix::randomize()
{
  Ref<SCElementOp> op = new SCElementRandomize();
  this->element_op(op);
}

void
SCMatrix::assign_val(double a)
{
  Ref<SCElementOp> op = new SCElementAssign(a);
  this->element_op(op);
}

void
SCMatrix::scale(double a)
{
  Ref<SCElementOp> op = new SCElementScale(a);
  this->element_op(op);
}

void
SCMatrix::scale_diagonal(double a)
{
  Ref<SCElementOp> op = new SCElementScaleDiagonal(a);
  this->element_op(op);
}

void
SCMatrix::shift_diagonal(double a)
{
  Ref<SCElementOp> op = new SCElementShiftDiagonal(a);
  this->element_op(op);
}

void
SCMatrix::unit()
{
  this->assign(0.0);
  this->shift_diagonal(1.0);
}

void
SCMatrix::assign_r(SCMatrix*a)
{
  this->assign(0.0);
  this->accumulate(a);
}

void
SCMatrix::assign_p(const double*a)
{
  int i;
  int nr = nrow();
  int nc = ncol();
  // some compilers need the following cast
  const double **v = (const double**) new double*[nr];
  for (i=0; i<nr; i++) {
      v[i] = &a[i*nc];
    }
  assign(v);
  delete[] v;
}

void
SCMatrix::assign_pp(const double**a)
{
  int i;
  int j;
  int nr;
  int nc;
  nr = nrow();
  nc = ncol();
  for (i=0; i<nr; i++) {
      for (j=0; j<nc; j++) {
          set_element(i,j,a[i][j]);
        }
    }
}

void
SCMatrix::convert(SCMatrix*a)
{
  assign(0.0);
  convert_accumulate(a);
}

void
SCMatrix::convert_accumulate(SCMatrix*a)
{
  Ref<SCElementOp> op = new SCElementAccumulateSCMatrix(a);
  element_op(op);
}

void
SCMatrix::convert(double*a) const
{
  int i;
  int nr = nrow();
  int nc = ncol();
  double **v = new double*[nr];
  for (i=0; i<nr; i++) {
      v[i] = &a[i*nc];
    }
  convert(v);
  delete[] v;
}

void
SCMatrix::convert(double**a) const
{
  int i, j;
  int nr, nc;
  nr = nrow();
  nc = ncol();
  for (i=0; i<nr; i++) {
      for (j=0; j<nc; j++) {
          a[i][j] = get_element(i,j);
        }
    }
}

void
SCMatrix::accumulate_product_sr(SymmSCMatrix*a,SCMatrix*b)
{
  SCMatrix *t = b->copy();
  t->transpose_this();
  SCMatrix *t2 = this->copy();
  t2->transpose_this();
  t2->accumulate_product(t,a);
  delete t;
  t2->transpose_this();
  assign(t2);
  delete t2;
}

void
SCMatrix::accumulate_product_dr(DiagSCMatrix*a,SCMatrix*b)
{
  SCMatrix *t = b->copy();
  t->transpose_this();
  SCMatrix *t2 = kit()->matrix(coldim(),rowdim());
  t2->assign(0.0);
  t2->accumulate_product(t,a);
  delete t;
  t2->transpose_this();
  accumulate(t2);
  delete t2;
}

void
SCMatrix::print(ostream&o) const
{
  vprint(0, o, 10);
}

void
SCMatrix::print(const char *t, ostream&o, int i) const
{
  vprint(t, o, i);
}

SCMatrix*
SCMatrix::clone()
{
  return kit()->matrix(rowdim(),coldim());
}

SCMatrix*
SCMatrix::copy()
{
  SCMatrix* result = clone();
  result->assign(this);
  return result;
}

void
SCMatrix::gen_invert_this()
{
  int i;

  // Compute the singular value decomposition of this
  RefSCMatrix U(rowdim(),rowdim(),kit());
  RefSCMatrix V(coldim(),coldim(),kit());
  RefSCDimension min;
  if (coldim().n()<rowdim().n()) min = coldim();
  else min = rowdim();
  int nmin = min.n();
  RefDiagSCMatrix sigma(min,kit());
  svd_this(U.pointer(),sigma.pointer(),V.pointer());

  // compute the epsilon rank of this
  int rank = 0;
  for (i=0; i<nmin; i++) {
      if (fabs(sigma(i)) > 0.0000001) rank++;
    }

  RefSCDimension drank = new SCDimension(rank);
  RefDiagSCMatrix sigma_i(drank,kit());
  for (i=0; i<rank; i++) {
      sigma_i(i) = 1.0/sigma(i);
    }
  RefSCMatrix Ur(rowdim(), drank, kit());
  RefSCMatrix Vr(coldim(), drank, kit());
  Ur.assign_subblock(U,0, rowdim().n()-1, 0, drank.n()-1, 0, 0);
  Vr.assign_subblock(V,0, coldim().n()-1, 0, drank.n()-1, 0, 0);
  assign((Vr * sigma_i * Ur.t()).t());
  transpose_this();
}

void
SCMatrix::svd_this(SCMatrix *U, DiagSCMatrix *sigma, SCMatrix *V)
{
  ExEnv::errn() << indent << class_name() << ": SVD not implemented\n";
  abort();
}

void
SCMatrix::accumulate_product_rs(SCMatrix*a,SymmSCMatrix*b)
{
  SCMatrix *brect = kit()->matrix(b->dim(),b->dim());
  brect->assign(0.0);
  brect->accumulate(b);
  accumulate_product(a,brect);
  delete brect;
}

void
SCMatrix::accumulate_product_ss(SymmSCMatrix*a,SymmSCMatrix*b)
{
  SCMatrix *arect = kit()->matrix(a->dim(),a->dim());
  arect->assign(0.0);
  arect->accumulate(a);
  SCMatrix *brect = kit()->matrix(b->dim(),b->dim());
  brect->assign(0.0);
  brect->accumulate(b);
  accumulate_product(arect,brect);
  delete arect;
  delete brect;
}

void
SCMatrix::accumulate_product_rd(SCMatrix*a,DiagSCMatrix*b)
{
  SCMatrix *brect = kit()->matrix(b->dim(),b->dim());
  brect->assign(0.0);
  brect->accumulate(b);
  accumulate_product(a,brect);
  delete brect;
}

Ref<MessageGrp>
SCMatrix::messagegrp() const
{
  return kit_->messagegrp();
}

/////////////////////////////////////////////////////////////////////////
// SymmSCMatrix member functions

static ClassDesc SymmSCMatrix_cd(
  typeid(SymmSCMatrix),"SymmSCMatrix",1,"public DescribedClass",
  0, 0, 0);

SymmSCMatrix::SymmSCMatrix(const RefSCDimension&a, SCMatrixKit *kit):
  d(a),
  kit_(kit)
{
}

void
SymmSCMatrix::save(StateOut&s)
{
  int nr = n();
  s.put(nr);
  for (int i=0; i<nr; i++) {
      for (int j=0; j<=i; j++) {
          s.put(get_element(i,j));
        }
    }
}

void
SymmSCMatrix::restore(StateIn& s)
{
  int nrt, nr = n();
  s.get(nrt);
  if (nrt != nr) {
      ExEnv::errn() << "SymmSCMatrix::restore(): bad dimension" << endl;
      abort();
    }
  for (int i=0; i<nr; i++) {
      for (int j=0; j<=i; j++) {
          double tmp;
          s.get(tmp);
          set_element(i,j, tmp);
        }
    }
}

double
SymmSCMatrix::maxabs() const
{
  Ref<SCElementMaxAbs> op = new SCElementMaxAbs();
  Ref<SCElementOp> abop = op.pointer();
  ((SymmSCMatrix*)this)->element_op(abop);
  return op->result();
}

void
SymmSCMatrix::randomize()
{
  Ref<SCElementOp> op = new SCElementRandomize();
  this->element_op(op);
}

void
SymmSCMatrix::assign_val(double a)
{
  Ref<SCElementOp> op = new SCElementAssign(a);
  this->element_op(op);
}

void
SymmSCMatrix::assign_p(const double*a)
{
  int i;
  int nr = n();
  // some compilers need the following cast
  const double **v = (const double **) new double*[nr];
  int ioff= 0;
  for (i=0; i<nr; i++) {
      ioff += i;
      v[i] = &a[ioff];
//    ioff += i;
    }
  assign(v);
  delete[] v;
}

void
SymmSCMatrix::assign_pp(const double**a)
{
  int i;
  int j;
  int nr;
  nr = n();
  for (i=0; i<nr; i++) {
      for (j=0; j<=i; j++) {
          set_element(i,j,a[i][j]);
        }
    }
}

void
SymmSCMatrix::convert(SymmSCMatrix*a)
{
  assign(0.0);
  convert_accumulate(a);
}

void
SymmSCMatrix::convert_accumulate(SymmSCMatrix*a)
{
  Ref<SCElementOp> op = new SCElementAccumulateSymmSCMatrix(a);
  element_op(op);
}

void
SymmSCMatrix::convert(double*a) const
{
  int i;
  int nr = n();
  double **v = new double*[nr];
  int ioff= 0;
  for (i=0; i<nr; i++) {
      ioff += i;
      v[i] = &a[ioff];