localrect.cc

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

//
// localrect.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>

using namespace std;
using namespace sc;

extern "C" {
    int sing_(double *q, int *lq, int *iq, double *s, double *p,
              int *lp, int *ip, double *a, int *la, int *m, int *n,
              double *w);
};

/////////////////////////////////////////////////////////////////////////////
// LocalSCMatrix member functions

static ClassDesc LocalSCMatrix_cd(
  typeid(LocalSCMatrix),"LocalSCMatrix",1,"public SCMatrix",
  0, 0, 0);

static double **
init_rect_rows(double *data, int ni,int nj)
{
  double** r = new double*[ni];
  int i;
  for (i=0; i<ni; i++) r[i] = &data[i*nj];
  return r;
}

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

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

int
LocalSCMatrix::compute_offset(int i,int j) const
{
  if (i<0 || j<0 || i>=d1->n() || j>=d2->n()) {
      ExEnv::errn() << indent << "LocalSCMatrix: index out of bounds\n";
      abort();
    }
  return i*(d2->n()) + j;
}

void
LocalSCMatrix::resize(int nr, int nc)
{
  block = new SCMatrixRectBlock(0,nr,0,nc);
  if (rows) delete[] rows;
  rows = init_rect_rows(block->data,nr,nc);
}

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

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

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

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

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

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

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

  RefSCDimension dncol;
  if (nscol==ncol()) dncol = coldim();
  else dncol = new SCDimension(nscol);

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

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

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

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

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

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

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

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

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

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

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

  for (int j=0; j < ncol(); j++)
    lv->set_element(j,rows[i][j]);
      
  return v;
}

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

  for (int j=0; j < ncol(); j++)
    rows[i][j] = lv->get_element(j);
}

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

  for (int j=0; j < ncol(); j++)
    rows[i][j] += lv->get_element(j);
}

SCVector *
LocalSCMatrix::get_column(int i)
{
  if (i >= ncol()) {
      ExEnv::errn() << indent <<
          "LocalSCMatrix::get_column: trying to get invalid column " <<
          i << " max " << ncol() << endl;
      abort();
    }
  
  SCVector * v = kit()->vector(rowdim());

  LocalSCVector *lv =
    require_dynamic_cast<LocalSCVector*>(v, "LocalSCMatrix::get_column");

  for (int j=0; j < nrow(); j++)
    lv->set_element(j,rows[j][i]);
      
  return v;
}

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

  for (int j=0; j < nrow(); j++)
    rows[j][i] = lv->get_element(j);
}

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

  for (int j=0; j < nrow(); j++)
    rows[j][i] += lv->get_element(j);
}

void
LocalSCMatrix::assign_val(double a)
{
  int n = d1->n() * d2->n();
  double *data = block->data;
  for (int i=0; i<n; i++) data[i] = a;
}

void
LocalSCMatrix::accumulate_product_rr(SCMatrix*a,SCMatrix*b)
{
  const char* name = "LocalSCMatrix::accumulate_product";
  // make sure that the arguments are of the correct type
  LocalSCMatrix* la = require_dynamic_cast<LocalSCMatrix*>(a,name);
  LocalSCMatrix* lb = require_dynamic_cast<LocalSCMatrix*>(b,name);

  // make sure that the dimensions match
  if (!rowdim()->equiv(la->rowdim()) || !coldim()->equiv(lb->coldim()) ||
      !la->coldim()->equiv(lb->rowdim())) {
      ExEnv::errn() << indent << "LocalSCMatrix::accumulate_product: bad dim" << endl;
      ExEnv::errn() << indent << "this row and col dimension:" << endl;
      rowdim()->print(ExEnv::errn());
      coldim()->print(ExEnv::errn());
      ExEnv::errn() << indent << "a row and col dimension:" << endl;
      a->rowdim()->print(ExEnv::errn());
      a->coldim()->print(ExEnv::errn());
      ExEnv::errn() << indent << "b row and col dimension:" << endl;
      b->rowdim()->print(ExEnv::errn());
      b->coldim()->print(ExEnv::errn());
      abort();
    }

  cmat_mxm(la->rows, 0,
           lb->rows, 0,
           rows, 0,
           nrow(), la->ncol(), this->ncol(),
           1);
}

// does the outer product a x b.  this must have rowdim() == a->dim() and
// coldim() == b->dim()
void
LocalSCMatrix::accumulate_outer_product(SCVector*a,SCVector*b)
{
  const char* name = "LocalSCMatrix::accumulate_outer_product";
  // make sure that the arguments are of the correct type
  LocalSCVector* la = require_dynamic_cast<LocalSCVector*>(a,name);
  LocalSCVector* lb = require_dynamic_cast<LocalSCVector*>(b,name);

  // make sure that the dimensions match
  if (!rowdim()->equiv(la->dim()) || !coldim()->equiv(lb->dim())) {
      ExEnv::errn() << indent <<
          "LocalSCMatrix::accumulate_outer_product(SCVector*a,SCVector*b): " <<
          "dimensions don't match" << endl;
      abort();
    }

  int nr = a->n();
  int nc = b->n();
  int i, j;
  double* adat = la->block->data;
  double* bdat = lb->block->data;
  double** thisdat = rows;
  for (i=0; i<nr; i++) {
      for (j=0; j<nc; j++) {
          thisdat[i][j] += adat[i] * bdat[j];
        }
    }
}

void
LocalSCMatrix::accumulate_product_rs(SCMatrix*a,SymmSCMatrix*b)
{
  const char* name = "LocalSCMatrix::accumulate_product";
  // make sure that the arguments are of the correct type
  LocalSCMatrix* la = require_dynamic_cast<LocalSCMatrix*>(a,name);
  LocalSymmSCMatrix* lb = require_dynamic_cast<LocalSymmSCMatrix*>(b,name);

  // make sure that the dimensions match
  if (!rowdim()->equiv(la->rowdim()) || !coldim()->equiv(lb->dim()) ||
      !la->coldim()->equiv(lb->dim())) {
      ExEnv::errn() << indent <<
          "LocalSCMatrix::accumulate_product_rs(SCMatrix*a,SymmSCMatrix*b): " <<
          "dimensions don't match" << endl;
      abort();
    }

  double **cd = rows;
  double **ad = la->rows;
  double **bd = lb->rows;
  int ni = a->rowdim().n();
  int njk = b->dim().n();
  int i, j, k;
  for (i=0; i<ni; i++) {