Index: ThirdParty/mpqc_open/src/lib/math/scmat/distdiag.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/scmat/distdiag.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/scmat/distdiag.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,359 @@
+//
+// distdiag.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// 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 <math.h>
+
+#include <util/misc/formio.h>
+#include <util/keyval/keyval.h>
+#include <math/scmat/dist.h>
+#include <math/scmat/cmatrix.h>
+#include <math/scmat/elemop.h>
+
+using namespace std;
+using namespace sc;
+
+#define DEBUG 0
+
+/////////////////////////////////////////////////////////////////////////////
+// DistDiagSCMatrix member functions
+
+static ClassDesc DistDiagSCMatrix_cd(
+  typeid(DistDiagSCMatrix),"DistDiagSCMatrix",1,"public DiagSCMatrix",
+  0, 0, 0);
+
+DistDiagSCMatrix::DistDiagSCMatrix(const RefSCDimension&a,DistSCMatrixKit*k):
+  DiagSCMatrix(a,k)
+{
+  init_blocklist();
+}
+
+int
+DistDiagSCMatrix::block_to_node(int i) const
+{
+  return (i)%messagegrp()->n();
+}
+
+Ref<SCMatrixBlock>
+DistDiagSCMatrix::block_to_block(int i) const
+{
+  int offset = i;
+  int nproc = messagegrp()->n();
+
+  if ((offset%nproc) != messagegrp()->me()) return 0;
+
+  SCMatrixBlockListIter I;
+  for (I=blocklist->begin(); I!=blocklist->end(); I++) {
+      if (I.block()->blocki == i)
+          return I.block();
+    }
+
+  ExEnv::errn() << indent
+       << "DistDiagSCMatrix::block_to_block: internal error" << endl;
+  abort();
+  return 0;
+}
+
+double *
+DistDiagSCMatrix::find_element(int i) const
+{
+  int bi, oi;
+  d->blocks()->elem_to_block(i, bi, oi);
+
+  if (DEBUG)
+      ExEnv::outn() << messagegrp()->me() << ": "
+                   << "find_element(" << i << "): "
+                   << "block = " << bi << ", "
+                   << "offset = " << oi
+                   << endl;
+
+  Ref<SCMatrixDiagBlock> blk; blk << block_to_block(bi);
+  if (blk.nonnull()) {
+      if (DEBUG)
+          ExEnv::outn() << messagegrp()->me() << ": ndat = " << blk->ndat() << endl;
+      if (oi >= blk->ndat()) {
+          ExEnv::errn() << messagegrp()->me() << ": DistDiagSCMatrix::find_element"
+               << ": internal error" << endl;
+          abort();
+        }
+      return &blk->dat()[oi];
+    }
+  else {
+      if (DEBUG)
+          ExEnv::outn() << messagegrp()->me() << ": can't find" << endl;
+      return 0;
+    }
+}
+
+int
+DistDiagSCMatrix::element_to_node(int i) const
+{
+  int bi, oi;
+  d->blocks()->elem_to_block(i, bi, oi);
+
+  return block_to_node(bi);
+}
+
+void
+DistDiagSCMatrix::init_blocklist()
+{
+  int i;
+  int nproc = messagegrp()->n();
+  int me = messagegrp()->me();
+  blocklist = new SCMatrixBlockList;
+  SCMatrixBlock *b;
+  for (i=0; i<d->blocks()->nblock(); i++) {
+      if (i%nproc != me) continue;
+      b = new SCMatrixDiagBlock(d->blocks()->start(i),
+                                d->blocks()->fence(i),
+                                d->blocks()->start(i));
+      b->blocki = i;
+      b->blockj = i;
+      blocklist->insert(b);
+    }
+}
+
+DistDiagSCMatrix::~DistDiagSCMatrix()
+{
+}
+
+double
+DistDiagSCMatrix::get_element(int i) const
+{
+  double res;
+  double *e = find_element(i);
+  if (e) {
+      res = *e;
+      messagegrp()->bcast(res, messagegrp()->me());
+    }
+  else {
+      messagegrp()->bcast(res, element_to_node(i));
+    }
+  return res;
+}
+
+void
+DistDiagSCMatrix::set_element(int i,double a)
+{
+  double *e = find_element(i);
+  if (e) {
+      *e = a;
+    }
+}
+
+void
+DistDiagSCMatrix::accumulate_element(int i,double a)
+{
+  double *e = find_element(i);
+  if (e) {
+      *e += a;
+    }
+}
+
+void
+DistDiagSCMatrix::accumulate(const DiagSCMatrix*a)
+{
+  // make sure that the argument is of the correct type
+  const DistDiagSCMatrix* la
+    = require_dynamic_cast<const DistDiagSCMatrix*>(a,"DistDiagSCMatrix::accumulate");
+
+  // make sure that the dimensions match
+  if (!dim()->equiv(la->dim())) {
+      ExEnv::errn() << indent << "DistDiagSCMatrix::accumulate(SCMatrix*a): "
+           << "dimensions don't match\n";
+      abort();
+    }
+
+  SCMatrixBlockListIter i1, i2;
+  for (i1=la->blocklist->begin(),i2=blocklist->begin();
+       i1!=la->blocklist->end() && i2!=blocklist->end();
+       i1++,i2++) {
+      int n = i1.block()->ndat();
+      if (n != i2.block()->ndat()) {
+          ExEnv::errn() << indent << "DistDiagSCMatrix::accumulate "
+               << "mismatch: internal error" << endl;
+          abort();
+        }
+      double *dat1 = i1.block()->dat();
+      double *dat2 = i2.block()->dat();
+      for (int i=0; i<n; i++) {
+          dat2[i] += dat1[i];
+        }
+    }
+}
+
+double
+DistDiagSCMatrix::invert_this()
+{
+  Ref<SCMatrixSubblockIter> I = local_blocks(SCMatrixSubblockIter::Read);
+  double det = 1.0;
+  for (I->begin(); I->ready(); I->next()) {
+      int n = I->block()->ndat();
+      double *data = I->block()->dat();
+      for (int i=0; i<n; i++) {
+          det *= data[i];
+          data[i] = 1.0/data[i];
+        }
+    }
+  GrpProductReduce<double> gred;
+  messagegrp()->reduce(&det, 1, gred);
+  return det;
+}
+
+double
+DistDiagSCMatrix::determ_this()
+{
+  Ref<SCMatrixSubblockIter> I = local_blocks(SCMatrixSubblockIter::Read);
+  double det = 1.0;
+  for (I->begin(); I->ready(); I->next()) {
+      int n = I->block()->ndat();
+      double *data = I->block()->dat();
+      for (int i=0; i<n; i++) {
+          det *= data[i];
+        }
+    }
+  GrpProductReduce<double> gred;
+  messagegrp()->reduce(det, gred);
+  return det;
+}
+
+double
+DistDiagSCMatrix::trace()
+{
+  double ret=0.0;
+  Ref<SCMatrixSubblockIter> I = local_blocks(SCMatrixSubblockIter::Read);
+  for (I->begin(); I->ready(); I->next()) {
+      int n = I->block()->ndat();
+      double *data = I->block()->dat();
+      for (int i=0; i<n; i++) {
+          ret += data[i];
+        }
+    }
+  messagegrp()->sum(ret);
+  return ret;
+}
+
+void
+DistDiagSCMatrix::gen_invert_this()
+{
+  Ref<SCMatrixSubblockIter> I = local_blocks(SCMatrixSubblockIter::Read);
+  for (I->begin(); I->ready(); I->next()) {
+      int n = I->block()->ndat();
+      double *data = I->block()->dat();
+      for (int i=0; i<n; i++) {
+          if (fabs(data[i]) > 1.0e-8)
+              data[i] = 1.0/data[i];
+          else
+              data[i] = 0.0;
+        }
+    }
+}
+
+void
+DistDiagSCMatrix::element_op(const Ref<SCElementOp>& op)
+{
+  SCMatrixBlockListIter i;
+  for (i = blocklist->begin(); i != blocklist->end(); i++) {
+      op->process_base(i.block());
+    }
+  if (op->has_collect()) op->collect(messagegrp());
+}
+
+void
+DistDiagSCMatrix::element_op(const Ref<SCElementOp2>& op,
+                              DiagSCMatrix* m)
+{
+  DistDiagSCMatrix *lm
+      = require_dynamic_cast<DistDiagSCMatrix*>(m,"DistDiagSCMatrix::element_op");
+
+  if (!dim()->equiv(lm->dim())) {
+      ExEnv::errn() << indent << "DistDiagSCMatrix: bad element_op\n";
+      abort();
+    }
+  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_collect()) op->collect(messagegrp());
+}
+
+void
+DistDiagSCMatrix::element_op(const Ref<SCElementOp3>& op,
+                              DiagSCMatrix* m,DiagSCMatrix* n)
+{
+  DistDiagSCMatrix *lm
+      = require_dynamic_cast<DistDiagSCMatrix*>(m,"DistDiagSCMatrix::element_op");
+  DistDiagSCMatrix *ln
+      = require_dynamic_cast<DistDiagSCMatrix*>(n,"DistDiagSCMatrix::element_op");
+
+  if (!dim()->equiv(lm->dim()) || !dim()->equiv(ln->dim())) {
+      ExEnv::errn() << indent << "DistDiagSCMatrix: bad element_op\n";
+      abort();
+    }
+  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_collect()) op->collect(messagegrp());
+}
+
+Ref<SCMatrixSubblockIter>
+DistDiagSCMatrix::local_blocks(SCMatrixSubblockIter::Access access)
+{
+  return new SCMatrixListSubblockIter(access, blocklist);
+}
+
+Ref<SCMatrixSubblockIter>
+DistDiagSCMatrix::all_blocks(SCMatrixSubblockIter::Access access)
+{
+  return new DistSCMatrixListSubblockIter(access, blocklist, messagegrp());
+}
+
+void
+DistDiagSCMatrix::error(const char *msg)
+{
+  ExEnv::errn() << indent << "DistDiagSCMatrix: error: " << msg << endl;
+}
+
+Ref<DistSCMatrixKit>
+DistDiagSCMatrix::skit()
+{
+  return dynamic_cast<DistSCMatrixKit*>(kit().pointer());
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
