localrect.cc

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

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

  double **cd = rows;
  double **ad = la->rows;
  double *bd = lb->block->data;
  int ni = a->rowdim().n();
  int nj = b->dim().n();
  int i, j;
  for (i=0; i<ni; i++) {
      for (j=0; j<nj; j++) {
          cd[i][j] += ad[i][j]*bd[j];
        }
    }
}

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

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

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

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

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

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

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

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

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

void
LocalSCMatrix::accumulate(const SCVector*a)
{
  // make sure that the arguments is of the correct type
  const LocalSCVector* la
    = require_dynamic_cast<const LocalSCVector*>(a,"LocalSCVector::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 << "LocalSCMatrix::accumulate(SCVector*a): " <<
          "dimensions don't match" << endl;
      abort();
    }

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

void
LocalSCMatrix::transpose_this()
{
  cmat_transpose_matrix(rows,nrow(),ncol());
  delete[] rows;
  rows = new double*[ncol()];
  cmat_matrix_pointers(rows,block->data,ncol(),nrow());
  RefSCDimension tmp = d1;
  d1 = d2;
  d2 = tmp;

  int itmp = block->istart;
  block->istart = block->jstart;
  block->jstart = itmp;

  itmp = block->iend;
  block->iend = block->jend;
  block->jend = itmp;
}

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

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

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

void
LocalSCMatrix::svd_this(SCMatrix *U, DiagSCMatrix *sigma, SCMatrix *V)
{
  LocalSCMatrix* lU =
    require_dynamic_cast<LocalSCMatrix*>(U,"LocalSCMatrix::svd_this");
  LocalSCMatrix* lV =
    require_dynamic_cast<LocalSCMatrix*>(V,"LocalSCMatrix::svd_this");
  LocalDiagSCMatrix* lsigma =
    require_dynamic_cast<LocalDiagSCMatrix*>(sigma,"LocalSCMatrix::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 << "LocalSCMatrix: svd_this: dimension mismatch\n";
      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->block->data[i*n + 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->block->data[i*m + j] = dU[i + j*m];
        }
    }

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

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

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

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

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

void
LocalSCMatrix::schmidt_orthog(SymmSCMatrix *S, int nc)
{
  LocalSymmSCMatrix* lS =
    require_dynamic_cast<LocalSymmSCMatrix*>(S,"LocalSCMatrix::schmidt_orthog");
  
  // make sure that the dimensions match
  if (!rowdim()->equiv(lS->dim())) {
      ExEnv::errn() << indent << "LocalSCMatrix::schmidt_orthog(): " <<
          "dimensions don't match\n";
      abort();
    }

  cmat_schmidt(rows,lS->block->data,nrow(),nc);
}

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

  return cmat_schmidt_tol(rows,lS->block->data,nrow(),ncol(),tol,res);
}

void
LocalSCMatrix::element_op(const Ref<SCElementOp>& op)
{
  op->process_spec_rect(block.pointer());
}

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

  if (!rowdim()->equiv(lm->rowdim()) || !coldim()->equiv(lm->coldim())) {
      ExEnv::errn() << indent << "LocalSCMatrix: bad element_op\n";
      abort();
    }
  op->process_spec_rect(block.pointer(), lm->block.pointer());
}

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

  if (!rowdim()->equiv(lm->rowdim()) || !coldim()->equiv(lm->coldim()) ||
      !rowdim()->equiv(ln->rowdim()) || !coldim()->equiv(ln->coldim())) {
      ExEnv::errn() << indent << "LocalSCMatrix: bad element_op\n";
      abort();
    }
  op->process_spec_rect(block.pointer(),
                        lm->block.pointer(), ln->block.pointer());
}

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

  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));

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

  if (nrow()==0 || ncol()==0) {
    os << indent << "empty matrix\n";
    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 << scprintf("%*d",lwidth,i);
    os << endl;

    // print the rows
    for (i=0; i < nrow() ; i++) {
      os << indent << scprintf("%5d",i+1);
      for (j=ii-1; j < nn; j++)
        os << scprintf("%*.*f",lwidth,prec,rows[i][j]);
      os << endl;
    }
    os << endl;

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

    ii=kk;
  }
}

Ref<SCMatrixSubblockIter>
LocalSCMatrix::local_blocks(SCMatrixSubblockIter::Access access)
{
  if (messagegrp()->n() > 1) {
      ExEnv::errn() << indent
           << "LocalSCMatrix::local_blocks: not valid for local matrices"
           << endl;
      abort();
    }
  Ref<SCMatrixSubblockIter> iter
      = new SCMatrixSimpleSubblockIter(access, block.pointer());
  return iter;
}

Ref<SCMatrixSubblockIter>
LocalSCMatrix::all_blocks(SCMatrixSubblockIter::Access access)
{
  if (access == SCMatrixSubblockIter::Write) {
      ExEnv::errn() << indent << "LocalSCMatrix::all_blocks: "
           << "Write access permitted for local blocks only"
           << endl;
      abort();
    }
  return local_blocks(access);
}

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

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