localsymm.cc

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

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

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

#include <math/scmat/mops.h>

using namespace std;
using namespace sc;

/////////////////////////////////////////////////////////////////////////////
// LocalSymmSCMatrix member functions

static ClassDesc LocalSymmSCMatrix_cd(
  typeid(LocalSymmSCMatrix),"LocalSymmSCMatrix",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;
}

LocalSymmSCMatrix::LocalSymmSCMatrix(const RefSCDimension&a,
                                     LocalSCMatrixKit *kit):
  SymmSCMatrix(a,kit),
  rows(0)
{
  resize(a->n());
}

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

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

void
LocalSymmSCMatrix::resize(int n)
{
  block = new SCMatrixLTriBlock(0,n);
  rows = init_symm_rows(block->data,n);
}

double *
LocalSymmSCMatrix::get_data()
{
  return block->data;
}

double **
LocalSymmSCMatrix::get_rows()
{
  return rows;
}

double
LocalSymmSCMatrix::get_element(int i,int j) const
{
  return block->data[compute_offset(i,j)];
}

void
LocalSymmSCMatrix::set_element(int i,int j,double a)
{
  block->data[compute_offset(i,j)] = a;
}

void
LocalSymmSCMatrix::accumulate_element(int i,int j,double a)
{
  block->data[compute_offset(i,j)] += a;
}

SCMatrix *
LocalSymmSCMatrix::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
         << "LocalSymmSCMatrix::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);

  LocalSCMatrix *lsb =
    require_dynamic_cast<LocalSCMatrix*>(sb, "LocalSymmSCMatrix::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 *
LocalSymmSCMatrix::get_subblock(int br, int er)
{
  int nsrow = er-br+1;

  if (nsrow > n()) {
    ExEnv::errn() << indent
         << "LocalSymmSCMatrix::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);

  LocalSymmSCMatrix *lsb =
    require_dynamic_cast<LocalSymmSCMatrix*>(sb, "LocalSymmSCMatrix::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
LocalSymmSCMatrix::assign_subblock(SCMatrix*sb, int br, int er, int bc, int ec)
{
  LocalSCMatrix *lsb =
    require_dynamic_cast<LocalSCMatrix*>(sb, "LocalSCMatrix::assign_subblock");

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

  if (nsrow > n() || nscol > n()) {
    ExEnv::errn() << indent
         << "LocalSymmSCMatrix::assign_subblock: trying to assign too big a "
         << "subblock (" << nsrow << "," << nscol
         << ") from (" << 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
LocalSymmSCMatrix::assign_subblock(SymmSCMatrix*sb, int br, int er)
{
  LocalSymmSCMatrix *lsb = require_dynamic_cast<LocalSymmSCMatrix*>(sb,
                                        "LocalSymmSCMatrix::assign_subblock");

  int nsrow = er-br+1;

  if (nsrow > n()) {
    ExEnv::errn() << indent
         << "LocalSymmSCMatrix::assign_subblock: trying to assign too big a "
         << "subblock (" << nsrow << "," << nsrow
         << ") from (" << 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
LocalSymmSCMatrix::accumulate_subblock(SCMatrix*sb, int br, int er, int bc, int ec)
{
  LocalSCMatrix *lsb = require_dynamic_cast<LocalSCMatrix*>(sb,
                                  "LocalSymmSCMatrix::accumulate_subblock");

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

  if (nsrow > n() || nscol > n()) {
    ExEnv::errn() << indent
         << "LocalSymmSCMatrix::accumulate_subblock: trying to "
         << "accumulate too big a "
         << "subblock (" << nsrow << "," << nscol
         << ") from (" << 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
LocalSymmSCMatrix::accumulate_subblock(SymmSCMatrix*sb, int br, int er)
{
  LocalSCMatrix *lsb = require_dynamic_cast<LocalSCMatrix*>(sb,
                                  "LocalSymmSCMatrix::accumulate_subblock");

  int nsrow = er-br+1;

  if (nsrow > n()) {
    ExEnv::errn() << indent
         << "LocalSymmSCMatrix::accumulate_subblock: trying to "
         << "accumulate too big a "
         << "subblock (" << nsrow << "," << nsrow
         << ") from (" << 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 *
LocalSymmSCMatrix::get_row(int i)
{
  if (i >= n()) {
    ExEnv::errn() << indent
         << "LocalSymmSCMatrix::get_row: trying to get invalid row "
         << i << " max " << n() << endl;
    abort();
  }
  
  SCVector * v = kit()->vector(dim());

  LocalSCVector *lv =
    require_dynamic_cast<LocalSCVector*>(v, "LocalSymmSCMatrix::get_row");

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

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

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

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

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

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

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

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

double
LocalSymmSCMatrix::invert_this()
{
  return cmat_invert(rows,1,n());
}

double
LocalSymmSCMatrix::determ_this()
{
  return cmat_determ(rows,1,n());
}

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

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

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

void
LocalSymmSCMatrix::gen_invert_this()
{
  if (n() == 0) return;

  double *evals = new double[n()];
  double **evecs = cmat_new_square_matrix(n());
  
  cmat_diag(rows,evals,evecs,n(),1,1.0e-15);

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

  cmat_transform_diagonal_matrix(rows, n(), evals, n(), evecs, 0);