replrect.cc

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

  for (i=0; i<ni; i++) {
      for (j=0; j<nj; j++) {
          cd[i][j] += ad[i][j]*bd[j];
        }
    }
}

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

  // make sure that the dimensions match
  if (!rowdim()->equiv(la->rowdim()) || !coldim()->equiv(la->coldim())) {
      ExEnv::errn() << indent << "ReplSCMatrix::accumulate(SCMatrix*a): "
           << "dimensions don't match" << endl;
      abort();
    }

  int nelem = this->ncol() * this->nrow();
  int i;
  for (i=0; i<nelem; i++) matrix[i] += la->matrix[i];
}

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

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

  int n = this->ncol();
  double *dat = la->matrix;
  int i, j;
  for (i=0; i<n; i++) {
      for (j=0; j<i; j++) {
          double tmp = *dat;
          matrix[i*n+j] += tmp;
          matrix[j*n+i] += tmp;
          dat++;
        }
      matrix[i*n+i] += *dat++;
    }
}

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

  // make sure that the dimensions match
  if (!rowdim()->equiv(la->dim()) || !coldim()->equiv(la->dim())) {
      ExEnv::errn() << indent << "ReplSCMatrix::accumulate(DiagSCMatrix*a): "
           << "dimensions don't match" << endl;
      abort();
    }

  int n = this->ncol();
  double *dat = la->matrix;
  int i;
  for (i=0; i<n; i++) {
      matrix[i*n+i] += *dat++;
    }
}

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

  // make sure that the dimensions match
  if (!((rowdim()->equiv(la->dim()) && coldim()->n() == 1)
        || (coldim()->equiv(la->dim()) && rowdim()->n() == 1))) {
      ExEnv::errn() << indent << "ReplSCMatrix::accumulate(SCVector*a): "
           << "dimensions don't match" << endl;
      abort();
    }

  int n = this->ncol();
  int i;
  double *dat = la->vector;
  for (i=0; i<n; i++) {
      matrix[i*n+i] += dat[i];
    }
}

void
ReplSCMatrix::transpose_this()
{
  cmat_transpose_matrix(rows,nrow(),ncol());
  delete[] rows;
  rows = new double*[ncol()];
  cmat_matrix_pointers(rows,matrix,ncol(),nrow());
  RefSCDimension tmp = d1;
  d1 = d2;
  d2 = tmp;
  init_blocklist();
}

double
ReplSCMatrix::invert_this()
{
  if (nrow() != ncol()) {
      ExEnv::errn() << indent << "ReplSCMatrix::invert_this: matrix is not square" << endl;
      abort();
    }
  return cmat_invert(rows,0,nrow());
}

double
ReplSCMatrix::determ_this()
{
  if (nrow() != ncol()) {
    ExEnv::errn() << indent << "ReplSCMatrix::determ_this: matrix is not square" << endl;
    abort();
  }
  return cmat_determ(rows,0,nrow());
}

double
ReplSCMatrix::trace()
{
  if (nrow() != ncol()) {
    ExEnv::errn() << indent << "ReplSCMatrix::trace: matrix is not square" << endl;
    abort();
  }
  double ret=0;
  int i;
  for (i=0; i < nrow(); i++)
    ret += rows[i][i];
  return ret;
}

void
ReplSCMatrix::svd_this(SCMatrix *U, DiagSCMatrix *sigma, SCMatrix *V)
{
  ReplSCMatrix* lU =
    require_dynamic_cast<ReplSCMatrix*>(U,"ReplSCMatrix::svd_this");
  ReplSCMatrix* lV =
    require_dynamic_cast<ReplSCMatrix*>(V,"ReplSCMatrix::svd_this");
  ReplDiagSCMatrix* lsigma =
    require_dynamic_cast<ReplDiagSCMatrix*>(sigma,"ReplSCMatrix::svd_this");

  RefSCDimension mdim = rowdim();
  RefSCDimension ndim = coldim();
  int m = mdim.n();
  int n = ndim.n();

  RefSCDimension pdim;
  if (m == n && m == sigma->dim().n())
    pdim = sigma->dim();
  else if (m<n)
    pdim = mdim;
  else
    pdim = ndim;

  int p = pdim.n();

  if (!mdim->equiv(lU->rowdim()) ||
      !mdim->equiv(lU->coldim()) ||
      !ndim->equiv(lV->rowdim()) ||
      !ndim->equiv(lV->coldim()) ||
      !pdim->equiv(sigma->dim())) {
      ExEnv::errn() << indent << "ReplSCMatrix: svd_this: dimension mismatch" << endl;
      abort();
    }

  // form a fortran style matrix for the SVD routines
  double *dA = new double[m*n];
  double *dU = new double[m*m];
  double *dV = new double[n*n];
  double *dsigma = new double[n];
  double *w = new double[(3*p-1>m)?(3*p-1):m];

  int i,j;
  for (i=0; i<m; i++) {
      for (j=0; j<n; j++) {
          dA[i + j*m] = this->rows[i][j];
        }
    }

  int three = 3;

  sing_(dU, &m, &three, dsigma, dV, &n, &three, dA, &m, &m, &n, w);

  for (i=0; i<m; i++) {
      for (j=0; j<m; j++) {
          lU->rows[i][j] = dU[i + j*m];
        }
    }

  for (i=0; i<n; i++) {
      for (j=0; j<n; j++) {
          lV->rows[i][j] = dV[i + j*n];
        }
    }

  for (i=0; i<p; i++) {
      lsigma->matrix[i] = dsigma[i];
    }

  delete[] dA;
  delete[] dU;
  delete[] dV;
  delete[] dsigma;
  delete[] w;
}

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

  return cmat_solve_lin(rows,0,lv->vector,nrow());
}

void
ReplSCMatrix::schmidt_orthog(SymmSCMatrix *S, int nc)
{
  int i,j,ij;
  int m;

  ReplSymmSCMatrix* lS =
    require_dynamic_cast<ReplSymmSCMatrix*>(S,"ReplSCMatrix::schmidt_orthog");
  
  // make sure that the dimensions match
  if (!rowdim()->equiv(lS->dim())) {
      ExEnv::errn() << indent << "ReplSCMatrix::schmidt_orthog(): "
           << "dimensions don't match" << endl;
      abort();
    }

#if 0
  cmat_schmidt(rows,lS->matrix,nrow(),nc);
#else
  int me = messagegrp()->me();
  int nproc = messagegrp()->n();
  int nr = nrow();
  
  double vtmp;
  double *v = new double[nr];
  double *cm = new double[nr];

  double **sblock = cmat_new_square_matrix(D1);
  
  int mod = nc%nproc;
  int ncoli = nc/nproc + (mod <= me ? 0 : 1);
  int cstart = (nc/nproc)*me + (mod <= me ? mod : me);
  int cend = cstart+ncoli;

  // copy my columns to a rows of temp matrix
  double **cols = cmat_new_rect_matrix(ncoli, nr);
  for (i=cstart; i < cend; i++)
    for (j=0; j < nr; j++)
      cols[i-cstart][j] = rows[j][i];
    
  for (m=0; m < nc; m++) {
    // who has this column
    for (i=0; i < nproc; i++) {
      int ni = nc/nproc + (mod <= i ? 0 : 1);
      int csi = (nc/nproc)*i + (mod <= i ? mod : i);
      if (m >= csi && m < csi+ni) {
        if (i==me)
          memcpy(cm, cols[m-csi], sizeof(double)*nr);
        messagegrp()->bcast(cm, nr, i);
        break;
      }
    }
    
    memset(v, 0, sizeof(double)*nr);
    
    for (i=ij=0; i < nr; i += D1) {
      int ni = nr-i;
      if (ni > D1) ni = D1;
      
      for (j=0; j < nr; j += D1, ij++) {
        if (ij%nproc != me)
          continue;

        int nj = nr-j;
        if (nj > D1) nj = D1;
        
        copy_sym_block(sblock, lS->rows, i, ni, j, nj);
        
        for (int ii=0; ii < ni; ii++)
          for (int jj=0; jj < nj; jj++)
            v[i+ii] += cm[j+jj]*sblock[ii][jj];
      }
    }

    messagegrp()->sum(v, nr);

    for (i=0,vtmp=0.0; i < nr; i++)
      vtmp += v[i]*cm[i];

    if (!vtmp) {
      ExEnv::errn() << "cmat_schmidt: bogus" << endl;
      abort();
    }

    if (vtmp < 1.0e-15)
      vtmp = 1.0e-15;

    vtmp = 1.0/sqrt(vtmp);
    
    for (i=0; i < nr; i++) {
      v[i] *= vtmp;
      cm[i] *= vtmp;
    }

    if (m < nc-1) {
      for (i=m+1; i < nc; i++) {
        if (i < cstart)
          continue;
        if (i >= cend)
          break;
        
        double *ci = cols[i-cstart];
        
        for (j=0,vtmp=0.0; j < nr; j++)
          vtmp += v[j] * ci[j];
        for (j=0; j < nr; j++)
          ci[j] -= vtmp * cm[j];
      }
    }

    // if I own cm then put it back into cols
    if (m >= cstart && m < cend)
      memcpy(cols[m-cstart], cm, sizeof(double)*nr);
  }

  // now collect columns again
  for (i=0; i < nproc; i++) {
    int ni = nc/nproc + (mod <= i ? 0 : 1);
    int csi = (nc/nproc)*i + (mod <= i ? mod : i);
    for (j=0; j < ni; j++) {
      if (i==me) {
        messagegrp()->bcast(cols[j], nr, i);
        for (int k=0; k < nr; k++)
          rows[k][j+csi] = cols[j][k];
      }
      else {
        messagegrp()->bcast(cm, nr, i);
        for (int k=0; k < nr; k++)
          rows[k][j+csi] = cm[k];
      }
    }
  }

  cmat_delete_matrix(sblock);
  cmat_delete_matrix(cols);
  delete[] v;
  delete[] cm;
#endif
}

int
ReplSCMatrix::schmidt_orthog_tol(SymmSCMatrix *S, double tol, double *res)
{
  ReplSymmSCMatrix* lS =
    require_dynamic_cast<ReplSymmSCMatrix*>(S,"ReplSCMatrix::schmidt_orthog_tol");
  
  // make sure that the dimensions match
  if (!rowdim()->equiv(lS->dim())) {
      ExEnv::errn() << indent << "ReplSCMatrix::schmidt_orthog_tol(): " <<
          "dimensions don't match" << endl;
      abort();
    }

  int northog;

  if (messagegrp()->me() == 0) {
      northog = cmat_schmidt_tol(rows,lS->matrix,nrow(),ncol(),tol,res);
    }

  // make sure everybody ends up with the same data
  messagegrp()->bcast(northog);
  messagegrp()->bcast(*res);
  for (int i=0; i<nrow(); i++) {
      messagegrp()->bcast(rows[i],ncol());
    }

  return northog;
}

void
ReplSCMatrix::element_op(const Ref<SCElementOp>& op)
{
  if (op->has_side_effects()) before_elemop();
  SCMatrixBlockListIter i;
  for (i = blocklist->begin(); i != blocklist->end(); i++) {
      op->process_base(i.block());
    }
  if (op->has_side_effects()) after_elemop();
  if (op->has_collect()) op->collect(messagegrp());
}

void
ReplSCMatrix::element_op(const Ref<SCElementOp2>& op,
                          SCMatrix* m)
{
  ReplSCMatrix *lm
      = require_dynamic_cast<ReplSCMatrix*>(m,"ReplSCMatrix::element_op");

  if (!rowdim()->equiv(lm->rowdim()) || !coldim()->equiv(lm->coldim())) {
      ExEnv::errn() << indent << "ReplSCMatrix: bad element_op" << endl;
      abort();
    }
  if (op->has_side_effects()) before_elemop();
  if (op->has_side_effects_in_arg()) lm->before_elemop();
  SCMatrixBlockListIter i, j;
  for (i = blocklist->begin(), j = lm->blocklist->begin();
       i != blocklist->end();
       i++, j++) {
      op->process_base(i.block(), j.block());
    }
  if (op->has_side_effects()) after_elemop();
  if (op->has_side_effects_in_arg()) lm->after_elemop();
  if (op->has_collect()) op->collect(messagegrp());
}

void
ReplSCMatrix::element_op(const Ref<SCElementOp3>& op,
                          SCMatrix* m,SCMatrix* n)
{
  ReplSCMatrix *lm
      = require_dynamic_cast<ReplSCMatrix*>(m,"ReplSCMatrix::element_op");
  ReplSCMatrix *ln
      = require_dynamic_cast<ReplSCMatrix*>(n,"ReplSCMatrix::element_op");

  if (!rowdim()->equiv(lm->rowdim()) || !coldim()->equiv(lm->coldim()) ||
      !rowdim()->equiv(ln->rowdim()) || !coldim()->equiv(ln->coldim())) {
      ExEnv::errn() << indent << "ReplSCMatrix: bad element_op" << endl;
      abort();
    }
  if (op->has_side_effects()) before_elemop();
  if (op->has_side_effects_in_arg1()) lm->before_elemop();
  if (op->has_side_effects_in_arg2()) ln->before_elemop();
  SCMatrixBlockListIter i, j, k;
  for (i = blocklist->begin(),
           j = lm->blocklist->begin(),
           k = ln->blocklist->begin();
       i != blocklist->end();
       i++, j++, k++) {
      op->process_base(i.block(), j.block(), k.block());
    }
  if (op->has_side_effects()) after_elemop();
  if (op->has_side_effects_in_arg1()) lm->after_elemop();
  if (op->has_side_effects_in_arg2()) ln->after_elemop();
  if (op->has_collect()) op->collect(messagegrp());
}

// from Ed Seidl at the NIH
void
ReplSCMatrix::vprint(const char *title, ostream& os, int prec) const
{
  int ii,jj,kk,nn;
  int i,j;
  int lwidth,width;
  double max=this->maxabs();

  if (messagegrp()->me() != 0) return;

  max = (max==0.0) ? 1.0 : log10(max);
  if (max < 0.0) max=1.0;

  lwidth = prec + 5 + (int) max;
  width = 75/(lwidth+SCFormIO::getindent(os));

  os.setf(ios::fixed,ios::floatfield); os.precision(prec);
  os.setf(ios::right,ios::adjustfield);

  if (title)
    os << endl << indent << title << endl;
  else
    os << endl;

  if (nrow()==0 || ncol()==0) {
    os << indent << "empty matrix" << endl;
    return;
  }

  for (ii=jj=0;;) {
    ii++; jj++;
    kk=width*jj;
    nn = (ncol() > kk) ? kk : ncol();

    // print column indices
    os << indent;
    for (i=ii; i <= nn; i++)
      os << setw(lwidth) << i;
    os << endl;

    // print the rows
    for (i=0; i < nrow() ; i++) {
      os << setw(5) << i+1;
      for (j=ii-1; j < nn; j++)
        os << setw(lwidth) << rows[i][j];
      os << endl;
    }
    os << endl;

    if (ncol() <= kk) {
      os.flush();
      return;
    }

    ii=kk;
  }
}

Ref<SCMatrixSubblockIter>
ReplSCMatrix::local_blocks(SCMatrixSubblockIter::Access access)
{
  return new ReplSCMatrixListSubblockIter(access, blocklist,
                                          messagegrp(),
                                          matrix, d1->n()*d2->n());
}

Ref<SCMatrixSubblockIter>
ReplSCMatrix::all_blocks(SCMatrixSubblockIter::Access access)
{
  if (access == SCMatrixSubblockIter::Write) {
      ExEnv::errn() << indent << "ReplSCMatrix::all_blocks: "
           << "Write access permitted for local blocks only"
           << endl;
      abort();
    }
  Ref<SCMatrixBlockList> allblocklist = new SCMatrixBlockList();
  allblocklist->insert(new SCMatrixRectSubBlock(0, d1->n(), d1->n(),
                                                0, d2->n(), matrix));
  return new ReplSCMatrixListSubblockIter(access, allblocklist,
                                          messagegrp(),
                                          matrix, d1->n()*d2->n());
}

Ref<ReplSCMatrixKit>
ReplSCMatrix::skit()
{
  return dynamic_cast<ReplSCMatrixKit*>(kit().pointer());
}

/////////////////////////////////////////////////////////////////////////////

// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End: