Permalink
Browse files

tensor cleaned up but still with shallow copy

  • Loading branch information...
robertjharrison committed Dec 10, 2009
1 parent 6b6056d commit 16ba1bf93757109b20b2e421508c28c5e3da204e
Showing with 4,358 additions and 3,556 deletions.
  1. +62 −51 aclocal.m4
  2. +4 −3 src/apps/examples/csqrt.cc
  3. +18 −18 src/apps/examples/h2.cc
  4. +1 −0 src/apps/examples/hatom_energy.cc
  5. +1 −1 src/apps/examples/he.cc
  6. +4 −3 src/apps/examples/interior_dirichlet.cc
  7. +91 −112 src/apps/examples/tdse1d.cc
  8. +5 −7 src/apps/examples/test_gmres.cc
  9. +1 −1 src/apps/examples/vnucso.cc
  10. +51 −18 src/apps/moldft/input
  11. +28 −25 src/apps/moldft/moldft.cc
  12. +2 −2 src/apps/moldft/molecularbasis.h
  13. +9 −9 src/apps/tdse/input
  14. +48 −27 src/lib/linalg/lapack.cc
  15. +4 −4 src/lib/linalg/solvers.cc
  16. +1 −1 src/lib/linalg/solvers.h
  17. +4 −4 src/lib/linalg/test_solvers.cc
  18. +4 −16 src/lib/mra/convolution1d.h
  19. +2 −2 src/lib/mra/funcdefaults.h
  20. +38 −36 src/lib/mra/funcimpl.h
  21. +4 −4 src/lib/mra/gfit.cc
  22. +22 −20 src/lib/mra/mraimpl.h
  23. +4 −4 src/lib/mra/mraplot.cc
  24. +6 −5 src/lib/mra/operator.h
  25. +2 −2 src/lib/mra/systolic.cc
  26. +3 −1 src/lib/mra/testbsh.cc
  27. +1 −1 src/lib/mra/testperiodic.cc
  28. +15 −15 src/lib/mra/testsuite.cc
  29. +1 −0 src/lib/mra/twoscale.cc
  30. +5 −5 src/lib/mra/vmra.h
  31. +4 −1 src/lib/tensor/Makefile.am
  32. +81 −169 src/lib/tensor/basetensor.cc
  33. +94 −43 src/lib/tensor/basetensor.h
  34. +978 −0 src/lib/tensor/graveyard
  35. +23 −23 src/lib/tensor/jimkernel.cc
  36. +0 −2 src/lib/tensor/mtxmq.cc
  37. +896 −0 src/lib/tensor/oldtest.cc
  38. +10 −1,099 src/lib/tensor/tensor.cc
  39. +1,181 −531 src/lib/tensor/tensor.h
  40. +151 −147 src/lib/tensor/tensor_macros.h
  41. +2 −2 src/lib/tensor/tensorexcept.h
  42. +1 −268 src/lib/tensor/tensoriter.cc
  43. +279 −0 src/lib/tensor/tensoriter.h
  44. +184 −870 src/lib/tensor/test.cc
  45. +10 −1 src/lib/world/Makefile.am
  46. +3 −0 src/lib/world/archive.h
  47. +1 −0 src/lib/world/array.h
  48. +19 −3 src/lib/world/typestuff.h
View
@@ -1,4 +1,4 @@
# generated automatically by aclocal 1.10.1 -*- Autoconf -*-
# generated automatically by aclocal 1.10.2 -*- Autoconf -*-
# Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004,
# 2005, 2006, 2007, 2008 Free Software Foundation, Inc.
@@ -13,13 +13,13 @@
m4_ifndef([AC_AUTOCONF_VERSION],
[m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl
m4_if(AC_AUTOCONF_VERSION, [2.61],,
[m4_warning([this file was generated for autoconf 2.61.
m4_if(m4_defn([AC_AUTOCONF_VERSION]), [2.63],,
[m4_warning([this file was generated for autoconf 2.63.
You have another version of autoconf. It may work, but is not guaranteed to.
If you have problems, you may need to regenerate the build system entirely.
To do so, use the procedure documented by the package, typically `autoreconf'.])])
# Copyright (C) 2002, 2003, 2005, 2006, 2007 Free Software Foundation, Inc.
# Copyright (C) 2002, 2003, 2005, 2006, 2007, 2008 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
@@ -34,7 +34,7 @@ AC_DEFUN([AM_AUTOMAKE_VERSION],
[am__api_version='1.10'
dnl Some users find AM_AUTOMAKE_VERSION and mistake it for a way to
dnl require some minimum version. Point them to the right macro.
m4_if([$1], [1.10.1], [],
m4_if([$1], [1.10.2], [],
[AC_FATAL([Do not call $0, use AM_INIT_AUTOMAKE([$1]).])])dnl
])
@@ -48,12 +48,12 @@ m4_define([_AM_AUTOCONF_VERSION], [])
# AM_SET_CURRENT_AUTOMAKE_VERSION
# -------------------------------
# Call AM_AUTOMAKE_VERSION and AM_AUTOMAKE_VERSION so they can be traced.
# This function is AC_REQUIREd by AC_INIT_AUTOMAKE.
# This function is AC_REQUIREd by AM_INIT_AUTOMAKE.
AC_DEFUN([AM_SET_CURRENT_AUTOMAKE_VERSION],
[AM_AUTOMAKE_VERSION([1.10.1])dnl
[AM_AUTOMAKE_VERSION([1.10.2])dnl
m4_ifndef([AC_AUTOCONF_VERSION],
[m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl
_AM_AUTOCONF_VERSION(AC_AUTOCONF_VERSION)])
_AM_AUTOCONF_VERSION(m4_defn([AC_AUTOCONF_VERSION]))])
# Figure out how to run the assembler. -*- Autoconf -*-
@@ -325,57 +325,68 @@ _AM_SUBST_NOTMAKE([AMDEPBACKSLASH])dnl
# Generate code to set up dependency tracking. -*- Autoconf -*-
# Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005
# Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2008
# Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.
#serial 3
#serial 5
# _AM_OUTPUT_DEPENDENCY_COMMANDS
# ------------------------------
AC_DEFUN([_AM_OUTPUT_DEPENDENCY_COMMANDS],
[for mf in $CONFIG_FILES; do
# Strip MF so we end up with the name of the file.
mf=`echo "$mf" | sed -e 's/:.*$//'`
# Check whether this is an Automake generated Makefile or not.
# We used to match only the files named `Makefile.in', but
# some people rename them; so instead we look at the file content.
# Grep'ing the first line is not enough: some people post-process
# each Makefile.in and add a new line on top of each file to say so.
# Grep'ing the whole file is not good either: AIX grep has a line
# limit of 2048, but all sed's we know have understand at least 4000.
if sed -n 's,^#.*generated by automake.*,X,p' "$mf" | grep X >/dev/null 2>&1; then
dirpart=`AS_DIRNAME("$mf")`
else
continue
fi
# Extract the definition of DEPDIR, am__include, and am__quote
# from the Makefile without running `make'.
DEPDIR=`sed -n 's/^DEPDIR = //p' < "$mf"`
test -z "$DEPDIR" && continue
am__include=`sed -n 's/^am__include = //p' < "$mf"`
test -z "am__include" && continue
am__quote=`sed -n 's/^am__quote = //p' < "$mf"`
# When using ansi2knr, U may be empty or an underscore; expand it
U=`sed -n 's/^U = //p' < "$mf"`
# Find all dependency output files, they are included files with
# $(DEPDIR) in their names. We invoke sed twice because it is the
# simplest approach to changing $(DEPDIR) to its actual value in the
# expansion.
for file in `sed -n "
s/^$am__include $am__quote\(.*(DEPDIR).*\)$am__quote"'$/\1/p' <"$mf" | \
sed -e 's/\$(DEPDIR)/'"$DEPDIR"'/g' -e 's/\$U/'"$U"'/g'`; do
# Make sure the directory exists.
test -f "$dirpart/$file" && continue
fdir=`AS_DIRNAME(["$file"])`
AS_MKDIR_P([$dirpart/$fdir])
# echo "creating $dirpart/$file"
echo '# dummy' > "$dirpart/$file"
[{
# Autoconf 2.62 quotes --file arguments for eval, but not when files
# are listed without --file. Let's play safe and only enable the eval
# if we detect the quoting.
case $CONFIG_FILES in
*\'*) eval set x "$CONFIG_FILES" ;;
*) set x $CONFIG_FILES ;;
esac
shift
for mf
do
# Strip MF so we end up with the name of the file.
mf=`echo "$mf" | sed -e 's/:.*$//'`
# Check whether this is an Automake generated Makefile or not.
# We used to match only the files named `Makefile.in', but
# some people rename them; so instead we look at the file content.
# Grep'ing the first line is not enough: some people post-process
# each Makefile.in and add a new line on top of each file to say so.
# Grep'ing the whole file is not good either: AIX grep has a line
# limit of 2048, but all sed's we know have understand at least 4000.
if sed -n 's,^#.*generated by automake.*,X,p' "$mf" | grep X >/dev/null 2>&1; then
dirpart=`AS_DIRNAME("$mf")`
else
continue
fi
# Extract the definition of DEPDIR, am__include, and am__quote
# from the Makefile without running `make'.
DEPDIR=`sed -n 's/^DEPDIR = //p' < "$mf"`
test -z "$DEPDIR" && continue
am__include=`sed -n 's/^am__include = //p' < "$mf"`
test -z "am__include" && continue
am__quote=`sed -n 's/^am__quote = //p' < "$mf"`
# When using ansi2knr, U may be empty or an underscore; expand it
U=`sed -n 's/^U = //p' < "$mf"`
# Find all dependency output files, they are included files with
# $(DEPDIR) in their names. We invoke sed twice because it is the
# simplest approach to changing $(DEPDIR) to its actual value in the
# expansion.
for file in `sed -n "
s/^$am__include $am__quote\(.*(DEPDIR).*\)$am__quote"'$/\1/p' <"$mf" | \
sed -e 's/\$(DEPDIR)/'"$DEPDIR"'/g' -e 's/\$U/'"$U"'/g'`; do
# Make sure the directory exists.
test -f "$dirpart/$file" && continue
fdir=`AS_DIRNAME(["$file"])`
AS_MKDIR_P([$dirpart/$fdir])
# echo "creating $dirpart/$file"
echo '# dummy' > "$dirpart/$file"
done
done
done
}
])# _AM_OUTPUT_DEPENDENCY_COMMANDS
@@ -669,13 +680,13 @@ esac
# Helper functions for option handling. -*- Autoconf -*-
# Copyright (C) 2001, 2002, 2003, 2005 Free Software Foundation, Inc.
# Copyright (C) 2001, 2002, 2003, 2005, 2008 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.
# serial 3
# serial 4
# _AM_MANGLE_OPTION(NAME)
# -----------------------
@@ -692,7 +703,7 @@ AC_DEFUN([_AM_SET_OPTION],
# ----------------------------------
# OPTIONS is a space-separated list of Automake options.
AC_DEFUN([_AM_SET_OPTIONS],
[AC_FOREACH([_AM_Option], [$1], [_AM_SET_OPTION(_AM_Option)])])
[m4_foreach_w([_AM_Option], [$1], [_AM_SET_OPTION(_AM_Option)])])
# _AM_IF_OPTION(OPTION, IF-SET, [IF-NOT-SET])
# -------------------------------------------
@@ -24,8 +24,9 @@ static std::complex<double> f(const coordT& r)
return std::complex<double>(0.0,2.0);
}
inline static void csqrt_op(const Key<1>& key, Tensor< std::complex<double> >& t) {
UNARY_OPTIMIZED_ITERATOR(std::complex<double>, t, *_p0 = sqrt(*_p0));
template <typename T>
inline static void csqrt_op(const Key<1>& key, Tensor< T >& t) {
UNARY_OPTIMIZED_ITERATOR(T, t, *_p0 = sqrt(*_p0));
}
@@ -49,7 +50,7 @@ int main(int argc, char** argv)
boring.truncate();
cfunctionT sqrt_of_boring = copy(boring);
sqrt_of_boring.unaryop(&csqrt_op);
sqrt_of_boring.unaryop(&csqrt_op<double_complex>);
double error = (square(sqrt_of_boring) - boring).norm2();
print("error is ", error);
world.gop.fence();
View
@@ -15,28 +15,27 @@ typedef FunctionFactory<double,3> factoryT;
typedef SeparatedConvolution<double,3> operatorT;
static const double R = 1.4; // bond length
static const double L = 32.0*R; // box size
static const long k = 5; // wavelet order
static const double thresh = 1e-3; // precision
static const double thresh1 = thresh*0.1;
static const double L = 64.0*R; // box size
static const long k = 6; // wavelet order
static const double thresh = 1e-4; // precision
static double guess(const coordT& r) {
const double x=r[0], y=r[1], z=r[2];
return (exp(-sqrt(x*x+y*y+(z-R/2)*(z-R/2)))+
exp(-sqrt(x*x+y*y+(z+R/2)*(z+R/2))));
return (exp(-sqrt(x*x+y*y+(z-R/2)*(z-R/2)+1e-8))+
exp(-sqrt(x*x+y*y+(z+R/2)*(z+R/2)+1e-8)));
}
static double V(const coordT& r) {
const double x=r[0], y=r[1], z=r[2];
return -1.0/sqrt(x*x+y*y+(z-R/2)*(z-R/2))+
-1.0/sqrt(x*x+y*y+(z+R/2)*(z+R/2));
return -1.0/sqrt(x*x+y*y+(z-R/2)*(z-R/2)+1e-8)+
-1.0/sqrt(x*x+y*y+(z+R/2)*(z+R/2)+1e-8);
}
void iterate(World& world, functionT& V, functionT& psi, double& eps) {
operatorT op = BSHOperator3D<double>(world, sqrt(-2*eps), k, 0.001, 1e-6);
functionT Vpsi = (V*psi);
Vpsi.scale(-2.0).truncate(thresh1);
functionT tmp = apply(op,Vpsi).truncate(thresh1);
Vpsi.scale(-2.0).truncate();
functionT tmp = apply(op,Vpsi).truncate();
double norm = tmp.norm2();
functionT r = tmp-psi;
double rnorm = r.norm2();
@@ -61,8 +60,8 @@ int main(int argc, char** argv) {
FunctionDefaults<3>::set_k(k);
FunctionDefaults<3>::set_thresh(thresh);
FunctionDefaults<3>::set_refine(true);
FunctionDefaults<3>::set_initial_level(2);
FunctionDefaults<3>::set_truncate_mode(0);
FunctionDefaults<3>::set_initial_level(5);
FunctionDefaults<3>::set_truncate_mode(1);
FunctionDefaults<3>::set_cubic_cell(-L/2, L/2);
// for (int i=0; i<3; i++) {
// FunctionDefaults<3>::cell(i,0) = -L/2;
@@ -71,15 +70,15 @@ int main(int argc, char** argv) {
functionT Vnuc = factoryT(world).f(V);
functionT psi = factoryT(world).f(guess);
psi.truncate(thresh1);
psi.truncate();
psi.scale(1.0/psi.norm2());
operatorT op = CoulombOperator<double>(world, k, 0.001, 1e-6);
double eps = -0.6;
for (int iter=0; iter<10; iter++) {
functionT rho = square(psi).truncate(thresh1);
functionT potential = Vnuc + apply(op,rho).truncate(thresh1);
functionT rho = square(psi).truncate();
functionT potential = Vnuc + apply(op,rho).truncate();
iterate(world, potential, psi, eps);
}
@@ -90,20 +89,21 @@ int main(int argc, char** argv) {
kinetic_energy += inner(dpsi,dpsi);
}
functionT rho = square(psi).truncate(thresh1);
functionT rho = square(psi);
double two_electron_energy = inner(apply(op,rho),rho);
double nuclear_attraction_energy = 2.0*inner(Vnuc*psi,psi);
double nuclear_attraction_energy = 2.0*inner(Vnuc,rho);
double nuclear_repulsion_energy = 1.0/R;
double total_energy = kinetic_energy + two_electron_energy +
nuclear_attraction_energy + nuclear_repulsion_energy;
double virial = (nuclear_attraction_energy + two_electron_energy + nuclear_repulsion_energy) / kinetic_energy;
if (world.rank() == 0) {
print(" Kinetic energy ", kinetic_energy);
print(" Nuclear attraction energy ", nuclear_attraction_energy);
print(" Two-electron energy ", two_electron_energy);
print(" Nuclear repulsion energy ", nuclear_repulsion_energy);
print(" Total energy ", total_energy);
print(" Virial ", (nuclear_attraction_energy + two_electron_energy) / kinetic_energy);
print(" Virial ", virial);
}
world.gop.fence();
@@ -22,6 +22,7 @@ int main(int argc, char**argv) {
// Load info for MADNESS numerical routines
startup(world,argc,argv);
std::cout.precision(8);
// Use defaults for numerical functions except for user simulation volume
FunctionDefaults<3>::set_cubic_cell(-20,20);
View
@@ -15,7 +15,7 @@ typedef FunctionFactory<double,3> factoryT;
typedef SeparatedConvolution<double,3> operatorT;
static const double L = 32.0; // box size
static const long k = 6; // wavelet order
static const long k = 8; // wavelet order
static const double thresh = 1e-6; // precision
static const double thresh1 = thresh;
@@ -86,8 +86,9 @@ static double exact_sol(const coordT3d & pt) {
// gives the B(phi(x)) function == phi^2 (1-phi)^2
// this is a unary op, given the phi function, computes b pointwise
inline static void b_phi(const Key<3> &key, Tensor<double> & t) {
UNARY_OPTIMIZED_ITERATOR(double, t,
template <typename T>
inline static void b_phi(const Key<3> &key, Tensor<T> & t) {
UNARY_OPTIMIZED_ITERATOR(T, t,
*_p0 = (*_p0) * (*_p0) * (1.0-(*_p0)) * (1.0-(*_p0)));
}
@@ -180,7 +181,7 @@ int main(int argc, char **argv) {
functionT b = copy(phi);
phi.truncate();
b.unaryop(&b_phi);
b.unaryop(&b_phi<double>);
// apply the scale factor to b
// scale is (36 / epsilon), which normalizes b, -1 from the Green's function
Oops, something went wrong.

0 comments on commit 16ba1bf

Please sign in to comment.