replsymm.cc

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

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

#include <util/misc/formio.h>
#include <util/keyval/keyval.h>
#include <math/scmat/repl.h>
#include <math/scmat/cmatrix.h>
#include <math/scmat/elemop.h>
#include <math/scmat/disthql.h>
#include <math/scmat/offset.h>
#include <math/scmat/mops.h>
#include <math/scmat/util.h>

using namespace std;
using namespace sc;

/////////////////////////////////////////////////////////////////////////////
// ReplSymmSCMatrix member functions

static ClassDesc ReplSymmSCMatrix_cd(
  typeid(ReplSymmSCMatrix),"ReplSymmSCMatrix",1,"public SymmSCMatrix",
  0, 0, 0);

static double **
init_symm_rows(double *data, int n)
{
  double** r = new double*[n];
  for (int i=0; i<n; i++) r[i] = &data[(i*(i+1))/2];
  return r;
}

ReplSymmSCMatrix::ReplSymmSCMatrix(const RefSCDimension&a,ReplSCMatrixKit*k):
  SymmSCMatrix(a,k),
  rows(0)
{
  int n = d->n();

  matrix = new double[n*(n+1)>>1];
  rows = init_symm_rows(matrix,n);

  init_blocklist();
}

void
ReplSymmSCMatrix::before_elemop()
{
  // zero out the blocks not in my block list
  int i, j, index;
  int nproc = messagegrp()->n();
  int me = messagegrp()->me();
  for (i=0, index=0; i<d->blocks()->nblock(); i++) {
      for (j=0; j<=i; j++, index++) {
          if (index%nproc == me) continue;
          for (int ii=d->blocks()->start(i); ii<d->blocks()->fence(i); ii++) {
              for (int jj=d->blocks()->start(j);
                   jj < d->blocks()->fence(j) && jj <= ii;
                   jj++) {
                  matrix[(ii*(ii+1)>>1) + jj] = 0.0;
                }
            }
        }
    }
}

void
ReplSymmSCMatrix::after_elemop()
{
  messagegrp()->sum(matrix, d->n()*(d->n()+1)>>1);
}

void
ReplSymmSCMatrix::init_blocklist()
{
  int i, j, index;
  int nproc = messagegrp()->n();
  int me = messagegrp()->me();
  blocklist = new SCMatrixBlockList;
  for (i=0, index=0; i<d->blocks()->nblock(); i++) {
      for (j=0; j<=i; j++, index++) {
          if (index%nproc != me) continue;
          blocklist->insert(
              new SCMatrixLTriSubBlock(d->blocks()->start(i),
                                       d->blocks()->fence(i),
                                       d->blocks()->start(j),
                                       d->blocks()->fence(j),
                                       matrix));
        }
    }
}

ReplSymmSCMatrix::~ReplSymmSCMatrix()
{
  if (matrix) delete[] matrix;
  if (rows) delete[] rows;
}

int
ReplSymmSCMatrix::compute_offset(int i,int j) const
{
  if (i<0 || j<0 || i>=d->n() || j>=d->n()) {
      ExEnv::errn() << indent << "ReplSymmSCMatrix: index out of bounds\n";
      abort();
    }
  return ij_offset(i,j);
}

double
ReplSymmSCMatrix::get_element(int i,int j) const
{
  return matrix[compute_offset(i,j)];
}

void
ReplSymmSCMatrix::set_element(int i,int j,double a)
{
  matrix[compute_offset(i,j)] = a;
}

void
ReplSymmSCMatrix::accumulate_element(int i,int j,double a)
{
  matrix[compute_offset(i,j)] += a;
}

SCMatrix *
ReplSymmSCMatrix::get_subblock(int br, int er, int bc, int ec)
{
  int nsrow = er-br+1;
  int nscol = ec-bc+1;

  if (nsrow > n() || nscol > n()) {
    ExEnv::errn() << indent
         << "ReplSymmSCMatrix::get_subblock: trying to get too big a "
         << "subblock (" << nsrow << "," << nscol
         << ") from (" << n() << "," << n() << ")\n";
    abort();
  }
  
  RefSCDimension dnrow = (nsrow==n()) ? dim().pointer():new SCDimension(nsrow);
  RefSCDimension dncol = (nscol==n()) ? dim().pointer():new SCDimension(nscol);

  SCMatrix * sb = kit()->matrix(dnrow,dncol);
  sb->assign(0.0);

  ReplSCMatrix *lsb =
    require_dynamic_cast<ReplSCMatrix*>(sb, "ReplSymmSCMatrix::get_subblock");

  for (int i=0; i < nsrow; i++)
    for (int j=0; j < nscol; j++)
      lsb->rows[i][j] = get_element(i+br,j+bc);
      
  return sb;
}

SymmSCMatrix *
ReplSymmSCMatrix::get_subblock(int br, int er)
{
  int nsrow = er-br+1;

  if (nsrow > n()) {
    ExEnv::errn() << indent
         << "ReplSymmSCMatrix::get_subblock: trying to get too big a "
         << "subblock (" << nsrow << "," << nsrow
         << ") from (" << n() << "," << n() << ")\n";
    abort();
  }
  
  RefSCDimension dnrow = new SCDimension(nsrow);

  SymmSCMatrix * sb = kit()->symmmatrix(dnrow);
  sb->assign(0.0);

  ReplSymmSCMatrix *lsb =
    require_dynamic_cast<ReplSymmSCMatrix*>(sb, "ReplSymmSCMatrix::get_subblock");

  for (int i=0; i < nsrow; i++)
    for (int j=0; j <= i; j++)
      lsb->rows[i][j] = get_element(i+br,j+br);
      
  return sb;
}

void
ReplSymmSCMatrix::assign_subblock(SCMatrix*sb, int br, int er, int bc, int ec)
{
  ReplSCMatrix *lsb = require_dynamic_cast<ReplSCMatrix*>(sb,
                                      "ReplSCMatrix::assign_subblock");

  int nsrow = er-br+1;
  int nscol = ec-bc+1;

  if (nsrow > n() || nscol > n()) {
    ExEnv::errn() << indent
         << "ReplSymmSCMatrix::assign_subblock: trying to assign too big a "
         << "subblock (" << nsrow << "," << nscol
         << ") to (" << n() << "," << n() << ")\n";
    abort();
  }
  
  for (int i=0; i < nsrow; i++)
    for (int j=0; j < nscol; j++)
      set_element(i+br,j+bc,lsb->rows[i][j]);
}

void
ReplSymmSCMatrix::assign_subblock(SymmSCMatrix*sb, int br, int er)
{
  ReplSymmSCMatrix *lsb = require_dynamic_cast<ReplSymmSCMatrix*>(sb,
                                      "ReplSymmSCMatrix::assign_subblock");

  int nsrow = er-br+1;

  if (nsrow > n()) {
    ExEnv::errn() << indent
         << "ReplSymmSCMatrix::assign_subblock: trying to assign too big a "
         << "subblock (" << nsrow << "," << nsrow
         << ") to (" << n() << "," << n() << ")\n";
    abort();
  }
  
  for (int i=0; i < nsrow; i++)
    for (int j=0; j <= i; j++)
      set_element(i+br,j+br,lsb->rows[i][j]);
}

void
ReplSymmSCMatrix::accumulate_subblock(SCMatrix*sb, int br, int er, int bc, int ec)
{
  ReplSCMatrix *lsb = require_dynamic_cast<ReplSCMatrix*>(sb,
                                  "ReplSymmSCMatrix::accumulate_subblock");

  int nsrow = er-br+1;
  int nscol = ec-bc+1;

  if (nsrow > n() || nscol > n()) {
    ExEnv::errn() << indent << "ReplSymmSCMatrix::accumulate_subblock: "
         << "trying to accumulate too big a "
         << "subblock (" << nsrow << "," << nscol
         << ") to (" << n() << "," << n() << ")\n";
    abort();
  }
  
  for (int i=0; i < nsrow; i++)
    for (int j=0; j < nscol; j++)
      set_element(i+br,j+br,get_element(i+br,j+br)+lsb->rows[i][j]);
}

void
ReplSymmSCMatrix::accumulate_subblock(SymmSCMatrix*sb, int br, int er)
{
  ReplSCMatrix *lsb = require_dynamic_cast<ReplSCMatrix*>(sb,
                                  "ReplSymmSCMatrix::accumulate_subblock");

  int nsrow = er-br+1;

  if (nsrow > n()) {
    ExEnv::errn() << indent << "ReplSymmSCMatrix::accumulate_subblock: trying to "
         << "accumulate too big a "
         << "subblock (" << nsrow << "," << nsrow
         << ") to (" << n() << "," << n() << ")\n";
    abort();
  }
  
  for (int i=0; i < nsrow; i++)
    for (int j=0; j <= i; j++)
      set_element(i+br,j+br,get_element(i+br,j+br)+lsb->rows[i][j]);
}

SCVector *
ReplSymmSCMatrix::get_row(int i)
{
  if (i >= n()) {
    ExEnv::errn() << indent << "ReplSymmSCMatrix::get_row: trying to get invalid row "
         << i << " max " << n() << endl;
    abort();
  }
  
  SCVector * v = kit()->vector(dim());

  ReplSCVector *lv =
    require_dynamic_cast<ReplSCVector*>(v, "ReplSymmSCMatrix::get_row");

  for (int j=0; j < n(); j++)
    lv->set_element(j,get_element(i,j));
      
  return v;
}

void
ReplSymmSCMatrix::assign_row(SCVector *v, int i)
{
  if (i >= n()) {
    ExEnv::errn() << indent
         << "ReplSymmSCMatrix::assign_row: trying to assign invalid row "
         << i << " max " << n() << endl;
    abort();
  }
  
  if (v->n() != n()) {
    ExEnv::errn() << indent << "ReplSymmSCMatrix::assign_row: vector is wrong size, "
         << "is " << v->n() << ", should be " << n() << endl;
    abort();
  }
  
  ReplSCVector *lv =
    require_dynamic_cast<ReplSCVector*>(v, "ReplSymmSCMatrix::assign_row");

  for (int j=0; j < n(); j++)
    set_element(i,j,lv->get_element(j));
}

void
ReplSymmSCMatrix::accumulate_row(SCVector *v, int i)
{
  if (i >= n()) {
    ExEnv::errn() << indent
         << "ReplSymmSCMatrix::accumulate_row: trying to assign invalide row "
         << i << " max " << n() << endl;
    abort();
  }
  
  if (v->n() != n()) {
    ExEnv::errn() << indent
         << "ReplSymmSCMatrix::accumulate_row: vector is wrong size, "
         << "is " << v->n() << ", should be " << n() << endl;
    abort();
  }
  
  ReplSCVector *lv =
    require_dynamic_cast<ReplSCVector*>(v, "ReplSymmSCMatrix::accumulate_row");

  for (int j=0; j < n(); j++)
    set_element(i,j,get_element(i,j)+lv->get_element(j));
}

void
ReplSymmSCMatrix::assign_val(double val)
{
  int n = (d->n()*(d->n()+1))/2;
  for (int i=0; i<n; i++) matrix[i] = val;
}

void
ReplSymmSCMatrix::assign_s(SymmSCMatrix*m)
{
  ReplSymmSCMatrix* lm = dynamic_cast<ReplSymmSCMatrix*>(m);
  if (lm && dim()->equiv(lm->dim())) {
      int d = i_offset(n());
      memcpy(matrix, lm->matrix, sizeof(double)*d);
    }
  else
      SymmSCMatrix::assign(m);
}

void
ReplSymmSCMatrix::assign_p(const double*m)
{
  int d = i_offset(n());
  memcpy(matrix, m, sizeof(double)*d);
}

void
ReplSymmSCMatrix::assign_pp(const double**m)
{
  for (int i=0; i < n(); i++)
      for (int j=0; j <= i; j++)
          rows[i][j] = m[i][j];
}

void
ReplSymmSCMatrix::scale(double s)
{
  int nelem = i_offset(n());
  for (int i=0; i < nelem; i++) matrix[i] *= s;
}

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

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

  int nelem = (this->n() * (this->n() + 1))/2;
  for (int i=0; i<nelem; i++) matrix[i] += la->matrix[i];
}

double
ReplSymmSCMatrix::invert_this()
{
  if (messagegrp()->n() == 1)
      return cmat_invert(rows,1,n());
  else {
      RefDiagSCMatrix refa = kit()->diagmatrix(d);
      RefSCMatrix refb = kit()->matrix(d,d);
      diagonalize(refa.pointer(),refb.pointer());
      double determ = 1.0;
      for (int i=0; i<dim()->n(); i++) {
          double val = refa->get_element(i);
          determ *= val;
        }
      Ref<SCElementOp> op = new SCElementInvert(1.0e-12);
      refa->element_op(op.pointer());
      assign(0.0);
      accumulate_transform(refb.pointer(), refa.pointer());
      return determ;
    }
}

double
ReplSymmSCMatrix::determ_this()
{
  if (messagegrp()->n() == 1)
      return cmat_determ(rows,1,n());
  else {
      RefDiagSCMatrix refa = kit()->diagmatrix(d);
      RefSCMatrix refb = kit()->matrix(d,d);
      diagonalize(refa.pointer(),refb.pointer());
      double determ = 1.0;
      for (int i=0; i<dim()->n(); i++) {
          double val = refa->get_element(i);
          determ *= val;
        }
      return determ;
    }
}

double
ReplSymmSCMatrix::trace()
{
  double ret=0;
  for (int i=0; i < n(); i++) ret += rows[i][i];
  return ret;
}

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

  return cmat_solve_lin(rows,1,lv->vector,n());
}

void
ReplSymmSCMatrix::gen_invert_this()
{
  RefSCMatrix evecs = kit()->matrix(dim(),dim());
  RefDiagSCMatrix evals = kit()->diagmatrix(dim());

  const char *name = "ReplSymmSCMatrix::gen_invert_this";
  ReplDiagSCMatrix *levals = require_dynamic_cast<ReplDiagSCMatrix*>(evals,name);
  ReplSCMatrix *levecs = require_dynamic_cast<ReplSCMatrix*>(evecs,name);

  this->diagonalize(evals.pointer(), evecs.pointer());

  for (int i=0; i < n(); i++) {
      if (fabs(levals->matrix[i]) > 1.0e-8)
          levals->matrix[i] = 1.0/levals->matrix[i];
      else
          levals->matrix[i] = 0;
    }

  assign(0.0);
  accumulate_transform(levecs, levals);
}

void
ReplSymmSCMatrix::diagonalize(DiagSCMatrix*a,SCMatrix*b)
{
  int i;
  
  const char* name = "ReplSymmSCMatrix::diagonalize";
  // make sure that the arguments is of the correct type
  ReplDiagSCMatrix* la = require_dynamic_cast<ReplDiagSCMatrix*>(a,name);
  ReplSCMatrix* lb = require_dynamic_cast<ReplSCMatrix*>(b,name);

  if (!dim()->equiv(la->dim()) ||
      !dim()->equiv(lb->coldim()) || !dim()->equiv(lb->rowdim())) {
      ExEnv::errn() << indent
           << "ReplSymmSCMatrix::diagonalize(DiagSCMatrix*a,SCMatrix*b): "
           << "bad dims\n";
      abort();
    }

  // This sets up the index list of columns to be stored on this node
  int n = dim()->n();
  int me = messagegrp()->me();
  int nproc = messagegrp()->n();

  // if there are fewer vectors than processors, do serial diagonalization
  if (n < nproc || nproc==1) {
      double *eigvals = la->matrix;
      double **eigvecs = lb->rows;
      cmat_diag(rows, eigvals, eigvecs, n, 1, 1.0e-15);
    }
  else {
      int nvec = n/nproc + (me<(n%nproc)?1:0);
      int mvec = n/nproc + ((n%nproc) ? 1 : 0);
      
      int *ivec = new int[nvec];
      for (i=0; i<nvec; i++)
          ivec[i] = i*nproc + me;

      double *eigvals = new double[n];
      double **eigvecs = cmat_new_rect_matrix(mvec, n);
      double **rect = cmat_new_rect_matrix(mvec, n);

      lb->assign(0.0);

      for (i=0; i < nvec; i++) {
          int c = ivec[i];
          int j;
          for (j=0; j <= c; j++)