localsymm.cc

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

  delete[] evals;
  cmat_delete_matrix(evecs);  
}

void
LocalSymmSCMatrix::diagonalize(DiagSCMatrix*a,SCMatrix*b)
{
  if (n() == 0) return;

  const char* name = "LocalSymmSCMatrix::diagonalize";
  // make sure that the arguments is of the correct type
  LocalDiagSCMatrix* la = require_dynamic_cast<LocalDiagSCMatrix*>(a,name);
  LocalSCMatrix* lb = require_dynamic_cast<LocalSCMatrix*>(b,name);

  if (!dim()->equiv(la->dim()) ||
      !dim()->equiv(lb->coldim()) || !dim()->equiv(lb->rowdim())) {
      ExEnv::errn() << indent
           << "LocalSymmSCMatrix::"
           << "diagonalize(DiagSCMatrix*a,SCMatrix*b): bad dims";
      abort();
    }

  double *eigvals;
  double **eigvecs;
  if (!la) {
      eigvals = new double[n()];
    }
  else {
      eigvals = la->block->data;
    }

  if (!lb) {
      eigvecs = cmat_new_square_matrix(n());
    }
  else {
      eigvecs = lb->rows;
    }

  cmat_diag(rows,eigvals,eigvecs,n(),1,1.0e-15);

  if (!la) delete[] eigvals;
  if (!lb) cmat_delete_matrix(eigvecs);
}

// computes this += a * a.t
void
LocalSymmSCMatrix::accumulate_symmetric_product(SCMatrix*a)
{
  // make sure that the argument is of the correct type
  LocalSCMatrix* la
    = require_dynamic_cast<LocalSCMatrix*>(a,"LocalSymmSCMatrix::"
                                          "accumulate_symmetric_product");

  if (!dim()->equiv(la->rowdim())) {
      ExEnv::errn() << indent
           << "LocalSymmSCMatrix::"
           << "accumulate_symmetric_product(SCMatrix*a): bad dim";
      abort();
    }

  cmat_symmetric_mxm(rows,n(),la->rows,la->ncol(),1);
}

// computes this += a + a.t
void
LocalSymmSCMatrix::accumulate_symmetric_sum(SCMatrix*a)
{
  // make sure that the argument is of the correct type
  LocalSCMatrix* la
    = require_dynamic_cast<LocalSCMatrix*>(a,"LocalSymmSCMatrix::"
                                          "accumulate_symmetric_sum");

  if (!dim()->equiv(la->rowdim()) || !dim()->equiv(la->coldim())) {
      ExEnv::errn() << indent
           << "LocalSymmSCMatrix::"
           << "accumulate_symmetric_sum(SCMatrix*a): bad dim";
      abort();
    }

  int n = dim().n();
  double** tdat = this->rows;
  double** adat = la->rows;
  for (int i=0; i<n; i++) {
      for (int j=0; j<=i; j++) {
          tdat[i][j] += adat[i][j] + adat[j][i];
        }
    }
}

void
LocalSymmSCMatrix::accumulate_symmetric_outer_product(SCVector*a)
{
  // make sure that the argument is of the correct type
  LocalSCVector* la
    = require_dynamic_cast<LocalSCVector*>(a,"LocalSymmSCMatrix::"
                                      "accumulate_symmetric_outer_product");

  if (!dim()->equiv(la->dim())) {
      ExEnv::errn() << indent
           << "LocalSymmSCMatrix::"
           << "accumulate_symmetric_outer_product(SCMatrix*a): bad dim";
      abort();
    }

  int n = dim().n();
  double** tdat = this->rows;
  double* adat = la->block->data;
  for (int i=0; i<n; i++) {
      for (int j=0; j<=i; j++) {
          tdat[i][j] += adat[i]*adat[j];
        }
    }
}

// this += a * b * transpose(a)
void
LocalSymmSCMatrix::accumulate_transform(SCMatrix*a,SymmSCMatrix*b,
                                       SCMatrix::Transform t)
{
  int i,j,k;
  int ii,jj;
  int nc, nr;

  // do the necessary castdowns
  LocalSCMatrix*la
    = require_dynamic_cast<LocalSCMatrix*>(a,"%s::accumulate_transform",
                                      class_name());
  LocalSymmSCMatrix*lb = require_dynamic_cast<LocalSymmSCMatrix*>(
      b,"%s::accumulate_transform", class_name());

  // check the dimensions
  if (t == SCMatrix::NormalTransform) {
    if (!dim()->equiv(la->rowdim()) || !lb->dim()->equiv(la->coldim())) {
      ExEnv::errn() << indent << "LocalSymmSCMatrix::accumulate_transform: bad dim\n";
      abort();
    }

    nc = lb->n();
    nr = la->nrow();
  } else {
    if (!dim()->equiv(la->coldim()) || !lb->dim()->equiv(la->rowdim())) {
      ExEnv::errn() << indent << "LocalSymmSCMatrix::accumulate_transform: bad dim\n";
      abort();
    }

    nc = lb->n();
    nr = la->ncol();
  }

  if (nr==0 || nc==0)
    return;
  
  int nproc = messagegrp()->n();

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

  double **temp = cmat_new_rect_matrix(D1,nc);

  for (i=0; i < nr; i += D1) {
      int ni = nr-i;
      if (ni > D1) ni = D1;

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

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

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

              if (t == SCMatrix::NormalTransform)
                  copy_block(ablock, la->rows, i, ni, k, nk);
              else
                  copy_trans_block(ablock, la->rows, i, ni, k, nk);
          
              copy_sym_block(bblock, lb->rows, j, nj, k, nk);
              copy_block(cblock, temp, 0, ni, j, nj);
              mult_block(ablock, bblock, cblock, ni, nj, nk);
              return_block(temp, cblock, 0, ni, j, nj);
            }
        }

      // now do ab * a~
      for (j=0; j <= i; j+= D1) {
          int nj = nr-j;
          if (nj > D1) nj = D1;

          memset(cblock[0], 0, sizeof(double)*D1*D1);
      
          for (k=0; k < nc; k += D1) {
        
              int nk = nc-k;
              if (nk > D1) nk = D1;

              copy_block(ablock, temp, 0, ni, k, nk);
              if (t == SCMatrix::NormalTransform)
                  copy_block(bblock, la->rows, j, nj, k, nk);
              else
                  copy_trans_block(bblock, la->rows, j, nj, k, nk);
          
              mult_block(ablock, bblock, cblock, ni, nj, nk);
            }

          // copy cblock(i,j) into result
          if (j==i) {
              for (ii=0; ii < ni; ii++)
                  for (jj=0; jj <= ii; jj++)
                      rows[i+ii][j+jj] += cblock[ii][jj];
            } else {
                for (ii=0; ii < ni; ii++)
                    for (jj=0; jj < nj; jj++)
                        rows[i+ii][j+jj] += cblock[ii][jj];
              }
        }
    }

  cmat_delete_matrix(temp);

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

// this += a * b * transpose(a)
void
LocalSymmSCMatrix::accumulate_transform(SCMatrix*a,DiagSCMatrix*b,
                                        SCMatrix::Transform t)
{
  // do the necessary castdowns
  LocalSCMatrix*la
    = require_dynamic_cast<LocalSCMatrix*>(a,"%s::accumulate_transform",
                                      class_name());
  LocalDiagSCMatrix*lb
    = require_dynamic_cast<LocalDiagSCMatrix*>(b,"%s::accumulate_transform",
                                          class_name());

  // check the dimensions
  if (!dim()->equiv(la->rowdim()) || !lb->dim()->equiv(la->coldim())) {
      ExEnv::errn() << indent
           << "LocalSymmSCMatrix::accumulate_transform: bad dim\n";
      abort();
    }

  cmat_transform_diagonal_matrix(rows,n(),lb->block->data,lb->n(),la->rows,1);
}

void
LocalSymmSCMatrix::accumulate_transform(SymmSCMatrix*a,SymmSCMatrix*b)
{
  SymmSCMatrix::accumulate_transform(a,b);
}

double
LocalSymmSCMatrix::scalar_product(SCVector*a)
{
  // make sure that the argument is of the correct type
  LocalSCVector* la
    = require_dynamic_cast<LocalSCVector*>(a,"LocalSCVector::scalar_product");

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

  int nelem = n();
  double* adat = la->block->data;
  double result = 0.0;
  for (int i=0; i<nelem; i++) {
      for (int j=0; j<i; j++) {
          result += 2.0 * rows[i][j] * adat[i] * adat[j];
        }
      result += rows[i][i] * adat[i] * adat[i];
    }
  return result;
}

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

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

  if (!dim()->equiv(lm->dim())) {
      ExEnv::errn() << indent << "LocalSymmSCMatrix: bad element_op\n";
      abort();
    }
  op->process_spec_ltri(block.pointer(), lm->block.pointer());
}

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

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

// from Ed Seidl at the NIH (with a bit of hacking)
void
LocalSymmSCMatrix::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 (n()==0) {
    os << indent << "empty matrix\n";
    return;
  }

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

    // print column indices
    os << indent;
    for (i=ii; i <= nn; i++)
      os << scprintf("%*d",lwidth,i);
    os << endl;

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

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

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

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

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

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