replsymm.cc

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

              rect[i][j] = rows[c][j];
          for (; j < n; j++)
              rect[i][j] = rows[j][c];
        }
  
      dist_diagonalize(n, nvec, rect[0], eigvals, eigvecs[0], messagegrp());

      la->assign(eigvals);
      delete[] eigvals;

      int *tivec = new int [mvec];
      for (i=0; i < nproc; i++) {
          int tnvec;
          
          if (i==me) {
              messagegrp()->bcast(nvec, me);
              messagegrp()->bcast(eigvecs[0], n*nvec, me);
              messagegrp()->bcast(ivec, nvec, me);
              tnvec = nvec;
              memcpy(tivec, ivec, sizeof(int)*nvec);
              memcpy(rect[0], eigvecs[0], sizeof(double)*n*nvec);
            }
          else {
              messagegrp()->bcast(tnvec, i);
              messagegrp()->bcast(rect[0], n*tnvec, i);
              messagegrp()->bcast(tivec, tnvec, i);
            }

          for (int j=0; j < tnvec; j++) {
              int c = tivec[j];
              for (int k=0; k < n; k++)
                  lb->rows[k][c] = rect[j][k];
            }
        }

      delete[] ivec;
      delete[] tivec;
      cmat_delete_matrix(eigvecs);
      cmat_delete_matrix(rect);
    }
}

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

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

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

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

  if (!dim()->equiv(la->rowdim()) || !dim()->equiv(la->coldim())) {
      ExEnv::errn() << indent << "ReplSymmSCMatrix::"
           << "accumulate_symmetric_sum(SCMatrix*a): bad dim\n";
      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
ReplSymmSCMatrix::accumulate_symmetric_outer_product(SCVector*a)
{
  // make sure that the argument is of the correct type
  ReplSCVector* la
    = require_dynamic_cast<ReplSCVector*>(a,"ReplSymmSCMatrix::"
                                      "accumulate_symmetric_outer_product");

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

  int n = dim().n();
  double** tdat = this->rows;
  double* adat = la->vector;
  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
ReplSymmSCMatrix::accumulate_transform(SCMatrix*a,SymmSCMatrix*b,
                                       SCMatrix::Transform t)
{
  int i,j,k;
  int ii,jj;
  int nc, nr;

  // do the necessary castdowns
  ReplSCMatrix*la
    = require_dynamic_cast<ReplSCMatrix*>(a,"%s::accumulate_transform",
                                      class_name());
  ReplSymmSCMatrix*lb = require_dynamic_cast<ReplSymmSCMatrix*>(
      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 << "ReplSymmSCMatrix::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 << "ReplSymmSCMatrix::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);

  // if one processor then minimize the amount of memory used
  if (nproc == 1) {
    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);
  }

  // this version requires a full temp matrix be kept
  else {
    int me = messagegrp()->me();
    int mod = nr%nproc;
    int njrow = nr/nproc + ((mod <= me) ? 0 : 1);
    int jstart = (nr/nproc)*me + ((mod <= me) ? mod : me);
    int jend = jstart+njrow;

    double **temp = cmat_new_rect_matrix(nr,nc);
    memset(temp[0], 0, sizeof(double)*nr*nc);
    
    for (i=0; i < nc; i += D1) {
      int ni = nc-i;
      if (ni > D1) ni = D1;

      for (k=0; k < nc; k += D1) {
        int nk = nc-k;
        if (nk > D1) nk = D1;
          
        copy_sym_block(ablock, lb->rows, i, ni, k, nk);

        for (j=jstart; j < jend; j += D1) {
          int nj = jend-j;
          if (nj > D1) nj = D1;
          
          if (t == SCMatrix::NormalTransform)
            copy_block(bblock, la->rows, j, nj, k, nk);
          else
            copy_trans_block(bblock, la->rows, j, nj, k, nk);

          memset(cblock[0], 0, sizeof(double)*D1*D1);
          mult_block(ablock, bblock, cblock, ni, nj, nk);

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

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

    int ind=0;
    for (i=0; i < nr; i += D1) {
      int ni = nr-i;
      if (ni > D1) ni = D1;
        
      for (j=0; j <= i; j += D1, ind++) {
        if (ind%nproc != me)
          continue;
        
        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;
            
          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_block(bblock, temp, j, nj, k, nk);
          mult_block(ablock, bblock, cblock, ni, nj, nk);
        }

        if (i==j) {
          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];
        }
      }
    }
    
    ind=0;
    for (i=0; i < nr; i += D1) {
      int ni = nr-i;
      if (ni > D1) ni = D1;

      for (j=0; j <= i; j += D1, ind++) {
        int nj = nr-j;
        if (nj > D1) nj = D1;

        int proc = ind%nproc;

        if (proc==me)
          copy_sym_block(ablock, rows, i, ni, j, nj);
          
        messagegrp()->bcast(ablock[0], D1*D1, proc);

        if (i==j) {
          for (ii=0; ii < ni; ii++)
            for (jj=0; jj <= ii; jj++)
              rows[i+ii][j+jj] = ablock[ii][jj];
        } else {
          for (ii=0; ii < ni; ii++)
            for (jj=0; jj < nj; jj++)
              rows[i+ii][j+jj] = ablock[ii][jj];
        }
      }
    }

    cmat_delete_matrix(temp);
  }

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

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

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

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

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

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

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

  int nelem = n();
  double* adat = la->vector;
  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
ReplSymmSCMatrix::element_op(const Ref<SCElementOp>& op)
{
  if (op->has_side_effects()) before_elemop();
  scmat_perform_op_on_blocks(op, blocklist);
  if (op->has_side_effects()) after_elemop();
  if (op->has_collect()) op->collect(messagegrp());
}

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

  if (!dim()->equiv(lm->dim())) {
      ExEnv::errn() << indent << "ReplSymmSCMatrix: bad element_op\n";
      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
ReplSymmSCMatrix::element_op(const Ref<SCElementOp3>& op,
                              SymmSCMatrix* m,SymmSCMatrix* n)
{
  ReplSymmSCMatrix *lm
      = require_dynamic_cast<ReplSymmSCMatrix*>(m,"ReplSymSCMatrix::element_op");
  ReplSymmSCMatrix *ln
      = require_dynamic_cast<ReplSymmSCMatrix*>(n,"ReplSymSCMatrix::element_op");

  if (!dim()->equiv(lm->dim()) || !dim()->equiv(ln->dim())) {
      ExEnv::errn() << indent << "ReplSymmSCMatrix: bad element_op\n";
      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 (with a bit of hacking)
void
ReplSymmSCMatrix::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 (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 << setw(lwidth) << i;
    os << endl;

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

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

    ii=kk;
  }
}

Ref<SCMatrixSubblockIter>
ReplSymmSCMatrix::local_blocks(SCMatrixSubblockIter::Access access)
{
  return new ReplSCMatrixListSubblockIter(access, blocklist,
                                          messagegrp(),
                                          matrix, (d->n()*(d->n()+1))/2);
}

Ref<SCMatrixSubblockIter>
ReplSymmSCMatrix::all_blocks(SCMatrixSubblockIter::Access access)
{
  if (access == SCMatrixSubblockIter::Write) {
      ExEnv::errn() << indent << "ReplSymmSCMatrix::all_blocks: "
           << "Write access permitted for local blocks only"
           << endl;
      abort();
    }
  Ref<SCMatrixBlockList> allblocklist = new SCMatrixBlockList();
  allblocklist->insert(new SCMatrixLTriSubBlock(0, d->n(),
                                                0, d->n(), matrix));
  return new ReplSCMatrixListSubblockIter(access, allblocklist,
                                          messagegrp(),
                                          matrix, (d->n()*(d->n()+1))/2);
}

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

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

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