matrix.cc

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

//
// matrix.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 <util/state/stateio.h>
#include <math/scmat/block.h>
#include <math/scmat/blkiter.h>
#include <math/scmat/matrix.h>
#include <math/scmat/blocked.h>

using namespace std;

namespace sc {

/////////////////////////////////////////////////////////////////////////////
// SCDimension reference member functions

RefSCDimension::RefSCDimension() {}
             
RefSCDimension::RefSCDimension (const RefSCDimension & o):
  Ref<SCDimension> (o) {}
             
RefSCDimension::RefSCDimension (SCDimension * o): Ref<SCDimension> (o) {}

RefSCDimension::~RefSCDimension () {}

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

RefSCDimension&
RefSCDimension::operator<<(const RefBase & c)
{
  Ref<SCDimension>::operator<<(c);
  return *this;
}

RefSCDimension&
RefSCDimension::operator<<(RefCount*a)
{
  Ref<SCDimension>::operator<<(a);
  return *this;
}

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

int
RefSCDimension::n() const
{
  int result;
  if (null()) result = 0;
  else result = pointer()->n();
  return result;
}

RefSCDimension::operator int() const
{
  if (null()) return 0;
  return pointer()->n();
}

/////////////////////////////////////////////////////////////////////////////
// SCMatrix reference member functions
RefSCMatrix::RefSCMatrix() {}
             
RefSCMatrix::RefSCMatrix (const RefSCMatrix & o): Ref<SCMatrix> (o) {}
             
RefSCMatrix::RefSCMatrix (SCMatrix * o): Ref<SCMatrix> (o) {}

RefSCMatrix::~RefSCMatrix () {}

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

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

RefSCMatrix::RefSCMatrix(const RefSCDimension&a,const RefSCDimension&b,
                         const Ref<SCMatrixKit>&k)
{
  assign_pointer(k->matrix(a,b));
}

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

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

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

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

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

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

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

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

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

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

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

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

  RefSCMatrix ret(rowdim(),coldim(),kit());
  
  ret->assign(pointer());
  ret->accumulate(a.pointer());

  return ret;
}

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

  RefSCMatrix ret(rowdim(),coldim(),kit());
  
  ret->assign(a.pointer());
  ret->scale(-1.0);
  ret->accumulate(pointer());

  return ret;
}

RefSCMatrix
RefSCMatrix::t() const
{
  require_nonnull();
  
  RefSCMatrix ret;
  ret = clone();
  ret->assign(pointer());
  ret->transpose_this();
  return ret;
}

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

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

int
RefSCMatrix::nrow() const
{
  if (null()) return 0;
  else return pointer()->nrow();
}

int
RefSCMatrix::ncol() const
{
  if (null()) return 0;
  else return pointer()->ncol();
}

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

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

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

SCMatrixdouble
RefSCMatrix::operator()(int i,int j)  const
{
  return SCMatrixdouble(pointer(),i,j);
}

RefSCMatrix
RefSCMatrix::clone() const
{
  RefSCMatrix r = kit()->matrix(rowdim(),coldim());
  return r;
}

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

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

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

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

RefSCVector
RefSCMatrix::get_column(int i) const
{
  require_nonnull();
  
  RefSCVector ret = pointer()->get_column(i);
  return ret;
}

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

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

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

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

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

RefSCMatrix
RefSCMatrix::copy() const
{
  if (null()) return 0;
  RefSCMatrix v = kit()->matrix(rowdim(),coldim());
  v.assign(*this);
  return v;
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

void
RefSCMatrix::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
RefSCMatrix::operator *(double a) const
{
  RefSCMatrix r(copy());
  r.scale(a);
  return r;
}

void
RefSCMatrix::svd(const RefSCMatrix &U,
                 const RefDiagSCMatrix &sigma,
                 const RefSCMatrix &V)
{
  require_nonnull();
  RefSCMatrix c = clone();
  c->assign(pointer());
  c->svd_this(U.pointer(), sigma.pointer(), V.pointer());
}

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

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

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

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

void
RefSCMatrix::accumulate_outer_product(const RefSCVector& v1,
                                      const RefSCVector&v2) const
{
  require_nonnull();
  pointer()->accumulate_outer_product(v1.pointer(),v2.pointer());
}

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

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

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