Index: ThirdParty/mpqc_open/src/lib/math/scmat/localvect.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/scmat/localvect.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/scmat/localvect.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,353 @@
+//
+// localvect.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 <math.h>
+
+#include <util/misc/formio.h>
+#include <util/keyval/keyval.h>
+#include <math/scmat/local.h>
+#include <math/scmat/cmatrix.h>
+#include <math/scmat/elemop.h>
+
+using namespace std;
+using namespace sc;
+
+/////////////////////////////////////////////////////////////////////////////
+// LocalSCVector member functions
+
+static ClassDesc LocalSCVector_cd(
+  typeid(LocalSCVector),"LocalSCVector",1,"public SCVector",
+  0, 0, 0);
+
+LocalSCVector::LocalSCVector(const RefSCDimension&a,
+                             LocalSCMatrixKit *kit):
+  SCVector(a,kit)
+{
+  resize(a->n());
+}
+
+void
+LocalSCVector::resize(int n)
+{
+  block = new SCVectorSimpleBlock(0,n);
+}
+
+LocalSCVector::~LocalSCVector()
+{
+}
+
+double *
+LocalSCVector::get_data()
+{
+  return block->data;
+}
+
+double
+LocalSCVector::get_element(int i) const
+{
+  int size = block->iend - block->istart;
+  if (i < 0 || i >= size) {
+      ExEnv::errn() << indent << "LocalSCVector::get_element: bad index\n";
+      abort();
+    }
+  return block->data[i];
+}
+
+void
+LocalSCVector::set_element(int i,double a)
+{
+  int size = block->iend - block->istart;
+  if (i < 0 || i >= size) {
+      ExEnv::errn() << indent << "LocalSCVector::set_element: bad index\n";
+      abort();
+    }
+  block->data[i] = a;
+}
+
+void
+LocalSCVector::accumulate_element(int i,double a)
+{
+  int size = block->iend - block->istart;
+  if (i < 0 || i >= size) {
+      ExEnv::errn() << indent << "LocalSCVector::accumulate_element: bad index\n";
+      abort();
+    }
+  block->data[i] += a;
+}
+
+void
+LocalSCVector::accumulate_product_rv(SCMatrix*a,SCVector*b)
+{
+  const char* name = "LocalSCVector::accumulate_product";
+  // make sure that the arguments are of the correct type
+  LocalSCMatrix* la = require_dynamic_cast<LocalSCMatrix*>(a,name);
+  LocalSCVector* lb = require_dynamic_cast<LocalSCVector*>(b,name);
+
+  // make sure that the dimensions match
+  if (!dim()->equiv(la->rowdim()) || !la->coldim()->equiv(lb->dim())) {
+      ExEnv::errn() << indent
+           << "LocalSCVector:: accumulate_product(SCMatrix*a,SCVector*b): "
+           << "dimensions don't match\n";
+      abort();
+    }
+
+  cmat_mxm(la->rows, 0,
+           &(lb->block->data), 1,
+           &(block->data), 1,
+           n(), la->ncol(), 1,
+           1);
+}
+
+void
+LocalSCVector::accumulate_product_sv(SymmSCMatrix*a,SCVector*b)
+{
+  const char* name = "LocalSCVector::accumulate_product";
+  // make sure that the arguments are of the correct type
+  LocalSymmSCMatrix* la = require_dynamic_cast<LocalSymmSCMatrix*>(a,name);
+  LocalSCVector* lb = require_dynamic_cast<LocalSCVector*>(b,name);
+
+  // make sure that the dimensions match
+  if (!dim()->equiv(la->dim()) || !la->dim()->equiv(lb->dim())) {
+      ExEnv::errn() << indent
+           << "LocalSCVector:: accumulate_product(SymmSCMatrix*a,SCVector*b): "
+           << "dimensions don't match\n";
+      abort();
+    }
+
+  double* thisdat = block->data;
+  double** adat = la->rows;
+  double* bdat = lb->block->data;
+  double tmp;
+  int n = dim()->n();
+  int i, j;
+  for (i=0; i<n; i++) {
+      tmp = 0.0;
+      for (j=0; j<=i; j++) {
+          tmp += adat[i][j] * bdat[j];
+        }
+      for (; j<n; j++) {
+          tmp += adat[j][i] * bdat[j];
+        }
+      thisdat[i] += tmp;
+    }
+}
+
+void
+LocalSCVector::accumulate(const SCVector*a)
+{
+  // make sure that the argument is of the correct type
+  const LocalSCVector* la
+    = require_dynamic_cast<const LocalSCVector*>(a,"LocalSCVector::accumulate");
+
+  // make sure that the dimensions match
+  if (!dim()->equiv(la->dim())) {
+      ExEnv::errn() << indent << "LocalSCVector::accumulate(SCVector*a): "
+           << "dimensions don't match\n";
+      abort();
+    }
+
+  int nelem = d->n();
+  int i;
+  for (i=0; i<nelem; i++) block->data[i] += la->block->data[i];
+}
+
+void
+LocalSCVector::accumulate(const SCMatrix*a)
+{
+  // make sure that the argument is of the correct type
+  const LocalSCMatrix* la
+    = require_dynamic_cast<const LocalSCMatrix*>(a,"LocalSCVector::accumulate");
+
+  // make sure that the dimensions match
+  if (!((la->rowdim()->equiv(dim()) && la->coldim()->n() == 1)
+        || (la->coldim()->equiv(dim()) && la->rowdim()->n() == 1))) {
+      ExEnv::errn() << indent << "LocalSCVector::accumulate(SCMatrix*a): "
+           << "dimensions don't match\n";
+      abort();
+    }
+
+  int nelem = d->n();
+  int i;
+  for (i=0; i<nelem; i++) block->data[i] += la->block->data[i];
+}
+
+void
+LocalSCVector::assign_val(double a)
+{
+  int nelem = d->n();
+  int i;
+  for (i=0; i<nelem; i++) block->data[i] = a;
+}
+
+void
+LocalSCVector::assign_v(SCVector*a)
+{
+  // make sure that the argument is of the correct type
+  LocalSCVector* la
+    = require_dynamic_cast<LocalSCVector*>(a,"LocalSCVector::assign_v");
+
+  // make sure that the dimensions match
+  if (!dim()->equiv(la->dim())) {
+      ExEnv::errn() << indent << "LocalSCVector::assign_v(SCVector*a): "
+           << "dimensions don't match\n";
+      abort();
+    }
+
+  int nelem = d->n();
+  int i;
+  for (i=0; i<nelem; i++) block->data[i] = la->block->data[i];
+}
+
+void
+LocalSCVector::assign_p(const double*a)
+{
+  int nelem = d->n();
+  int i;
+  for (i=0; i<nelem; i++) block->data[i] = a[i];
+}
+
+double
+LocalSCVector::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 = d->n();
+  int i;
+  double result = 0.0;
+  for (i=0; i<nelem; i++) result += block->data[i] * la->block->data[i];
+  return result;
+}
+
+void
+LocalSCVector::element_op(const Ref<SCElementOp>& op)
+{
+  op->process_spec_vsimp(block.pointer());
+}
+
+void
+LocalSCVector::element_op(const Ref<SCElementOp2>& op,
+                          SCVector* m)
+{
+  LocalSCVector *lm
+      = require_dynamic_cast<LocalSCVector*>(m, "LocalSCVector::element_op");
+
+  if (!dim()->equiv(lm->dim())) {
+      ExEnv::errn() << indent << "LocalSCVector: bad element_op\n";
+      abort();
+    }
+  op->process_spec_vsimp(block.pointer(), lm->block.pointer());
+}
+
+void
+LocalSCVector::element_op(const Ref<SCElementOp3>& op,
+                          SCVector* m,SCVector* n)
+{
+  LocalSCVector *lm
+      = require_dynamic_cast<LocalSCVector*>(m, "LocalSCVector::element_op");
+  LocalSCVector *ln
+      = require_dynamic_cast<LocalSCVector*>(n, "LocalSCVector::element_op");
+
+  if (!dim()->equiv(lm->dim()) || !dim()->equiv(ln->dim())) {
+      ExEnv::errn() << indent << "LocalSCVector: bad element_op\n";
+      abort();
+    }
+  op->process_spec_vsimp(block.pointer(),
+                         lm->block.pointer(), ln->block.pointer());
+}
+
+// from Ed Seidl at the NIH (with a bit of hacking)
+void
+LocalSCVector::vprint(const char *title, ostream& os, int prec) const
+{
+  int i;
+  int lwidth;
+  double max=this->maxabs();
+
+  max = (max==0.0) ? 1.0 : log10(max);
+  if (max < 0.0) max=1.0;
+
+  lwidth = prec + 5 + (int) max;
+
+  if (title)
+    os << endl << indent << title << endl;
+  else
+    os << endl;
+
+  if (n()==0) {
+    os << indent << "empty vector\n";
+    return;
+  }
+
+  for (i=0; i<n(); i++)
+    os << indent
+       << scprintf("%5d %*.*f\n",i+1,lwidth,prec,block->data[i]);
+  os << endl;
+
+  os.flush();
+}
+
+Ref<SCMatrixSubblockIter>
+LocalSCVector::local_blocks(SCMatrixSubblockIter::Access access)
+{
+  if (messagegrp()->n() > 1) {
+      ExEnv::errn() << indent
+           << "LocalSCVector::local_blocks: not valid for local matrices"
+           << endl;
+      abort();
+    }
+  Ref<SCMatrixSubblockIter> iter
+      = new SCMatrixSimpleSubblockIter(access, block.pointer());
+  return iter;
+}
+
+Ref<SCMatrixSubblockIter>
+LocalSCVector::all_blocks(SCMatrixSubblockIter::Access access)
+{
+  if (access == SCMatrixSubblockIter::Write) {
+      ExEnv::errn() << indent << "LocalVectorSCMatrix::all_blocks: "
+           << "Write access permitted for local blocks only"
+           << endl;
+      abort();
+    }
+  return local_blocks(access);
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
