replrect.cc

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

//
// replrect.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/mops.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);
};

/////////////////////////////////////////////////////////////////////////////
// ReplSCMatrix member functions

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

ReplSCMatrix::ReplSCMatrix(const RefSCDimension&a,const RefSCDimension&b,
                           ReplSCMatrixKit*k):
  SCMatrix(a,b,k)
{
  int nr = a->n();
  int nc = b->n();

  matrix = new double[nr*nc];

  rows = init_rect_rows(matrix,nr,nc);

  init_blocklist();
}

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

void
ReplSCMatrix::after_elemop()
{
  messagegrp()->sum(matrix, d1->n()*d2->n());
}

void
ReplSCMatrix::init_blocklist()
{
  int i, j, index;
  int nc = d2->n();
  int nproc = messagegrp()->n();
  int me = messagegrp()->me();
  blocklist = new SCMatrixBlockList;
  for (i=0, index=0; i<d1->blocks()->nblock(); i++) {
      for (j=0; j<d2->blocks()->nblock(); j++, index++) {
          if (index%nproc == me) {
              blocklist->insert(
                  new SCMatrixRectSubBlock(d1->blocks()->start(i),
                                           d1->blocks()->fence(i),
                                           nc,
                                           d2->blocks()->start(j),
                                           d2->blocks()->fence(j),
                                           matrix));
            }
        }
    }
}

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

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

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

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

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

SCMatrix *
ReplSCMatrix::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 << "ReplSCMatrix::get_subblock: trying to get too big a"
         << "subblock (" << nsrow << "," << nscol
         << ") from (" << nrow() << "," << ncol() << ")" << endl;
    abort();
  }
  
  RefSCDimension dnrow = (nsrow==nrow()) ? rowdim().pointer()
                                         : new SCDimension(nsrow);
  RefSCDimension dncol = (nscol==ncol()) ? coldim().pointer()
                                         : new SCDimension(nscol);

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

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

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

  if (nsrow > nrow() || nscol > ncol()) {
    ExEnv::errn() << indent
         << "ReplSCMatrix::assign_subblock: trying to assign too big a"
         << "subblock (" << nsrow << "," << nscol
         << ") to (" << nrow() << "," << ncol() << ")" << endl;;
    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
ReplSCMatrix::accumulate_subblock(SCMatrix*sb, int br, int er, int bc, int ec,
                                  int source_br, int source_bc)
{
  ReplSCMatrix *lsb = require_dynamic_cast<ReplSCMatrix*>(sb,
                                      "ReplSCMatrix::accumulate_subblock");

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

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

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

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

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

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

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

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

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

  ReplSCVector *lv =
    require_dynamic_cast<ReplSCVector*>(v, "ReplSCMatrix::get_column");

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

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

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

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

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

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

void
ReplSCMatrix::accumulate_product_rr(SCMatrix*a,SCMatrix*b)
{
  const char* name = "ReplSCMatrix::accumulate_product_rr";
  // make sure that the arguments are of the correct type
  ReplSCMatrix* la = require_dynamic_cast<ReplSCMatrix*>(a,name);
  ReplSCMatrix* lb = require_dynamic_cast<ReplSCMatrix*>(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
           << "ReplSCMatrix::accumulate_product_rr(SCMatrix*a,SCMatrix*b): "
           << "dimensions don't match" << endl;
      abort();
    }

#if 0
  cmat_transpose_matrix(lb->rows, la->ncol(), this->ncol());
  double** btrans;
  btrans = new double*[this->ncol()];
  btrans[0] = lb->rows[0];
  cmat_matrix_pointers(btrans,btrans[0],this->ncol(),la->ncol());

  Ref<SCElementOp> op = new SCElementDot(la->rows, btrans, la->ncol());
  element_op(op);

  cmat_transpose_matrix(btrans,this->ncol(),la->ncol());
  delete[] btrans;
#else
  int i,j,k;
  int ii,jj;

  int nr = la->nrow();
  int nc = lb->ncol();
  int ncc = la->ncol();

  if (nr==0 || nc==0 || ncc==0)
    return;
  
  int nproc = messagegrp()->n();
  int me = messagegrp()->me();
  int mod = nr%nproc;
  int nirow = nr/nproc + ((mod <= me) ? 0 : 1);
  int istart = (nr/nproc)*me + ((mod <= me) ? mod : me);
  int iend = istart+nirow;

  double ** ablock = cmat_new_square_matrix(D1);
  double ** bblock = cmat_new_square_matrix(D1);
  double ** cblock = cmat_new_square_matrix(D1);

  for (i=istart; i < iend; i += D1) {
    int ni = iend-i;
    if (ni > D1) ni = D1;
    
    for (j=0; j < nc; j += D1) {
      int nj = nc-j;
      if (nj > D1) nj = D1;

      memset(cblock[0], 0, sizeof(double)*D1*D1);

      for (k=0; k < ncc; k += D1) {
        int nk = ncc-k;
        if (nk > D1) nk = D1;

        copy_block(ablock, la->rows, i, ni, k, nk);
        copy_trans_block(bblock, lb->rows, j, nj, k, nk);
        mult_block(ablock, bblock, cblock, ni, nj, nk);
      }

      for (ii=0; ii < ni; ii++)
        for (jj=0; jj < nj; jj++)
          rows[i+ii][j+jj] += cblock[ii][jj];
    }
  }

  for (i=0; i < nproc; i++) {
    nirow = nr/nproc + ((mod <= i) ? 0 : 1);
    istart = (nr/nproc)*i + ((mod <= i) ? mod : i);
    if (!nirow)
      break;
    messagegrp()->bcast(rows[istart], nirow*nc, i);
  }

  cmat_delete_matrix(ablock);
  cmat_delete_matrix(bblock);
  cmat_delete_matrix(cblock);
        
#endif  
}

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

  // make sure that the dimensions match
  if (!rowdim()->equiv(la->dim()) || !coldim()->equiv(lb->dim())) {
      ExEnv::errn() << indent
           << "ReplSCMatrix::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->vector;
  double* bdat = lb->vector;
  double** thisdat = rows;
  for (i=0; i<nr; i++) {
      for (j=0; j<nc; j++) {
          thisdat[i][j] += adat[i] * bdat[j];
        }
    }
}

void
ReplSCMatrix::accumulate_product_rs(SCMatrix*a,SymmSCMatrix*b)
{
  const char* name = "ReplSCMatrix::accumulate_product_rs";
  // make sure that the arguments are of the correct type
  ReplSCMatrix* la = require_dynamic_cast<ReplSCMatrix*>(a,name);
  ReplSymmSCMatrix* lb = require_dynamic_cast<ReplSymmSCMatrix*>(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
           << "ReplSCMatrix::accumulate_product_rs(SCMatrix*a,SymmSCMatrix*b): "
           << "dimensions don't match" << endl;
      ExEnv::err0() << indent << "rowdim():" << endl;
      rowdim().print();
      ExEnv::err0() << indent << "coldim():" << endl;
      coldim().print();
      ExEnv::err0() << indent << "la->rowdim():" << endl;
      la->rowdim().print();
      ExEnv::err0() << indent << "la->coldim():" << endl;
      la->coldim().print();
      ExEnv::err0() << indent << "lb->dim():" << endl;
      lb->dim().print();
      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++) {
      for (j=0; j<njk; j++) {
          for (k=0; k<=j; k++) {
              cd[i][k] += ad[i][j]*bd[j][k];
            }
          for (; k<njk; k++) {
              cd[i][k] += ad[i][j]*bd[k][j];
            }
        }
    }
}

void
ReplSCMatrix::accumulate_product_rd(SCMatrix*a,DiagSCMatrix*b)
{
  const char* name = "ReplSCMatrix::accumulate_product_rd";
  // make sure that the arguments are of the correct type
  ReplSCMatrix* la = require_dynamic_cast<ReplSCMatrix*>(a,name);
  ReplDiagSCMatrix* lb = require_dynamic_cast<ReplDiagSCMatrix*>(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
           << "ReplSCMatrix::accumulate_product_rd(SCMatrix*a,DiagSCMatrix*b): "
           << "dimensions don't match" << endl;
      abort();
    }

  double **cd = rows;
  double **ad = la->rows;
  double *bd = lb->matrix;
  int ni = a->rowdim().n();
  int nj = b->dim().n();
  int i, j;