matrix.cc

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

RefSCMatrix
RefSCMatrix::block(int i) const
{
  BlockedSCMatrix *b = dynamic_cast<BlockedSCMatrix*>(pointer());
  if (b) return b->block(i);
  return *this;
}

///////////////////////////////////////////////////////////////////
// RefSymmSCMatrix members

RefSymmSCMatrix::RefSymmSCMatrix()
{
}
             
RefSymmSCMatrix::RefSymmSCMatrix (const RefSymmSCMatrix & o):
  Ref<SymmSCMatrix> (o)
{
}
             
RefSymmSCMatrix::RefSymmSCMatrix (SymmSCMatrix * o):
  Ref<SymmSCMatrix> (o)
{
}

RefSymmSCMatrix::~RefSymmSCMatrix ()
{
}

RefSymmSCMatrix&
RefSymmSCMatrix::operator=(SymmSCMatrix* cr)
{
  Ref<SymmSCMatrix>::operator=(cr);
  return *this;
}

RefSymmSCMatrix&
RefSymmSCMatrix::operator=(const RefSymmSCMatrix & c)
{
  Ref<SymmSCMatrix>::operator=(c);
  return *this;
}

RefSymmSCMatrix::RefSymmSCMatrix(const RefSCDimension&a,
                                 const Ref<SCMatrixKit>&k)
{
  assign_pointer(k->symmmatrix(a));
}

void
RefSymmSCMatrix::set_element(int i, int j, double a) const
{
  require_nonnull();
  pointer()->set_element(i,j,a);
}

void
RefSymmSCMatrix::accumulate_element(int i, int j, double a) const
{
  require_nonnull();
  pointer()->accumulate_element(i,j,a);
}

double
RefSymmSCMatrix::get_element(int i, int j) const
{
  require_nonnull();
  return pointer()->get_element(i,j);
}

RefSCMatrix
RefSymmSCMatrix::get_subblock(int br, int er, int bc, int ec)
{
  require_nonnull();
  
  RefSCMatrix ret = pointer()->get_subblock(br,er,bc,ec);
  return ret;
}

RefSymmSCMatrix
RefSymmSCMatrix::get_subblock(int br, int er)
{
  require_nonnull();
  
  RefSymmSCMatrix ret = pointer()->get_subblock(br,er);
  return ret;
}

void
RefSymmSCMatrix::assign_subblock(const RefSCMatrix& sb,
                                 int br, int er, int bc, int ec)
{
  require_nonnull();
  sb.require_nonnull();
  pointer()->assign_subblock(sb.pointer(),br,er,bc,ec);
}

void
RefSymmSCMatrix::assign_subblock(const RefSymmSCMatrix& sb, int br, int er)
{
  require_nonnull();
  sb.require_nonnull();
  pointer()->assign_subblock(sb.pointer(),br,er);
}

void
RefSymmSCMatrix::accumulate_subblock(const RefSCMatrix& sb,
                                     int br, int er, int bc, int ec)
{
  require_nonnull();
  sb.require_nonnull();
  pointer()->accumulate_subblock(sb.pointer(),br,er,bc,ec);
}

void
RefSymmSCMatrix::accumulate_subblock(const RefSymmSCMatrix& sb, int br, int er)
{
  require_nonnull();
  sb.require_nonnull();
  pointer()->accumulate_subblock(sb.pointer(),br,er);
}

RefSCVector
RefSymmSCMatrix::get_row(int i)
{
  require_nonnull();
  
  RefSCVector ret = pointer()->get_row(i);
  return ret;
}

void
RefSymmSCMatrix::assign_row(const RefSCVector& v, int i)
{
  require_nonnull();
  v.require_nonnull();
  pointer()->assign_row(v.pointer(),i);
}

void
RefSymmSCMatrix::accumulate_row(const RefSCVector& v, int i)
{
  require_nonnull();
  v.require_nonnull();
  pointer()->accumulate_row(v.pointer(),i);
}

void
RefSymmSCMatrix::accumulate_symmetric_product(const RefSCMatrix& a) const
{
  require_nonnull();
  pointer()->accumulate_symmetric_product(a.pointer());
}

void
RefSymmSCMatrix::accumulate_symmetric_sum(const RefSCMatrix& a) const
{
  require_nonnull();
  pointer()->accumulate_symmetric_sum(a.pointer());
}

void
RefSymmSCMatrix::accumulate_transform(const RefSCMatrix& a,
                                      const RefSymmSCMatrix&b,
                                      SCMatrix::Transform t) const
{
  require_nonnull();
  pointer()->accumulate_transform(a.pointer(),b.pointer(),t);
}

void
RefSymmSCMatrix::accumulate_transform(const RefSCMatrix& a,
                                      const RefDiagSCMatrix&b,
                                      SCMatrix::Transform t) const
{
  require_nonnull();
  pointer()->accumulate_transform(a.pointer(),b.pointer(),t);
}

void
RefSymmSCMatrix::accumulate_transform(const RefSymmSCMatrix& a,
                                      const RefSymmSCMatrix&b) const
{
  require_nonnull();
  pointer()->accumulate_transform(a.pointer(),b.pointer());
}

RefSymmSCMatrix
RefSymmSCMatrix::operator+(const RefSymmSCMatrix&a) const
{
  require_nonnull();
  a.require_nonnull();

  RefSymmSCMatrix ret(dim(),kit());
  
  ret->assign(pointer());
  ret->accumulate(a.pointer());

  return ret;
}

RefSymmSCMatrix
RefSymmSCMatrix::operator-(const RefSymmSCMatrix&a) const
{
  require_nonnull();
  a.require_nonnull();

  RefSymmSCMatrix ret(dim(),kit());
  
  ret->assign(a.pointer());
  ret->scale(-1.0);
  ret->accumulate(pointer());

  return ret;
}

RefSymmSCMatrix
RefSymmSCMatrix::i() const
{
  require_nonnull();
  
  RefSymmSCMatrix ret;
  ret = clone();
  ret->assign(pointer());
  ret->invert_this();
  return ret;
}

RefSymmSCMatrix
RefSymmSCMatrix::gi() const
{
  require_nonnull();
  
  RefSymmSCMatrix ret;
  ret = clone();
  ret->assign(pointer());
  ret->gen_invert_this();
  return ret;
}

int
RefSymmSCMatrix::n() const
{
  if (null()) return 0;
  else return pointer()->dim()->n();
}

RefSCDimension
RefSymmSCMatrix::dim() const
{
  if (null()) return 0;
  else return pointer()->dim();
}

Ref<SCMatrixKit>
RefSymmSCMatrix::kit() const
{
  if (null()) return 0;
  else return pointer()->kit();
}

SymmSCMatrixdouble
RefSymmSCMatrix::operator()(int i,int j) const
{
  return SymmSCMatrixdouble(pointer(),i,j);
}

RefSymmSCMatrix
RefSymmSCMatrix::clone() const
{
  RefSymmSCMatrix r = kit()->symmmatrix(dim());
  return r;
}

RefDiagSCMatrix
RefSymmSCMatrix::eigvals() const
{
  if (null()) return 0;
  RefDiagSCMatrix vals = kit()->diagmatrix(dim());
  RefSCMatrix vecs = kit()->matrix(dim(),dim());
  diagonalize(vals,vecs);
  return vals;
}

RefSCMatrix
RefSymmSCMatrix::eigvecs() const
{
  if (null()) return 0;
  RefDiagSCMatrix vals = kit()->diagmatrix(dim());
  RefSCMatrix vecs = kit()->matrix(dim(),dim());
  diagonalize(vals,vecs);
  return vecs;
}

double
RefSymmSCMatrix::solve_lin(const RefSCVector& v) const
{
  require_nonnull();
  
  RefSymmSCMatrix ret = clone();
  ret->assign(pointer());
  return ret->solve_this(v.pointer());
}

double
RefSymmSCMatrix::determ() const
{
  require_nonnull();
  
  RefSymmSCMatrix ret = clone();
  ret->assign(pointer());
  return ret->determ_this();
}

double
RefSymmSCMatrix::trace() const
{
  require_nonnull();
  return pointer()->trace();
}

void
RefSymmSCMatrix::diagonalize(const RefDiagSCMatrix& vals,
                             const RefSCMatrix& vecs) const
{
  require_nonnull();
  pointer()->diagonalize(vals.pointer(),vecs.pointer());
}

RefSymmSCMatrix
RefSymmSCMatrix::copy() const
{
  if (null()) return 0;
  RefSymmSCMatrix v = kit()->symmmatrix(dim());
  v.assign(*this);
  return v;
}

void
RefSymmSCMatrix::randomize() const
{
  require_nonnull();
  pointer()->randomize();
}

void
RefSymmSCMatrix::assign(const RefSymmSCMatrix&a) const
{
  require_nonnull();
  pointer()->assign(a.pointer());
}

void
RefSymmSCMatrix::assign(const double*v) const
{
  require_nonnull();
  pointer()->assign(v);
}

void
RefSymmSCMatrix::assign(const double**v) const
{
  require_nonnull();
  pointer()->assign(v);
}

void
RefSymmSCMatrix::convert(double*v) const
{
  require_nonnull();
  pointer()->convert(v);
}

void
RefSymmSCMatrix::convert(double**v) const
{
  require_nonnull();
  pointer()->convert(v);
}

void
RefSymmSCMatrix::scale(double a) const
{
  require_nonnull();
  pointer()->scale(a);
}

void
RefSymmSCMatrix::assign(double a) const
{
  require_nonnull();
  pointer()->assign(a);
}

void
RefSymmSCMatrix::accumulate(const RefSymmSCMatrix&a) const
{
  require_nonnull();
  pointer()->accumulate(a.pointer());
}

void
RefSymmSCMatrix::element_op(const Ref<SCElementOp>&op) const
{
  if (nonnull()) pointer()->element_op(op);
}

void
RefSymmSCMatrix::element_op(const Ref<SCElementOp2>&op,
                            const RefSymmSCMatrix&m) const
{
  if (nonnull()) pointer()->element_op(op,m.pointer());
}

void
RefSymmSCMatrix::element_op(const Ref<SCElementOp3>&op,
                            const RefSymmSCMatrix&m,
                            const RefSymmSCMatrix&n) const
{
  if (nonnull()) pointer()->element_op(op,m.pointer(),n.pointer());
}

void
RefSymmSCMatrix::print(ostream& out) const
{
  print(0,out);
}

void
RefSymmSCMatrix::print(const char*title,ostream&out, int precision) const
{
  if (nonnull()) {
      pointer()->print(title,out,precision);
    }
  else {
      if (title) out << endl << title << endl;
      out << "null matrix" << endl;
    }
}

RefSCMatrix
RefSymmSCMatrix::operator*(const RefSCMatrix&a) const
{
  require_nonnull();
  a.require_nonnull();

  RefSCMatrix r = kit()->matrix(dim(),a->coldim());
  r->assign(0.0);
  r->accumulate_product(pointer(),a.pointer());
  return r;
}

RefSCMatrix
RefSymmSCMatrix::operator*(const RefSymmSCMatrix&a) const
{
  require_nonnull();
  a.require_nonnull();

  RefSCMatrix r = kit()->matrix(dim(),a->dim());
  r->assign(0.0);
  r->accumulate_product(pointer(),a.pointer());
  return r;
}

RefSCVector
RefSymmSCMatrix::operator*(const RefSCVector&a) const
{
  require_nonnull();
  a.require_nonnull();

  RefSCVector r = kit()->vector(dim());
  r->assign(0.0);
  r->accumulate_product(pointer(),a.pointer());
  return r;
}

RefSymmSCMatrix
RefSymmSCMatrix::operator *(double a) const
{
  RefSymmSCMatrix r(copy());
  r.scale(a);
  return r;
}

RefSymmSCMatrix
operator *(double a, const RefSymmSCMatrix& v)
{
  return v*a;
}

void
RefSymmSCMatrix::accumulate_symmetric_outer_product(const RefSCVector&v) const
{
  require_nonnull();
  pointer()->accumulate_symmetric_outer_product(v.pointer());
}

double
RefSymmSCMatrix::scalar_product(const RefSCVector&v) const
{
  if (null()) return 0.0;
  return pointer()->scalar_product(v.pointer());
}

void
RefSymmSCMatrix::save(StateOut&s)
{
  if (null()) s.put(0);
  else {
      s.put(1);
      pointer()->save(s);
    }
}

void
RefSymmSCMatrix::restore(StateIn&s)
{
  int have_matrix;
  s.get(have_matrix);
  if (have_matrix && nonnull()) {
      pointer()->restore(s);
    }
  else if (have_matrix) {
      ExEnv::errn() << "RefSymmSCMatrix::restore: "
           << "matrix not properly initialized" << endl;
      abort();
    }
  else {
      clear();
    }
}

int
RefSymmSCMatrix::nblock() const
{
  BlockedSymmSCMatrix *b = dynamic_cast<BlockedSymmSCMatrix*>(pointer());
  if (b) return b->nblocks();
  return 1;
}

RefSymmSCMatrix
RefSymmSCMatrix::block(int i) const
{
  BlockedSymmSCMatrix *b = dynamic_cast<BlockedSymmSCMatrix*>(pointer());
  if (b) return b->block(i);
  return *this;
}


///////////////////////////////////////////////////////////////////
// RefDiagSCMatrix members

RefDiagSCMatrix::RefDiagSCMatrix()
{
}
             
RefDiagSCMatrix::RefDiagSCMatrix (const RefDiagSCMatrix & o):
  Ref<DiagSCMatrix> (o)
{
}
             
RefDiagSCMatrix::RefDiagSCMatrix (DiagSCMatrix * o):
  Ref<DiagSCMatrix> (o)
{
}

RefDiagSCMatrix::~RefDiagSCMatrix ()
{
}

RefDiagSCMatrix&
RefDiagSCMatrix::operator=(DiagSCMatrix* cr)
{
  Ref<DiagSCMatrix>::operator=(cr);
  return *this;
}

RefDiagSCMatrix&
RefDiagSCMatrix::operator=(const RefDiagSCMatrix & c)
{
  Ref<DiagSCMatrix>::operator=(c);
  return *this;
}

RefDiagSCMatrix::RefDiagSCMatrix(const RefSCDimension&a,
                                 const Ref<SCMatrixKit>&k)
{
  a.require_nonnull();
  assign_pointer(k->diagmatrix(a));
}

void
RefDiagSCMatrix::set_element(int i, double a) const
{
  require_nonnull();
  pointer()->set_element(i,a);
}

void
RefDiagSCMatrix::accumulate_element(int i, double a) const
{
  require_nonnull();
  pointer()->accumulate_element(i,a);
}

double
RefDiagSCMatrix::get_element(int i) const
{
  require_nonnull();