Initial revision

darcs-hash:20010301114010-f6438-12fb8a9ffcc80e86c0a97e37b5b0dae0dbc59b79.gz

17 files changed, 2663 insertions, 0 deletions

diff --git a/CarpetDev/COPYING b/CarpetDev/COPYING
new file mode 100644
index 000000000..1942c4334
--- /dev/null
+++ b/CarpetDev/COPYING
@@ -0,0 +1,341 @@
+		    GNU GENERAL PUBLIC LICENSE
+		       Version 2, June 1991
+
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.
+                          59 Temple Place - Suite 330
+                          Boston, MA 02111-1307, USA.
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+			    Preamble
+
+  The licenses for most software are designed to take away your
+freedom to share and change it.  By contrast, the GNU General Public
+License is intended to guarantee your freedom to share and change free
+software--to make sure the software is free for all its users.  This
+General Public License applies to most of the Free Software
+Foundation's software and to any other program whose authors commit to
+using it.  (Some other Free Software Foundation software is covered by
+the GNU Library General Public License instead.)  You can apply it to
+your programs, too.
+
+  When we speak of free software, we are referring to freedom, not
+price.  Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+this service if you wish), that you receive source code or can get it
+if you want it, that you can change the software or use pieces of it
+in new free programs; and that you know you can do these things.
+
+  To protect your rights, we need to make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the rights.
+These restrictions translate to certain responsibilities for you if you
+distribute copies of the software, or if you modify it.
+
+  For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must give the recipients all the rights that
+you have.  You must make sure that they, too, receive or can get the
+source code.  And you must show them these terms so they know their
+rights.
+
+  We protect your rights with two steps: (1) copyright the software, and
+(2) offer you this license which gives you legal permission to copy,
+distribute and/or modify the software.
+
+  Also, for each author's protection and ours, we want to make certain
+that everyone understands that there is no warranty for this free
+software.  If the software is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original, so
+that any problems introduced by others will not reflect on the original
+authors' reputations.
+
+  Finally, any free program is threatened constantly by software
+patents.  We wish to avoid the danger that redistributors of a free
+program will individually obtain patent licenses, in effect making the
+program proprietary.  To prevent this, we have made it clear that any
+patent must be licensed for everyone's free use or not licensed at all.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.
+
+		    GNU GENERAL PUBLIC LICENSE
+   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+  0. This License applies to any program or other work which contains
+a notice placed by the copyright holder saying it may be distributed
+under the terms of this General Public License.  The "Program", below,
+refers to any such program or work, and a "work based on the Program"
+means either the Program or any derivative work under copyright law:
+that is to say, a work containing the Program or a portion of it,
+either verbatim or with modifications and/or translated into another
+language.  (Hereinafter, translation is included without limitation in
+the term "modification".)  Each licensee is addressed as "you".
+
+Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope.  The act of
+running the Program is not restricted, and the output from the Program
+is covered only if its contents constitute a work based on the
+Program (independent of having been made by running the Program).
+Whether that is true depends on what the Program does.
+
+  1. You may copy and distribute verbatim copies of the Program's
+source code as you receive it, in any medium, provided that you
+conspicuously and appropriately publish on each copy an appropriate
+copyright notice and disclaimer of warranty; keep intact all the
+notices that refer to this License and to the absence of any warranty;
+and give any other recipients of the Program a copy of this License
+along with the Program.
+
+You may charge a fee for the physical act of transferring a copy, and
+you may at your option offer warranty protection in exchange for a fee.
+
+  2. You may modify your copy or copies of the Program or any portion
+of it, thus forming a work based on the Program, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+    a) You must cause the modified files to carry prominent notices
+    stating that you changed the files and the date of any change.
+
+    b) You must cause any work that you distribute or publish, that in
+    whole or in part contains or is derived from the Program or any
+    part thereof, to be licensed as a whole at no charge to all third
+    parties under the terms of this License.
+
+    c) If the modified program normally reads commands interactively
+    when run, you must cause it, when started running for such
+    interactive use in the most ordinary way, to print or display an
+    announcement including an appropriate copyright notice and a
+    notice that there is no warranty (or else, saying that you provide
+    a warranty) and that users may redistribute the program under
+    these conditions, and telling the user how to view a copy of this
+    License.  (Exception: if the Program itself is interactive but
+    does not normally print such an announcement, your work based on
+    the Program is not required to print an announcement.)
+
+These requirements apply to the modified work as a whole.  If
+identifiable sections of that work are not derived from the Program,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works.  But when you
+distribute the same sections as part of a whole which is a work based
+on the Program, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Program.
+
+In addition, mere aggregation of another work not based on the Program
+with the Program (or with a work based on the Program) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+  3. You may copy and distribute the Program (or a work based on it,
+under Section 2) in object code or executable form under the terms of
+Sections 1 and 2 above provided that you also do one of the following:
+
+    a) Accompany it with the complete corresponding machine-readable
+    source code, which must be distributed under the terms of Sections
+    1 and 2 above on a medium customarily used for software interchange; or,
+
+    b) Accompany it with a written offer, valid for at least three
+    years, to give any third party, for a charge no more than your
+    cost of physically performing source distribution, a complete
+    machine-readable copy of the corresponding source code, to be
+    distributed under the terms of Sections 1 and 2 above on a medium
+    customarily used for software interchange; or,
+
+    c) Accompany it with the information you received as to the offer
+    to distribute corresponding source code.  (This alternative is
+    allowed only for noncommercial distribution and only if you
+    received the program in object code or executable form with such
+    an offer, in accord with Subsection b above.)
+
+The source code for a work means the preferred form of the work for
+making modifications to it.  For an executable work, complete source
+code means all the source code for all modules it contains, plus any
+associated interface definition files, plus the scripts used to
+control compilation and installation of the executable.  However, as a
+special exception, the source code distributed need not include
+anything that is normally distributed (in either source or binary
+form) with the major components (compiler, kernel, and so on) of the
+operating system on which the executable runs, unless that component
+itself accompanies the executable.
+
+If distribution of executable or object code is made by offering
+access to copy from a designated place, then offering equivalent
+access to copy the source code from the same place counts as
+distribution of the source code, even though third parties are not
+compelled to copy the source along with the object code.
+
+  4. You may not copy, modify, sublicense, or distribute the Program
+except as expressly provided under this License.  Any attempt
+otherwise to copy, modify, sublicense or distribute the Program is
+void, and will automatically terminate your rights under this License.
+However, parties who have received copies, or rights, from you under
+this License will not have their licenses terminated so long as such
+parties remain in full compliance.
+
+  5. You are not required to accept this License, since you have not
+signed it.  However, nothing else grants you permission to modify or
+distribute the Program or its derivative works.  These actions are
+prohibited by law if you do not accept this License.  Therefore, by
+modifying or distributing the Program (or any work based on the
+Program), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Program or works based on it.
+
+  6. Each time you redistribute the Program (or any work based on the
+Program), the recipient automatically receives a license from the
+original licensor to copy, distribute or modify the Program subject to
+these terms and conditions.  You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties to
+this License.
+
+  7. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Program at all.  For example, if a patent
+license would not permit royalty-free redistribution of the Program by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Program.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system, which is
+implemented by public license practices.  Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+
+  8. If the distribution and/or use of the Program is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Program under this License
+may add an explicit geographical distribution limitation excluding
+those countries, so that distribution is permitted only in or among
+countries not thus excluded.  In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+  9. The Free Software Foundation may publish revised and/or new versions
+of the General Public License from time to time.  Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+Each version is given a distinguishing version number.  If the Program
+specifies a version number of this License which applies to it and "any
+later version", you have the option of following the terms and conditions
+either of that version or of any later version published by the Free
+Software Foundation.  If the Program does not specify a version number of
+this License, you may choose any version ever published by the Free Software
+Foundation.
+
+  10. If you wish to incorporate parts of the Program into other free
+programs whose distribution conditions are different, write to the author
+to ask for permission.  For software which is copyrighted by the Free
+Software Foundation, write to the Free Software Foundation; we sometimes
+make exceptions for this.  Our decision will be guided by the two goals
+of preserving the free status of all derivatives of our free software and
+of promoting the sharing and reuse of software generally.
+
+			    NO WARRANTY
+
+  11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+  12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+
+		     END OF TERMS AND CONDITIONS
+
+	    How to Apply These Terms to Your New Programs
+
+  If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+  To do so, attach the following notices to the program.  It is safest
+to attach them to the start of each source file to most effectively
+convey the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+    <one line to give the program's name and a brief idea of what it does.>
+    Copyright (C) 19yy  <name of author>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program 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 General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program; see the file COPYING.  If not, write to
+    the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+    Boston, MA 02111-1307, USA.
+
+Also add information on how to contact you by electronic and paper mail.
+
+If the program is interactive, make it output a short notice like this
+when it starts in an interactive mode:
+
+    Gnomovision version 69, Copyright (C) 19yy name of author
+    Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+    This is free software, and you are welcome to redistribute it
+    under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License.  Of course, the commands you use may
+be called something other than `show w' and `show c'; they could even be
+mouse-clicks or menu items--whatever suits your program.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the program, if
+necessary.  Here is a sample; alter the names:
+
+  Yoyodyne, Inc., hereby disclaims all copyright interest in the program
+  `Gnomovision' (which makes passes at compilers) written by James Hacker.
+
+  <signature of Ty Coon>, 1 April 1989
+  Ty Coon, President of Vice
+
+This General Public License does not permit incorporating your program into
+proprietary programs.  If your program is a subroutine library, you may
+consider it more useful to permit linking proprietary applications with the
+library.  If this is what you want to do, use the GNU Library General
+Public License instead of this License.
diff --git a/CarpetDev/CarpetCG/README b/CarpetDev/CarpetCG/README
new file mode 100644
index 000000000..5c55ee24d
--- /dev/null
+++ b/CarpetDev/CarpetCG/README
@@ -0,0 +1,14 @@
+CVS info   : $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetCG/README,v 1.1.1.1 2003/11/19 10:40:47 hawke Exp $
+
+Cactus Code Thorn CarpetCG
+Thorn Author(s)     : Ian Hawke <hawke@aei.mpg.de>
+                    : Erik Schnetter <eschnett@aei.mpg.de>
+Thorn Maintainer(s) : Ian Hawke <hawke@aei.mpg.de>
+                    : Erik Schnetter <eschnett@aei.mpg.de>
+--------------------------------------------------------------------------
+
+Purpose of the thorn:
+
+A conjugate gradient elliptic solver for Carpet using the TAT
+interface. Essentially a straight modification of TATCG using
+TATJacobi as a template.
diff --git a/CarpetDev/CarpetCG/doc/documentation.tex b/CarpetDev/CarpetCG/doc/documentation.tex
new file mode 100644
index 000000000..8f45346ef
--- /dev/null
+++ b/CarpetDev/CarpetCG/doc/documentation.tex
@@ -0,0 +1,144 @@
+% *======================================================================*
+%  Cactus Thorn template for ThornGuide documentation
+%  Author: Ian Kelley
+%  Date: Sun Jun 02, 2002
+%  $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetCG/doc/documentation.tex,v 1.1.1.1 2003/11/19 10:40:47 hawke Exp $                                                             
+%
+%  Thorn documentation in the latex file doc/documentation.tex 
+%  will be included in ThornGuides built with the Cactus make system.
+%  The scripts employed by the make system automatically include 
+%  pages about variables, parameters and scheduling parsed from the 
+%  relevant thorn CCL files.
+%  
+%  This template contains guidelines which help to assure that your     
+%  documentation will be correctly added to ThornGuides. More 
+%  information is available in the Cactus UsersGuide.
+%                                                    
+%  Guidelines:
+%   - Do not change anything before the line
+%       % START CACTUS THORNGUIDE",
+%     except for filling in the title, author, date, etc. fields.
+%        - Each of these fields should only be on ONE line.
+%        - Author names should be separated with a \\ or a comma.
+%   - You can define your own macros, but they must appear after
+%     the START CACTUS THORNGUIDE line, and must not redefine standard 
+%     latex commands.
+%   - To avoid name clashes with other thorns, 'labels', 'citations', 
+%     'references', and 'image' names should conform to the following 
+%     convention:          
+%       ARRANGEMENT_THORN_LABEL
+%     For example, an image wave.eps in the arrangement CactusWave and 
+%     thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps
+%   - Graphics should only be included using the graphicx package. 
+%     More specifically, with the "\includegraphics" command.  Do
+%     not specify any graphic file extensions in your .tex file. This 
+%     will allow us to create a PDF version of the ThornGuide
+%     via pdflatex.
+%   - References should be included with the latex "\bibitem" command. 
+%   - Use \begin{abstract}...\end{abstract} instead of \abstract{...}
+%   - Do not use \appendix, instead include any appendices you need as 
+%     standard sections. 
+%   - For the benefit of our Perl scripts, and for future extensions, 
+%     please use simple latex.     
+%
+% *======================================================================* 
+% 
+% Example of including a graphic image:
+%    \begin{figure}[ht]
+% 	\begin{center}
+%    	   \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure}
+% 	\end{center}
+% 	\caption{Illustration of this and that}
+% 	\label{MyArrangement_MyThorn_MyLabel}
+%    \end{figure}
+%
+% Example of using a label:
+%   \label{MyArrangement_MyThorn_MyLabel}
+%
+% Example of a citation:
+%    \cite{MyArrangement_MyThorn_Author99}
+%
+% Example of including a reference
+%   \bibitem{MyArrangement_MyThorn_Author99}
+%   {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999), 
+%   1--16. {\tt http://www.nowhere.com/}}
+%
+% *======================================================================* 
+
+% If you are using CVS use this line to give version information
+% $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetCG/doc/documentation.tex,v 1.1.1.1 2003/11/19 10:40:47 hawke Exp $
+
+\documentclass{article}
+
+% Use the Cactus ThornGuide style file
+% (Automatically used from Cactus distribution, if you have a 
+%  thorn without the Cactus Flesh download this from the Cactus
+%  homepage at www.cactuscode.org)
+\usepackage{../../../../doc/latex/cactus}
+
+\begin{document}
+
+% The author of the documentation
+\author{Ian Hawke \textless hawke@aei.mpg.de\textgreater \\ Erik Schnetter \textless eschnett@aei.mpg.de\textgreater}
+
+% The title of the document (not necessarily the name of the Thorn)
+\title{}
+
+% the date your document was last changed, if your document is in CVS, 
+% please use:
+%    \date{$ $Date: 2003/11/19 10:40:47 $ $}
+\date{November 16 2003}
+
+\maketitle
+
+% Do not delete next line
+% START CACTUS THORNGUIDE
+
+% Add all definitions used in this documentation here 
+%   \def\mydef etc
+
+% Add an abstract for this thorn's documentation
+\begin{abstract}
+
+\end{abstract}
+
+% The following sections are suggestive only.
+% Remove them or add your own.
+
+\section{Introduction}
+
+\section{Physical System}
+
+\section{Numerical Implementation}
+
+\section{Using This Thorn}
+
+\subsection{Obtaining This Thorn}
+
+\subsection{Basic Usage}
+
+\subsection{Special Behaviour}
+
+\subsection{Interaction With Other Thorns}
+
+\subsection{Examples}
+
+\subsection{Support and Feedback}
+
+\section{History}
+
+\subsection{Thorn Source Code}
+
+\subsection{Thorn Documentation}
+
+\subsection{Acknowledgements}
+
+
+\begin{thebibliography}{9}
+
+\end{thebibliography}
+
+% Do not delete next line
+% END CACTUS THORNGUIDE
+
+\end{document}
diff --git a/CarpetDev/CarpetCG/interface.ccl b/CarpetDev/CarpetCG/interface.ccl
new file mode 100644
index 000000000..0f1eaa6e4
--- /dev/null
+++ b/CarpetDev/CarpetCG/interface.ccl
@@ -0,0 +1,12 @@
+# Interface definition for thorn CarpetCG
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetCG/interface.ccl,v 1.1.1.1 2003/11/19 10:40:47 hawke Exp $
+
+implements: CarpetCG
+inherits: TATelliptic
+
+USES INCLUDE HEADER: carpet.hh
+USES INCLUDE HEADER: TATelliptic.h
+
+real prsvars[cg_maxsolvevars] TYPE=GF Timelevels=1 tags='Prolongation="None"'
+
+real dirvars[cg_maxsolvevars] TYPE=GF Timelevels=1 tags='Prolongation="None"'
diff --git a/CarpetDev/CarpetCG/param.ccl b/CarpetDev/CarpetCG/param.ccl
new file mode 100644
index 000000000..93c0c4ff2
--- /dev/null
+++ b/CarpetDev/CarpetCG/param.ccl
@@ -0,0 +1,25 @@
+# Parameter definitions for thorn CarpetCG
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetCG/param.ccl,v 1.1.1.1 2003/11/19 10:40:47 hawke Exp $
+
+BOOLEAN verbose "Produce log output while running" STEERABLE=always
+{
+} "no"
+
+BOOLEAN veryverbose "Produce much log output while running" STEERABLE=always
+{
+} "no"
+
+CCTK_INT outevery "Produce output every n iterations" STEERABLE=always
+{
+  1:* :: ""
+} 1
+
+CCTK_INT outeveryseconds "Produce output every n seconds" STEERABLE=always
+{
+  0:* :: ""
+} 10
+
+CCTK_INT cg_maxsolvevars "Maximum number of variables to solve"
+{
+  0:* :: "Controls the temporary memory allocation"
+} 4
diff --git a/CarpetDev/CarpetCG/schedule.ccl b/CarpetDev/CarpetCG/schedule.ccl
new file mode 100644
index 000000000..dd5c9505b
--- /dev/null
+++ b/CarpetDev/CarpetCG/schedule.ccl
@@ -0,0 +1,10 @@
+# Schedule definitions for thorn CarpetCG
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetCG/schedule.ccl,v 1.1.1.1 2003/11/19 10:40:47 hawke Exp $
+
+STORAGE: prsvars
+STORAGE: dirvars
+
+SCHEDULE CarpetCG_register AT startup
+{
+  LANG: C
+} "Register the elliptic solver"
diff --git a/CarpetDev/CarpetCG/src/CG.cc b/CarpetDev/CarpetCG/src/CG.cc
new file mode 100644
index 000000000..e27e465be
--- /dev/null
+++ b/CarpetDev/CarpetCG/src/CG.cc
@@ -0,0 +1,1240 @@
+/* $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetCG/src/CG.cc,v 1.4 2004/01/25 14:57:28 schnetter Exp $ */
+
+#include <cassert>
+#include <cmath>
+#include <cstdio>
+#include <ctime>
+#include <iostream>
+#include <vector>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
+#include "carpet.hh"
+#include "TATelliptic.h"
+
+namespace Carpet {
+  // TODO: fix this
+  void Restrict (const cGH* cgh);
+};
+
+
+
+namespace CarpetCG {
+  
+  using namespace std;
+  using namespace Carpet;
+  
+    
+  extern "C" {
+    void CarpetCG_register ();
+  }
+  
+  int CarpetCG_solve (const cGH * const cctkGH,
+                      const int * const var,
+                      const int * const res,
+                      const int nvars,
+                      const int options_table,
+                      const calcfunc calcres,
+                      const calcfunc applybnds,
+                      void * const userdata);
+  
+  namespace common {
+    
+    const int * var;
+    const int * res;
+    int nvars;
+    int options_table;
+    calcfunc calcres;
+    calcfunc applybnds;
+    void * userdata;
+    
+    int ierr;
+
+    vector<int> nboundaryzones(2*dim);
+
+    CCTK_REAL factor;
+    vector<CCTK_REAL> factors;
+
+    vector<int> fromindex; // Can't pass things through CallLocalFunction.
+    vector<int> toindex;   // Instead assume everything is going from a GF to a GF.
+
+    CCTK_REAL realconstant; // With only one required constant
+    CCTK_REAL realoutput; // And one output
+    CCTK_REAL realoutput_count; // And one _more_ output
+    
+    void call_calcres (cGH * const cctkGH)
+    {
+      if (ierr) return;
+      ierr = calcres (cctkGH, options_table, userdata);
+    }
+    
+    void call_applybnds (cGH * const cctkGH)
+    {
+      if (ierr) return;
+      ierr = applybnds (cctkGH, options_table, userdata);
+    }
+
+    void
+    global_sum (cGH const * restrict const cctkGH,
+                CCTK_REAL * restrict const var)
+    {
+      static int initialised = 0;
+      static int sum_handle;
+      
+      CCTK_REAL local, global;
+      
+      int ierr;
+      
+      if (! initialised) {
+        initialised = 1;
+        sum_handle = CCTK_ReductionArrayHandle ("sum");
+        assert (sum_handle >= 0);
+      }
+      
+      local = *var;
+      ierr = CCTK_ReduceLocScalar
+        (cctkGH, -1, sum_handle, &local, &global, CCTK_VARIABLE_REAL);
+      assert (!ierr);
+      *var = global;
+    }
+
+    // Copy 
+    void call_copy (cGH * const cctkGH)
+    {
+      
+      cGroup groupdata;
+      ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[0]), &groupdata);
+      assert (!ierr);
+      cGroupDynamicData groupdyndata;
+      ierr = CCTK_GroupDynamicData (cctkGH, CCTK_GroupIndexFromVarI(var[0]),
+                                    &groupdyndata);
+      assert (!ierr);
+      const int thedim = groupdata.dim;
+      assert (thedim>=0 && thedim<=dim);
+      assert (thedim == groupdyndata.dim);
+      int lsh[dim], nghostzones[dim];
+      for (int d=0; d<thedim; ++d) {
+        lsh[d] = groupdyndata.lsh[d];
+        nghostzones[d] = groupdyndata.nghostzones[d];
+      }
+      for (int d=thedim; d<dim; ++d) {
+        lsh[d] = 1;
+        nghostzones[d] = 0;
+      }
+      for (int d=0; d<dim; ++d) {
+        assert (lsh[d]>=0);
+        assert (nghostzones[d] >= 0 && 2*nghostzones[d] <= lsh[d]);
+      }      
+
+      int const lsize = lsh[0] * lsh[1] * lsh[2];
+      int n;
+      double chksum=0;
+      double chksum2=0;
+      
+      for (n=0; n<nvars; ++n) {
+        CCTK_REAL * dst = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, toindex[n]);
+        CCTK_REAL * src = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, fromindex[n]);
+        //        cout << "Copy indices " << toindex[n] << " (" << CCTK_VarName(toindex[n]) << ") " << fromindex[n] << " (" << CCTK_VarName(fromindex[n]) << ") " << endl;
+        
+        memcpy (dst, src, lsize * sizeof (*dst));
+        for (int k=0; k<lsh[2];++k){
+          for (int j=0; j<lsh[2];++j){
+            for (int i=0; i<lsh[2];++i){
+              int ind = i + lsh[0] * (j + lsh[1] * k);
+              chksum+=dst[ind];
+              chksum2+=src[ind];
+            }
+          }
+        }
+        //        cout << "copied " << chksum << " " << chksum2 << endl;
+          
+      }
+    }
+
+    //Correct residual sign
+    void call_correct_residual_sign (cGH * const cctkGH)
+    {
+      
+      cGroup groupdata;
+      ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[0]), &groupdata);
+      assert (!ierr);
+      cGroupDynamicData groupdyndata;
+      ierr = CCTK_GroupDynamicData (cctkGH, CCTK_GroupIndexFromVarI(var[0]),
+                                    &groupdyndata);
+      assert (!ierr);
+      const int thedim = groupdata.dim;
+      assert (thedim>=0 && thedim<=dim);
+      assert (thedim == groupdyndata.dim);
+      int lsh[dim], nghostzones[dim], bbox[2*dim];
+      int lbnd[dim], ubnd[dim];
+      for (int d=0; d<thedim; ++d) {
+        lsh[d] = groupdyndata.lsh[d];
+        nghostzones[d] = groupdyndata.nghostzones[d];
+        bbox[2*d] = groupdyndata.bbox[2*d];
+        bbox[2*d+1] = groupdyndata.bbox[2*d+1];
+      }
+      for (int d=thedim; d<dim; ++d) {
+        lsh[d] = 1;
+        nghostzones[d] = 0;
+        bbox[2*d] = groupdyndata.bbox[2*d];
+        bbox[2*d+1] = groupdyndata.bbox[2*d+1];
+      }
+      for (int d=0; d<dim; ++d) {
+        assert (lsh[d]>=0);
+        assert (nghostzones[d] >= 0 && 2*nghostzones[d] <= lsh[d]);
+      }
+      for (int d=0; d<dim; ++d) {
+        lbnd[d] = 
+          (bbox[2*d  ] ? 
+           ( (nboundaryzones[2*d]>=0)   ? nboundaryzones[2*d]   : nghostzones[d])
+             : nghostzones[d]);
+        ubnd[d] = lsh[d] - 
+          (bbox[2*d+1] ? 
+           ( (nboundaryzones[2*d+1]>=0) ? nboundaryzones[2*d+1] : nghostzones[d])
+             : nghostzones[d]);
+      }
+
+      int const lsize = lsh[0] * lsh[1] * lsh[2];
+      int n;
+      int i, j, k;
+      int ind;
+
+      for (n=0; n<nvars; ++n) {
+        CCTK_REAL * src = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, fromindex[n]);
+        for (k=lbnd[2]; k<ubnd[2]; ++k) {
+          for (j=lbnd[1]; j<ubnd[1]; ++j) {
+            for (i=lbnd[0]; i<ubnd[0]; ++i) {
+              ind = i + lsh[0] * (j + lsh[1] * k);
+              src[ind] *= copysign(1.0, factors[n]);
+            }
+          }
+        }
+      }
+    }
+
+
+    // ``Precondition''. Ha ha ha.
+    void call_apply_preconditioner (cGH * const cctkGH)
+    {
+      
+      cGroup groupdata;
+      ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[0]), &groupdata);
+      assert (!ierr);
+      cGroupDynamicData groupdyndata;
+      ierr = CCTK_GroupDynamicData (cctkGH, CCTK_GroupIndexFromVarI(var[0]),
+                                    &groupdyndata);
+      assert (!ierr);
+      const int thedim = groupdata.dim;
+      assert (thedim>=0 && thedim<=dim);
+      assert (thedim == groupdyndata.dim);
+      int lsh[dim], nghostzones[dim], bbox[2*dim];
+      int lbnd[dim], ubnd[dim];
+      for (int d=0; d<thedim; ++d) {
+        lsh[d] = groupdyndata.lsh[d];
+        nghostzones[d] = groupdyndata.nghostzones[d];
+        bbox[2*d] = groupdyndata.bbox[2*d];
+        bbox[2*d+1] = groupdyndata.bbox[2*d+1];
+      }
+      for (int d=thedim; d<dim; ++d) {
+        lsh[d] = 1;
+        nghostzones[d] = 0;
+        bbox[2*d] = groupdyndata.bbox[2*d];
+        bbox[2*d+1] = groupdyndata.bbox[2*d+1];
+      }
+      for (int d=0; d<dim; ++d) {
+        assert (lsh[d]>=0);
+        assert (nghostzones[d] >= 0 && 2*nghostzones[d] <= lsh[d]);
+      }
+      for (int d=0; d<dim; ++d) {
+        lbnd[d] = 
+          (bbox[2*d  ] ? 
+           ( (nboundaryzones[2*d]>=0)   ? nboundaryzones[2*d]   : nghostzones[d])
+             : nghostzones[d]);
+        ubnd[d] = lsh[d] - 
+          (bbox[2*d+1] ? 
+           ( (nboundaryzones[2*d+1]>=0) ? nboundaryzones[2*d+1] : nghostzones[d])
+             : nghostzones[d]);
+      }
+
+      int const lsize = lsh[0] * lsh[1] * lsh[2];
+      int n;
+      int i, j, k;
+      int ind;
+
+      for (n=0; n<nvars; ++n) {
+        CCTK_REAL * dst = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, toindex[n]);
+        CCTK_REAL * src = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, fromindex[n]);
+        for (k=lbnd[2]; k<ubnd[2]; ++k) {
+          for (j=lbnd[1]; j<ubnd[1]; ++j) {
+            for (i=lbnd[0]; i<ubnd[0]; ++i) {
+              ind = i + lsh[0] * (j + lsh[1] * k);
+              dst[ind] = factor * fabs(factors[n]) * src[ind];
+            }
+          }
+        }
+      }
+    }
+
+
+    // Add scaled 
+    void call_add_scaled (cGH * const cctkGH)
+    {
+      
+      cGroup groupdata;
+      ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[0]), &groupdata);
+      assert (!ierr);
+      cGroupDynamicData groupdyndata;
+      ierr = CCTK_GroupDynamicData (cctkGH, CCTK_GroupIndexFromVarI(var[0]),
+                                    &groupdyndata);
+      assert (!ierr);
+      const int thedim = groupdata.dim;
+      assert (thedim>=0 && thedim<=dim);
+      assert (thedim == groupdyndata.dim);
+      int lsh[dim], nghostzones[dim], bbox[2*dim];
+      int lbnd[dim], ubnd[dim];
+      for (int d=0; d<thedim; ++d) {
+        lsh[d] = groupdyndata.lsh[d];
+        nghostzones[d] = groupdyndata.nghostzones[d];
+        bbox[2*d] = groupdyndata.bbox[2*d];
+        bbox[2*d+1] = groupdyndata.bbox[2*d+1];
+      }
+      for (int d=thedim; d<dim; ++d) {
+        lsh[d] = 1;
+        nghostzones[d] = 0;
+        bbox[2*d] = groupdyndata.bbox[2*d];
+        bbox[2*d+1] = groupdyndata.bbox[2*d+1];
+      }
+      for (int d=thedim; d<dim; ++d) {
+        lsh[d] = 1;
+        nghostzones[d] = 0;
+        bbox[2*d] = groupdyndata.bbox[2*d];
+        bbox[2*d+1] = groupdyndata.bbox[2*d+1];
+      }
+      for (int d=0; d<dim; ++d) {
+        assert (lsh[d]>=0);
+        assert (nghostzones[d] >= 0 && 2*nghostzones[d] <= lsh[d]);
+      }
+      for (int d=0; d<dim; ++d) {
+        lbnd[d] = 
+          (bbox[2*d  ] ? 
+           ( (nboundaryzones[2*d]>=0)   ? nboundaryzones[2*d]   : nghostzones[d])
+             : nghostzones[d]);
+        ubnd[d] = lsh[d] - 
+          (bbox[2*d+1] ? 
+           ( (nboundaryzones[2*d+1]>=0) ? nboundaryzones[2*d+1] : nghostzones[d])
+             : nghostzones[d]);
+      }
+      
+      int n;
+      int i, j, k;
+      int ind;
+
+      CCTK_REAL const alpha = realconstant;
+      double chksum=0;
+      double chksum2=0;
+      
+      for (n=0; n<nvars; ++n) {
+        CCTK_REAL * src = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, fromindex[n]);
+        CCTK_REAL * dst = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, toindex[n]);
+        //        cout << "Indices " << toindex[n] << " (" << CCTK_VarName(toindex[n]) << ") " << fromindex[n] << " (" << CCTK_VarName(fromindex[n]) << ") " << endl;
+        for (k=lbnd[2]; k<ubnd[2]; ++k) {
+          for (j=lbnd[1]; j<ubnd[1]; ++j) {
+            for (i=lbnd[0]; i<ubnd[0]; ++i) {
+              ind = i + lsh[0] * (j + lsh[1] * k);
+              dst[ind] += alpha * src[ind];
+              chksum+= dst[ind];
+              chksum2+= src[ind];
+              
+            }
+          }
+        }
+      }
+      //      cout << "add scaled " << alpha << " " << chksum << " " << chksum2 << endl;
+      
+    }
+    
+
+
+    // Scale and add 
+    void call_scale_and_add (cGH * const cctkGH)
+    {
+      
+      cGroup groupdata;
+      ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[0]), &groupdata);
+      assert (!ierr);
+      cGroupDynamicData groupdyndata;
+      ierr = CCTK_GroupDynamicData (cctkGH, CCTK_GroupIndexFromVarI(var[0]),
+                                    &groupdyndata);
+      assert (!ierr);
+      const int thedim = groupdata.dim;
+      assert (thedim>=0 && thedim<=dim);
+      assert (thedim == groupdyndata.dim);
+      int lsh[dim], nghostzones[dim], bbox[2*dim];
+      int lbnd[dim], ubnd[dim];
+      for (int d=0; d<thedim; ++d) {
+        lsh[d] = groupdyndata.lsh[d];
+        nghostzones[d] = groupdyndata.nghostzones[d];
+        bbox[2*d] = groupdyndata.bbox[2*d];
+        bbox[2*d+1] = groupdyndata.bbox[2*d+1];
+      }
+      for (int d=thedim; d<dim; ++d) {
+        lsh[d] = 1;
+        nghostzones[d] = 0;
+        bbox[2*d] = groupdyndata.bbox[2*d];
+        bbox[2*d+1] = groupdyndata.bbox[2*d+1];
+      }
+      for (int d=0; d<dim; ++d) {
+        assert (lsh[d]>=0);
+        assert (nghostzones[d] >= 0 && 2*nghostzones[d] <= lsh[d]);
+      }      
+      for (int d=0; d<dim; ++d) {
+        lbnd[d] = 
+          (bbox[2*d  ] ? 
+           ( (nboundaryzones[2*d]>=0)   ? nboundaryzones[2*d]   : nghostzones[d])
+             : nghostzones[d]);
+        ubnd[d] = lsh[d] - 
+          (bbox[2*d+1] ? 
+           ( (nboundaryzones[2*d+1]>=0) ? nboundaryzones[2*d+1] : nghostzones[d])
+             : nghostzones[d]);
+      }
+
+
+      int n;
+      int i, j, k;
+      int ind;
+
+      CCTK_REAL const alpha = realconstant;
+      
+      for (n=0; n<nvars; ++n) {
+        CCTK_REAL * src = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, fromindex[n]);
+        CCTK_REAL * dst = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, toindex[n]);
+        for (k=lbnd[2]; k<ubnd[2]; ++k) {
+          for (j=lbnd[1]; j<ubnd[1]; ++j) {
+            for (i=lbnd[0]; i<ubnd[0]; ++i) {
+              ind = i + lsh[0] * (j + lsh[1] * k);
+              dst[ind] = alpha * dst[ind] + src[ind];
+            }
+          }
+        }
+      }
+    }
+
+    // Dot product 
+    void call_dot_product (cGH * const cctkGH)
+    {
+      
+      cGroup groupdata;
+      ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[0]), &groupdata);
+      assert (!ierr);
+      cGroupDynamicData groupdyndata;
+      ierr = CCTK_GroupDynamicData (cctkGH, CCTK_GroupIndexFromVarI(var[0]),
+                                    &groupdyndata);
+      assert (!ierr);
+      const int thedim = groupdata.dim;
+      assert (thedim>=0 && thedim<=dim);
+      assert (thedim == groupdyndata.dim);
+      int lsh[dim], nghostzones[dim], bbox[2*dim];
+      int lbnd[dim], ubnd[dim];
+      for (int d=0; d<thedim; ++d) {
+        lsh[d] = groupdyndata.lsh[d];
+        nghostzones[d] = groupdyndata.nghostzones[d];
+        bbox[2*d] = groupdyndata.bbox[2*d];
+        bbox[2*d+1] = groupdyndata.bbox[2*d+1];
+      }
+      for (int d=thedim; d<dim; ++d) {
+        lsh[d] = 1;
+        nghostzones[d] = 0;
+        bbox[2*d] = groupdyndata.bbox[2*d];
+        bbox[2*d+1] = groupdyndata.bbox[2*d+1];
+      }
+      for (int d=0; d<dim; ++d) {
+        assert (lsh[d]>=0);
+        assert (nghostzones[d] >= 0 && 2*nghostzones[d] <= lsh[d]);
+      }      
+      for (int d=0; d<dim; ++d) {
+        lbnd[d] = 
+          (bbox[2*d  ] ? 
+           ( (nboundaryzones[2*d]>=0)   ? nboundaryzones[2*d]   : nghostzones[d])
+             : nghostzones[d]);
+        ubnd[d] = lsh[d] - 
+          (bbox[2*d+1] ? 
+           ( (nboundaryzones[2*d+1]>=0) ? nboundaryzones[2*d+1] : nghostzones[d])
+             : nghostzones[d]);
+      }
+
+      int n;
+      int i, j, k;
+      int ind;
+      CCTK_REAL res;
+
+      double chksum=0;
+      double chksum2=0;
+      
+      res = 0;
+      realoutput_count=0;
+      
+      // Shift the count to here instead - what a waste of space
+      for (k=0; k<lsh[2]; ++k) {
+        for (j=0; j<lsh[1]; ++j) {
+          for (i=0; i<lsh[0]; ++i) {
+            realoutput_count++;
+          }
+        }
+      }
+      
+      for (n=0; n<nvars; ++n) {
+        CCTK_REAL * a = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, toindex[n]);
+        CCTK_REAL * b = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, fromindex[n]);
+        //        cout << "Indices " << toindex[n] << " (" << CCTK_VarName(toindex[n]) << ") " << fromindex[n] << " (" << CCTK_VarName(fromindex[n]) << ") " << endl;
+        //        cout << "Pointers " << a << " " << b << endl;
+        
+        for (k=lbnd[2]; k<ubnd[2]; ++k) {
+          for (j=lbnd[1]; j<ubnd[1]; ++j) {
+            for (i=lbnd[0]; i<ubnd[0]; ++i) {
+              ind = i + lsh[0] * (j + lsh[1] * k);
+              res += (a[ind]) * (b[ind]);
+              chksum += a[ind];
+              chksum2 += b[ind];
+            }
+          }
+        }
+      }
+//       cout << "before (" 
+//            << lbnd[0] << "," 
+//            << lbnd[1] << "," 
+//            << lbnd[2] << "), (" 
+//            << ubnd[0] << "," 
+//            << ubnd[1] << "," 
+//            << ubnd[2] << "): " 
+//            << res << "   " << chksum << "," << chksum2 << endl;
+      global_sum (cctkGH, &res);
+      //      cout << "after: " << res << endl;
+      realoutput = res;
+      
+    }
+
+    
+  } // namespace common
+
+
+  // Register this solver
+  void CarpetCG_register ()
+  {
+    int const ierr = TATelliptic_RegisterSolver 
+      (CarpetCG_solve, "CarpetCG");
+    assert (!ierr);
+  }
+
+
+
+  // Solve 
+  int
+  CarpetCG_solve (cGH const * restrict const cctkGH,
+                  int const * restrict const var,
+                  int const * restrict const res,
+                  int const nvars,
+                  int const options_table,
+                  calcfunc const calcres,
+                  calcfunc const applybounds,
+                  void * const userdata)
+  {
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+      
+    // Domain descriptors 
+    cGroup groupdata;
+    cGroupDynamicData groupdyndata;
+    int gsh[dim], lsh[dim], nghostzones[dim], bbox[2*dim];
+    int lbnd[dim], ubnd[dim];
+  
+    // Options 
+
+    CCTK_INT maxiters;
+    CCTK_INT maxiters2;
+    CCTK_INT smaxiters;
+    CCTK_REAL stepsize;
+    CCTK_REAL maxstepsize;
+    CCTK_REAL minerror;
+    CCTK_REAL sminerror;
+      
+    int iter;
+    int iter2;
+    int siter;
+      
+    int do_abort;
+      
+    CCTK_REAL alpha;              // secant step size 
+    CCTK_REAL alpha_old, sum_alpha;
+    CCTK_REAL beta;               // CG step size
+    CCTK_REAL delta_0;            // initial A-norm of residual
+    CCTK_REAL delta_new, delta_d, delta_old, delta_mid;
+    CCTK_REAL epsilon;            // norm of residual 
+    CCTK_REAL eta, eta_prev;
+  
+    time_t nexttime, currenttime;
+  
+    CCTK_REAL gsize;             // global grid size (no. points as real)
+  
+    int n, nn;
+    int d, f;
+      
+    int nelems;
+    int ierr;  
+  
+    // Security check 
+    if (! CCTK_IsThornActive(CCTK_THORNSTRING)) {
+      CCTK_WARN (0, "Thorn " CCTK_THORNSTRING " has not been activated.  It is therefore not possible to call CarpetCG_solve.");
+    }
+  
+  
+    // Check arguments 
+    assert (cctkGH);
+  
+    assert (nvars > 0);
+    assert (var);
+    assert (res);
+    for (n=0; n<nvars; ++n) {
+      assert (var[n] >= 0 && var[n] < CCTK_NumVars());
+      assert (CCTK_GroupTypeFromVarI(var[n]) == CCTK_GF
+              || CCTK_GroupTypeFromVarI(var[n]) == CCTK_ARRAY);
+      assert (CCTK_GroupDimFromVarI(var[n]) <= dim);
+      assert (CCTK_VarTypeI(var[n]) == CCTK_VARIABLE_REAL);
+      assert (CCTK_QueryGroupStorageI(cctkGH, CCTK_GroupIndexFromVarI(var[n])));
+    }
+    for (n=0; n<nvars; ++n) {
+      assert (res[n] >= 0 && res[n] < CCTK_NumVars());
+      assert (CCTK_GroupTypeFromVarI(res[n]) == CCTK_GF
+              || CCTK_GroupTypeFromVarI(res[n]) == CCTK_ARRAY);
+      assert (CCTK_GroupDimFromVarI(res[n]) <= dim);
+      assert (CCTK_VarTypeI(res[n]) == CCTK_VARIABLE_REAL);
+      assert (CCTK_QueryGroupStorageI(cctkGH, CCTK_GroupIndexFromVarI(res[n])));
+    }
+    for (n=0; n<nvars; ++n) {
+      assert (var[n] != res[n]);
+      for (nn=0; nn<n; ++nn) {
+        assert (var[nn] != var[n]);
+        assert (var[nn] != res[n]);
+        assert (res[nn] != var[n]);
+        assert (res[nn] != res[n]);
+      }
+    }      
+      
+    assert (options_table >= 0);
+      
+    assert (calcres);
+    assert (applybounds);
+
+    nelems = Util_TableGetIntArray
+      (options_table, 2*dim, &(common::nboundaryzones[0]), "nboundaryzones");
+    if (nelems == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+      for (d=0; d<dim; ++d) {
+        common::nboundaryzones[2*d  ] = -1;
+        common::nboundaryzones[2*d+1] = -1;
+      }
+    } else if (nelems != 2*dim) {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Options table key \"nboundaryzones\" is not an integer array of length %d", 2*dim);
+      return -1;      
+    }
+
+    (common::fromindex).reserve(cg_maxsolvevars);
+    (common::toindex).reserve(cg_maxsolvevars);
+    (common::factors).reserve(nvars);
+    common::factor = 0.5;
+    for (int i = 0; i < nvars; i++)
+    {
+      common::factors[i] = 1.0;
+    }
+  
+    maxiters = 1000;
+    maxiters2 = 1000;
+    minerror = 1.0e-8;
+    smaxiters = 10;
+    sminerror = 1.0e-4;
+    stepsize = 1.0e-6;
+    maxstepsize = 1.0;
+      
+    ierr = Util_TableGetInt (options_table, &maxiters, "maxiters");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetInt (options_table, &maxiters2, "maxiters2");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetInt (options_table, &smaxiters, "smaxiters");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetReal (options_table, &minerror, "minerror");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetReal (options_table, &sminerror, "sminerror");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetReal (options_table, &stepsize, "stepsize");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetReal (options_table, &maxstepsize, "maxstepsize");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+      
+    assert (maxiters >= 0);
+    assert (maxiters2 >= 0);
+    assert (smaxiters >= 0);
+    assert (minerror >= 0);
+    assert (sminerror >= 0);
+    assert (stepsize > 0);
+
+    // Switch to global mode
+    BEGIN_GLOBAL_MODE(cctkGH) {
+      
+      // Fill common block
+      common::var = var;
+      common::res = res;
+      common::nvars = nvars;
+      common::options_table = options_table;
+      common::calcres = calcres;
+      common::applybnds = applybounds;
+      common::userdata = userdata;
+      
+      if (verbose || veryverbose) {
+        CCTK_VInfo (CCTK_THORNSTRING,
+                    "Solving through nonlinear conjugate gradient iterations");
+      }
+      
+      
+      
+      /*
+       * Literature:
+       *
+       * Jonathan Richard Shewchuk, "An Introduction to the Conjugate
+       * Gradient Method Without the Agonizing Pain", Technical Report
+       * CMU-CS-94-125, Aug. 1994
+       *
+       * Available from http://www-2.cs.cmu.edu/~jrs/jrspapers.html
+       */
+      
+      /*
+       * Notation:
+       *
+       *    here        there   meaning
+       *                
+       *    iter        i       current cg iteration
+       *    siter       j       current secant iteration
+       *    iter2       k       current cg iteration since last restart
+       *    maxiters2   n       PARA cg iterations before restart
+       *                alpha   secant step size
+       *                beta    cg step size
+       *                delta   cg distance
+       *                epsilon PARA secant error tolerance (straight epsilon)
+       *    minerror    vareps  PARA cg error tolerance (curved epsilon)
+       *                eta     secant distance
+       *    stepsize    sigma   PARA secand method step parameter
+       *                
+       *    dir         d       cg direction
+       *   -calcres     f'      PARA nonlinear operator
+       *    res         r       cg residual
+       *    prs         s       preconditioned cg residual
+       *    var         x       INIT unknown variable
+       *                
+       *   -wgt         M       preconditioning matrix (not in this version)
+       */
+      
+      /* 
+       * Algorithm:
+       * (Preconditioned nonlinear conjugate gradients with secant and
+       * Polak-Ribière)
+       *
+       *    01. i <= 0
+       *    02. k <= 0
+       *    03. r <= - f'(x)
+       *    04. calculate M \approx f''(x)
+       *    05. s <= M^-1 r
+       *    06. d <= s
+       *    07. delta_new <= r^T d
+       *    08. delta_0 <= delta_new
+       *    09. WHILE i < i_max AND delta_new > vareps^2 delta_0 DO
+       *    10.    j <= 0
+       *    11.    delta_d <= d^T d
+       *    12.    alpha <= - sigma_0
+       *    13.    eta_prev <= [f'(x + sigma_0 d)]^T d
+       *    14.    DO
+       *    15.       eta <= [f'(x)]^T d
+       *    16.       alpha <= alpha   (eta) / (eta_prev - eta)
+       *    17.       x <= x + alpha d
+       *    18.       eta_prev <= eta
+       *    19.       j <= j + 1
+       *    20.    WHILE j < j_max AND alpha^2 delta_d > epsilon^2
+       *    21.    r <= - f'(x)
+       *    22.    delta_old <= delta_new
+       *    23.    delta_mid <= r^T s
+       *    24.    calculate M \approx f''(x)
+       *    25.    s <= M^-1 r
+       *    26.    delta_new <= r^T s
+       *    27.    beta <= (delta_new - delta_mid) / (delta_old)
+       *    28.    k <= k + 1
+       *    29.    IF k = n OR beta <= 0 THEN
+       *    30.       d <= s
+       *    31.       k <= 0
+       *    32.    ELSE
+       *    33.       d <= s + beta d
+       *    34.    i <= i + 1
+       */
+      
+
+      // Pointers to grid variables and temporary storage 
+      vector<int> varptrs(nvars);
+      vector<int> resptrs(nvars);
+      //    CCTK_REAL * restrict * restrict prsptrs;
+      //    CCTK_REAL * restrict * restrict dirptrs;  
+      
+      // Get storage pointers
+      for (n=0; n<nvars; ++n) {
+        varptrs[n] = var[n];
+        resptrs[n] = res[n];
+      }
+  
+      // Allocate temporary memory
+      //    prsptrs = malloc (nvars * sizeof *prsptrs);
+      //    dirptrs = malloc (nvars * sizeof *prsptrs);
+      assert (cg_maxsolvevars >= nvars);
+      vector<int> prsptrs(cg_maxsolvevars);
+      vector<int> dirptrs(cg_maxsolvevars);
+      for (n=0; n<nvars; ++n) {
+        char varname[1000];
+        sprintf(varname,"CarpetCG::prsvars[%d]",n);
+        prsptrs[n] = CCTK_VarIndex(varname);
+        assert ((prsptrs[n] >= 0)&&(prsptrs[n] < CCTK_NumVars()));
+        sprintf(varname,"CarpetCG::dirvars[%d]",n);
+        dirptrs[n] = CCTK_VarIndex(varname);
+        assert ((dirptrs[n] >= 0)&&(dirptrs[n] < CCTK_NumVars()));
+      }
+  
+  
+  
+      nexttime = time(0);
+  
+      /* 01. i <= 0 */
+      iter = 0;
+  
+      /* 02. k <= 0 */
+      iter2 = 0;
+  
+      /* 03. r <= - f'(x) */
+      /* 04. calculate M \approx f''(x) */
+      common::ierr = 0;
+      CallLocalFunction((cGH *)cctkGH, common::call_calcres);
+      ierr = common::ierr;
+      assert (! ierr);
+
+      for (int n = 0; n < nvars; n++) {
+        common::fromindex[n] = resptrs[n];
+      }
+      for (int n = nvars; n < cg_maxsolvevars; n++) {
+        common::fromindex[n] = -1;
+      }
+      CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign);
+    
+      /* 05. s <= M^-1 r */
+      // No preconditioning
+      for (int n = 0; n < nvars; n++) {
+        common::fromindex[n] = resptrs[n];
+        common::toindex[n] = prsptrs[n];
+      }
+      for (int n = nvars; n < cg_maxsolvevars; n++) {
+        common::fromindex[n] = -1;
+        common::toindex[n] = -1;
+      }
+    
+      //    cout << "dirptrs " << dirptrs[0] << endl;
+    
+      CallLocalFunction((cGH *)cctkGH, common::call_apply_preconditioner);
+    
+      //    cout << "dirptrs " << dirptrs[0] << endl;
+      
+      /* 06. d <= s */
+      for (int n = 0; n < nvars; n++) {
+        common::fromindex[n] = prsptrs[n];
+        common::toindex[n] = dirptrs[n];
+      }
+      for (int n = nvars; n < cg_maxsolvevars; n++) {
+        common::fromindex[n] = -1;
+        common::toindex[n] = -1;
+      }
+    
+      //    cout << "dirptrs " << dirptrs[0] << endl;
+    
+      CallLocalFunction((cGH *)cctkGH, common::call_copy);
+  
+      //    cout << "dirptrs " << dirptrs[0] << endl;
+    
+      /* 07. delta_new <= r^T d */
+      //    output (cctkGH, var, res, wgt, nvars, iter);
+      common::realoutput_count = 0;
+      for (int n = 0; n < nvars; n++) {
+        common::fromindex[n] = prsptrs[n];
+        common::toindex[n] = resptrs[n];
+      }
+      for (int n = nvars; n < cg_maxsolvevars; n++) {
+        common::fromindex[n] = -1;
+        common::toindex[n] = -1;
+      }
+    
+      //    cout << "dirptrs " << dirptrs[0] << endl;
+    
+      CallLocalFunction((cGH *)cctkGH, common::call_dot_product);
+      delta_new = common::realoutput;
+      gsize = common::realoutput_count;
+      //    cout << "delta_new " << delta_new << " gsize " << gsize << endl;
+    
+      for (int n = 0; n < nvars; n++) {
+        common::fromindex[n] = resptrs[n];
+        common::toindex[n] = resptrs[n];
+      }
+      for (int n = nvars; n < cg_maxsolvevars; n++) {
+        common::fromindex[n] = -1;
+        common::toindex[n] = -1;
+      }
+    
+      CallLocalFunction((cGH *)cctkGH, common::call_dot_product);
+      epsilon = common::realoutput;
+      //     cout << "epsilon " << epsilon << endl;
+    
+      /* 08. delta_0 <= delta_new */
+      delta_0 = delta_new;
+  
+      /* 09. WHILE i < i_max AND delta_new > vareps^2 delta_0 DO */
+      while (iter < maxiters && epsilon / (nvars * gsize) > pow(minerror,2)) {
+    
+        if (verbose || veryverbose) {
+          currenttime = time(0);
+          if (veryverbose || (iter % outevery == 0 && currenttime >= nexttime)) {
+            CCTK_VInfo (CCTK_THORNSTRING,
+                        "Iteration %d (%d since restart): residual is %g (%g,%d,%g)",
+                        iter, iter2, (double)sqrt(epsilon / (nvars * gsize)),epsilon,nvars,gsize);
+            if (outeveryseconds > 0) {
+              while (nexttime <= currenttime) nexttime += outeveryseconds;
+            }
+          }
+        }
+    
+        /* 10. j <= 0 */
+        siter = 0;
+
+        /* 11. delta_d <= d^T d */
+        for (int n = 0; n < nvars; n++) {
+          common::fromindex[n] = dirptrs[n];
+          common::toindex[n] = dirptrs[n];
+        }
+        for (int n = nvars; n < cg_maxsolvevars; n++) {
+          common::fromindex[n] = -1;
+          common::toindex[n] = -1;
+        }
+
+        CallLocalFunction((cGH *)cctkGH, common::call_dot_product);
+        delta_d = common::realoutput;
+
+        //       cout << "delta_d " << delta_d << endl;
+    
+        /* 12. alpha <= - sigma_0 */
+        alpha = - stepsize;
+        sum_alpha = alpha;
+        do_abort = 0;
+    
+        /* 13. eta_prev <= [f'(x + sigma_0 d)]^T d */
+        for (int n = 0; n < nvars; n++) {
+          common::fromindex[n] = dirptrs[n];
+          common::toindex[n] = varptrs[n];
+        }
+        for (int n = nvars; n < cg_maxsolvevars; n++) {
+          common::fromindex[n] = -1;
+          common::toindex[n] = -1;
+        }
+        common::realconstant = -alpha;
+        //       cout << "Add scaled, const " << -alpha << endl;
+
+        CallLocalFunction((cGH *)cctkGH, common::call_add_scaled);
+        common::ierr = 0;
+        CallLocalFunction((cGH *)cctkGH, common::call_applybnds);
+        assert (! ierr);
+    
+        ierr = 0;
+        CallLocalFunction((cGH *)cctkGH, common::call_calcres);
+        assert (! ierr);
+
+        for (int n = 0; n < nvars; n++) {
+          common::fromindex[n] = resptrs[n];
+        }
+        for (int n = nvars; n < cg_maxsolvevars; n++) {
+          common::fromindex[n] = -1;
+        }
+        CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign);
+    
+        for (int n = 0; n < nvars; n++) {
+          common::fromindex[n] = dirptrs[n];
+          common::toindex[n] = resptrs[n];
+        }
+        for (int n = nvars; n < cg_maxsolvevars; n++) {
+          common::fromindex[n] = -1;
+          common::toindex[n] = -1;
+        }
+        CallLocalFunction((cGH *)cctkGH, common::call_dot_product);
+        eta = - common::realoutput;
+
+        //       cout << "eta " << eta << endl;
+
+        for (int n = 0; n < nvars; n++) {
+          common::fromindex[n] = dirptrs[n];
+          common::toindex[n] = varptrs[n];
+        }
+        for (int n = nvars; n < cg_maxsolvevars; n++) {
+          common::fromindex[n] = -1;
+          common::toindex[n] = -1;
+        }
+        common::realconstant = alpha;
+        //       cout << "Add scaled, const " << alpha << endl;
+
+        CallLocalFunction((cGH *)cctkGH, common::call_add_scaled);
+        ierr = 0;
+        CallLocalFunction((cGH *)cctkGH, common::call_applybnds);
+        assert (! ierr);
+    
+        eta_prev = eta;
+    
+        /* 14. DO */
+        do {
+    
+          if (veryverbose) {
+            CCTK_VInfo
+              (CCTK_THORNSTRING,
+               "   Secant iteration %d: step size is %g, orthogonality is %g",
+               siter,
+               (double)alpha, (double)sqrt(fabs(eta_prev) / (nvars * gsize)));
+          }
+      
+          /* 15. eta <= [f'(x)]^T d */
+          ierr = 0;
+          CallLocalFunction((cGH *)cctkGH, common::call_calcres);
+          assert (! ierr);
+
+          for (int n = 0; n < nvars; n++) {
+            common::fromindex[n] = resptrs[n];
+          }
+          for (int n = nvars; n < cg_maxsolvevars; n++) {
+            common::fromindex[n] = -1;
+          }
+          CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign);
+
+          for (int n = 0; n < nvars; n++) {
+            common::fromindex[n] = dirptrs[n];
+            common::toindex[n] = resptrs[n];
+          }
+          for (int n = nvars; n < cg_maxsolvevars; n++) {
+            common::fromindex[n] = -1;
+            common::toindex[n] = -1;
+          }
+
+          CallLocalFunction((cGH *)cctkGH, common::call_dot_product);
+          eta = - common::realoutput;
+
+          //         cout << "eta " << eta << endl;
+      
+          /* 16. alpha <= alpha   (eta) / (eta_prev - eta) */
+          alpha_old = alpha;
+          alpha *= eta / (eta_prev - eta);
+          sum_alpha += alpha;
+          if (veryverbose) {
+            CCTK_VInfo
+              (CCTK_THORNSTRING,
+               "   Changing step size, iteration %d: was %g, now %g (%g, %g)", 
+               siter, alpha_old, alpha, eta, eta_prev);
+          }
+          assert (sum_alpha > - maxstepsize);
+          if (sum_alpha > maxstepsize) {
+            alpha = maxstepsize - alpha_old;
+            do_abort = 1;
+            if (verbose) {
+              CCTK_VInfo
+                (CCTK_THORNSTRING,
+                 "   Secant iteration %d: limiting total step size", siter);
+            }
+          }
+      
+          /* 17. x <= x + alpha d */
+
+          for (int n = 0; n < nvars; n++) {
+            common::fromindex[n] = dirptrs[n];
+            common::toindex[n] = varptrs[n];
+          }
+          for (int n = nvars; n < cg_maxsolvevars; n++) {
+            common::fromindex[n] = -1;
+            common::toindex[n] = -1;
+          }
+          common::realconstant = alpha;
+          //         cout << "Add scaled, const " << alpha << endl;
+
+          CallLocalFunction((cGH *)cctkGH, common::call_add_scaled);
+          ierr = 0;
+          CallLocalFunction((cGH *)cctkGH, common::call_applybnds);
+          assert (! ierr);
+      
+          /* 18. eta_prev <= eta */
+          eta_prev = eta;
+      
+          /* 19. j <= j + 1 */
+          ++ siter;
+      
+          /* 20. WHILE j < j_max AND alpha^2 delta_d > epsilon^2 */
+        } while (siter < smaxiters
+                 && pow(alpha,2) * delta_d > pow(sminerror,2)
+                 && !do_abort);
+    
+        if (veryverbose) {
+          CCTK_VInfo
+            (CCTK_THORNSTRING,
+             "   Secant iteration %d: step size is %g, orthogonality is %g",
+             siter,
+             (double)alpha, (double)sqrt(fabs(eta_prev) / (nvars * gsize)));
+        }
+    
+        /* 21. r <= - f'(x) */
+        /* 24. calculate M \approx f''(x) */
+        ierr = 0;
+        CallLocalFunction((cGH *)cctkGH, common::call_calcres);
+        assert (! ierr);
+
+        for (int n = 0; n < nvars; n++) {
+          common::fromindex[n] = resptrs[n];
+        }
+        for (int n = nvars; n < cg_maxsolvevars; n++) {
+          common::fromindex[n] = -1;
+        }
+        CallLocalFunction((cGH *)cctkGH, common::call_correct_residual_sign);
+    
+        /* 23. delta_mid <= r^T s */
+        for (int n = 0; n < nvars; n++) {
+          common::fromindex[n] = prsptrs[n];
+          common::toindex[n] = resptrs[n];
+        }
+        for (int n = nvars; n < cg_maxsolvevars; n++) {
+          common::fromindex[n] = -1;
+          common::toindex[n] = -1;
+        }
+
+        CallLocalFunction((cGH *)cctkGH, common::call_dot_product);
+        delta_mid = common::realoutput;
+        //       cout << "delta_mid " << delta_mid << endl;
+      
+        /* 25. s <= M^-1 r */
+        for (int n = 0; n < nvars; n++) {
+          common::fromindex[n] = resptrs[n];
+          common::toindex[n] = prsptrs[n];
+        }
+        for (int n = nvars; n < cg_maxsolvevars; n++) {
+          common::fromindex[n] = -1;
+          common::toindex[n] = -1;
+        }
+      
+        CallLocalFunction((cGH *)cctkGH, common::call_apply_preconditioner);
+
+        /* 22. delta_old <= delta_new */
+        delta_old = delta_new;
+    
+        /* 26. delta_new <= r^T s */
+        //       output (cctkGH, var, res, wgt, nvars, iter+1);
+        for (int n = 0; n < nvars; n++) {
+          common::fromindex[n] = prsptrs[n];
+          common::toindex[n] = resptrs[n];
+        }
+        for (int n = nvars; n < cg_maxsolvevars; n++) {
+          common::fromindex[n] = -1;
+          common::toindex[n] = -1;
+        }
+
+        CallLocalFunction((cGH *)cctkGH, common::call_dot_product);
+        delta_new = common::realoutput;
+        //       cout << "delta_new " << delta_new << endl;
+
+        for (int n = 0; n < nvars; n++) {
+          common::fromindex[n] = resptrs[n];
+          common::toindex[n] = resptrs[n];
+        }
+        for (int n = nvars; n < cg_maxsolvevars; n++) {
+          common::fromindex[n] = -1;
+          common::toindex[n] = -1;
+        }
+
+        CallLocalFunction((cGH *)cctkGH, common::call_dot_product);
+        epsilon = common::realoutput;
+        //       cout << "epsilon " << epsilon << endl;
+    
+        /* 27. beta <= (delta_new - delta_mid) / (delta_old) */
+        beta = (delta_new - delta_mid) / delta_old;
+    
+        /* 28. k <= k + 1 */
+        ++ iter2;
+    
+        /* 29. IF k = n OR beta <= 0 THEN */
+        if (iter2 >= maxiters2 || beta <= 0 || do_abort) {
+          /* Restart */
+      
+          /* 30. d <= s */
+          for (int n = 0; n < nvars; n++) {
+            common::fromindex[n] = prsptrs[n];
+            common::toindex[n] = dirptrs[n];
+          }
+          for (int n = nvars; n < cg_maxsolvevars; n++) {
+            common::fromindex[n] = -1;
+            common::toindex[n] = -1;
+          }
+          CallLocalFunction((cGH *)cctkGH, common::call_copy);
+      
+          /* 31. k <= 0 */
+          iter2 = 0;
+      
+          /* 32. ELSE */
+        } else {
+          /* Continue with CG iterations */
+      
+          /* 33. d <= s + beta d */
+          for (int n = 0; n < nvars; n++) {
+            common::fromindex[n] = prsptrs[n];
+            common::toindex[n] = dirptrs[n];
+          }
+          for (int n = nvars; n < cg_maxsolvevars; n++) {
+            common::fromindex[n] = -1;
+            common::toindex[n] = -1;
+          }
+          common::realconstant = beta;
+          //        cout << "beta " << beta << endl;
+        
+          CallLocalFunction((cGH *)cctkGH, common::call_scale_and_add);
+      
+        } /* if iter2 */
+    
+        /* 34. i <= i + 1 */
+        ++ iter;
+    
+      } /* for iter */
+  
+      if (verbose || veryverbose) {
+        CCTK_VInfo (CCTK_THORNSTRING,
+                    "Iteration %d (%d since restart): residual is %g",
+                    iter, iter2, (double)sqrt(epsilon / (nvars * gsize)));
+      }
+  
+    /* Free temporary memory */
+      //     for (n=0; n<nvars; ++n) {
+      //       free (prsptrs[n]);
+      //       free (dirptrs[n]);
+      //     }
+      
+      //     free (varptrs);
+      //     free (resptrs);
+      //     free (prsptrs);
+      //     free (dirptrs);
+      
+    } END_GLOBAL_MODE;
+    
+    return 0;
+  }
+    
+} // namespace CarpetCG
+  
+  
diff --git a/CarpetDev/CarpetCG/src/make.code.defn b/CarpetDev/CarpetCG/src/make.code.defn
new file mode 100644
index 000000000..3db692761
--- /dev/null
+++ b/CarpetDev/CarpetCG/src/make.code.defn
@@ -0,0 +1,9 @@
+# Main make.code.defn file for thorn CarpetCG
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetCG/src/make.code.defn,v 1.1.1.1 2003/11/19 10:40:47 hawke Exp $
+
+# Source files in this directory
+SRCS = CG.cc
+
+# Subdirectories containing source files
+SUBDIRS = 
+
diff --git a/CarpetDev/CarpetJacobi/README b/CarpetDev/CarpetJacobi/README
new file mode 100644
index 000000000..74967b13f
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/README
@@ -0,0 +1,11 @@
+CVS info   : $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/README,v 1.1 2003/09/02 14:35:57 schnetter Exp $
+
+Cactus Code Thorn CarpetJacobi
+Thorn Author(s)     : Erik Schnetter <schnetter@uni-tuebingen.de>
+Thorn Maintainer(s) : Erik Schnetter <schnetter@uni-tuebingen.de>
+--------------------------------------------------------------------------
+
+Purpose of the thorn:
+
+Solve an elliptic equation using Jacobi iterations with varying step
+size.
diff --git a/CarpetDev/CarpetJacobi/doc/documentation.tex b/CarpetDev/CarpetJacobi/doc/documentation.tex
new file mode 100644
index 000000000..11d9ffd8d
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/doc/documentation.tex
@@ -0,0 +1,144 @@
+% *======================================================================*
+%  Cactus Thorn template for ThornGuide documentation
+%  Author: Ian Kelley
+%  Date: Sun Jun 02, 2002
+%  $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/doc/documentation.tex,v 1.1 2003/09/02 14:35:57 schnetter Exp $                                                             
+%
+%  Thorn documentation in the latex file doc/documentation.tex 
+%  will be included in ThornGuides built with the Cactus make system.
+%  The scripts employed by the make system automatically include 
+%  pages about variables, parameters and scheduling parsed from the 
+%  relevant thorn CCL files.
+%  
+%  This template contains guidelines which help to assure that your     
+%  documentation will be correctly added to ThornGuides. More 
+%  information is available in the Cactus UsersGuide.
+%                                                    
+%  Guidelines:
+%   - Do not change anything before the line
+%       % START CACTUS THORNGUIDE",
+%     except for filling in the title, author, date etc. fields.
+%        - Each of these fields should only be on ONE line.
+%        - Author names should be separated with a \\ or a comma
+%   - You can define your own macros, but they must appear after
+%     the START CACTUS THORNGUIDE line, and must not redefine standard 
+%     latex commands.
+%   - To avoid name clashes with other thorns, 'labels', 'citations', 
+%     'references', and 'image' names should conform to the following 
+%     convention:          
+%       ARRANGEMENT_THORN_LABEL
+%     For example, an image wave.eps in the arrangement CactusWave and 
+%     thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps
+%   - Graphics should only be included using the graphicx package. 
+%     More specifically, with the "\includegraphics" command.  Do
+%     not specify any graphic file extensions in your .tex file. This 
+%     will allow us (later) to create a PDF version of the ThornGuide
+%     via pdflatex.
+%   - References should be included with the latex "\bibitem" command. 
+%   - Use \begin{abstract}...\end{abstract} instead of \abstract{...}
+%   - Do not use \appendix, instead include any appendices you need as 
+%     standard sections. 
+%   - For the benefit of our Perl scripts, and for future extensions, 
+%     please use simple latex.     
+%
+% *======================================================================* 
+% 
+% Example of including a graphic image:
+%    \begin{figure}[ht]
+% 	\begin{center}
+%    	   \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure}
+% 	\end{center}
+% 	\caption{Illustration of this and that}
+% 	\label{MyArrangement_MyThorn_MyLabel}
+%    \end{figure}
+%
+% Example of using a label:
+%   \label{MyArrangement_MyThorn_MyLabel}
+%
+% Example of a citation:
+%    \cite{MyArrangement_MyThorn_Author99}
+%
+% Example of including a reference
+%   \bibitem{MyArrangement_MyThorn_Author99}
+%   {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999), 
+%   1--16. {\tt http://www.nowhere.com/}}
+%
+% *======================================================================* 
+
+% If you are using CVS use this line to give version information
+% $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/doc/documentation.tex,v 1.1 2003/09/02 14:35:57 schnetter Exp $
+
+\documentclass{article}
+
+% Use the Cactus ThornGuide style file
+% (Automatically used from Cactus distribution, if you have a 
+%  thorn without the Cactus Flesh download this from the Cactus
+%  homepage at www.cactuscode.org)
+\usepackage{../../../../doc/latex/cactus}
+
+\begin{document}
+
+% The author of the documentation
+\author{Erik Schnetter \textless schnetter@uni-tuebingen.de\textgreater}
+
+% The title of the document (not necessarily the name of the Thorn)
+\title{}
+
+% the date your document was last changed, if your document is in CVS, 
+% please use:
+%    \date{$ $Date: 2003/09/02 14:35:57 $ $}
+\date{August 31 2003}
+
+\maketitle
+
+% Do not delete next line
+% START CACTUS THORNGUIDE
+
+% Add all definitions used in this documentation here 
+%   \def\mydef etc
+
+% Add an abstract for this thorn's documentation
+\begin{abstract}
+
+\end{abstract}
+
+% The following sections are suggestive only.
+% Remove them or add your own.
+
+\section{Introduction}
+
+\section{Physical System}
+
+\section{Numerical Implementation}
+
+\section{Using This Thorn}
+
+\subsection{Obtaining This Thorn}
+
+\subsection{Basic Usage}
+
+\subsection{Special Behaviour}
+
+\subsection{Interaction With Other Thorns}
+
+\subsection{Examples}
+
+\subsection{Support and Feedback}
+
+\section{History}
+
+\subsection{Thorn Source Code}
+
+\subsection{Thorn Documentation}
+
+\subsection{Acknowledgements}
+
+
+\begin{thebibliography}{9}
+
+\end{thebibliography}
+
+% Do not delete next line
+% END CACTUS THORNGUIDE
+
+\end{document}
diff --git a/CarpetDev/CarpetJacobi/interface.ccl b/CarpetDev/CarpetJacobi/interface.ccl
new file mode 100644
index 000000000..6b61da31a
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/interface.ccl
@@ -0,0 +1,9 @@
+# Interface definition for thorn CarpetJacobi
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/interface.ccl,v 1.1 2003/09/02 14:35:57 schnetter Exp $
+
+IMPLEMENTS: CarpetJacobi
+
+INHERITS: TATelliptic
+
+USES INCLUDE HEADER: carpet.hh
+USES INCLUDE HEADER: TATelliptic.h
diff --git a/CarpetDev/CarpetJacobi/par/charge_te_cjac.par b/CarpetDev/CarpetJacobi/par/charge_te_cjac.par
new file mode 100644
index 000000000..8db9c4cf6
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/par/charge_te_cjac.par
@@ -0,0 +1,60 @@
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/par/charge_te_cjac.par,v 1.4 2004/03/23 11:59:20 schnetter Exp $
+# Erik Schnetter <schnetter@uni-tuebingen.de>
+
+ActiveThorns = "Boundary CartGrid3D CoordBase SymBase CarpetIOASCII IOBasic IOUtil Time Carpet CarpetLib CarpetRegrid CarpetReduce NaNCatcher IDScalarWave WaveToyC IDSWTEsimple TATelliptic CarpetJacobi"
+
+!DESC "Charged sphere initial data, solved with TATelliptic/CarpetJacobi"
+
+cactus::cctk_full_warnings      = yes
+cactus::cctk_timer_output       = full
+
+#Cactus::cctk_itlast             = 10
+Cactus::cctk_itlast             = 0
+ 
+Time::dtfac                     = 0.5
+
+driver::global_nsize            = 21
+#driver::global_nsize           = 41
+
+Carpet::verbose                 = yes
+#Carpet::veryverbose             = yes
+
+Carpet::max_refinement_levels   = 2
+#Carpet::prolongation_order_space = 3
+Carpet::prolongation_order_time = 2
+Carpet::init_3_timelevels       = yes
+
+CarpetRegrid::refinement_levels = 1
+CarpetRegrid::refined_regions   = manual-gridpoint-list
+CarpetRegrid::gridpoints        = "[[ ([1,1,1]:[20,20,20]:[1,1,1]) ]]"
+CarpetRegrid::outerbounds       = "[[ [[1,0],[1,0],[1,0]] ]]"
+
+grid::type                      = byspacing
+grid::domain                    = octant
+grid::dxyz                      = 1.2
+grid::avoid_origin              = no
+
+IDScalarWave::initial_data      = charge-TATelliptic-simple
+IDSWTEsimple::solver            = CarpetJacobi
+#IDSWTEsimple::options           = "maxiters=10000 minerror=1e-4 factor=1e-2"
+#IDSWTEsimple::options           = "maxiters=100000 minerror=1e-8 factor=1e-2"
+IDSWTEsimple::options           = "maxiters=10000 minerror=1e-8"
+IDSWTEsimple::radius            = 5.5
+IDSWTEsimple::charge            = 1.0
+
+WaveToyC::bound                 = radiation
+
+CarpetJacobi::verbose           = yes
+
+IO::out_dir                     = "charge_te_cjac"
+IO::out_fileinfo                = none
+
+IOBasic::outInfo_every          = 1
+IOBasic::outInfo_vars           = "wavetoy::phi IDSWTEsimple::residual"
+
+IOBasic::outScalar_every        = 1
+IOBasic::outScalar_vars         = "wavetoy::phi IDSWTEsimple::residual"
+IOBasic::outScalar_style        = "gnuplot"
+
+CarpetIOASCII::out1D_every      = 1
+CarpetIOASCII::out1D_vars       = "wavetoy::phi IDSWTEsimple::residual"
diff --git a/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par b/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par
new file mode 100644
index 000000000..b26fe4d48
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par
@@ -0,0 +1,60 @@
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par,v 1.4 2004/03/23 11:59:20 schnetter Exp $
+# Erik Schnetter <schnetter@uni-tuebingen.de>
+
+ActiveThorns = "Boundary CartGrid3D CoordBase SymBase CarpetIOASCII IOBasic IOUtil Time Carpet CarpetLib CarpetRegrid CarpetReduce NaNCatcher IDScalarWave WaveToyC IDSWTEsimple TATelliptic CarpetJacobi"
+
+!DESC "Charged sphere initial data, solved with TATelliptic/CarpetJacobi"
+
+cactus::cctk_full_warnings      = yes
+cactus::cctk_timer_output       = full
+
+#Cactus::cctk_itlast             = 10
+Cactus::cctk_itlast             = 0
+ 
+Time::dtfac                     = 0.5
+
+driver::global_nsize            = 21
+#driver::global_nsize           = 41
+
+Carpet::verbose                 = yes
+#Carpet::veryverbose             = yes
+
+Carpet::max_refinement_levels   = 2
+#Carpet::prolongation_order_space = 3
+Carpet::prolongation_order_time = 2
+Carpet::init_3_timelevels       = yes
+
+CarpetRegrid::refinement_levels = 2
+CarpetRegrid::refined_regions   = manual-gridpoint-list
+CarpetRegrid::gridpoints        = "[[ ([1,1,1]:[20,20,20]:[1,1,1]) ]]"
+CarpetRegrid::outerbounds       = "[[ [[1,0],[1,0],[1,0]] ]]"
+
+grid::type                      = byspacing
+grid::domain                    = octant
+grid::dxyz                      = 1.2
+grid::avoid_origin              = no
+
+IDScalarWave::initial_data      = charge-TATelliptic-simple
+IDSWTEsimple::solver            = CarpetJacobi
+#IDSWTEsimple::options           = "maxiters=10000 minerror=1e-4 factor=1e-2"
+#IDSWTEsimple::options           = "maxiters=100000 minerror=1e-8 factor=1e-2"
+IDSWTEsimple::options           = "maxiters=10000 minerror=1e-8"
+IDSWTEsimple::radius            = 5.5
+IDSWTEsimple::charge            = 1.0
+
+WaveToyC::bound                 = radiation
+
+CarpetJacobi::verbose           = yes
+
+IO::out_dir                     = "charge_te_cjac_2"
+IO::out_fileinfo                = none
+
+IOBasic::outInfo_every          = 1
+IOBasic::outInfo_vars           = "wavetoy::phi IDSWTEsimple::residual"
+
+IOBasic::outScalar_every        = 1
+IOBasic::outScalar_vars         = "wavetoy::phi IDSWTEsimple::residual"
+IOBasic::outScalar_style        = "gnuplot"
+
+CarpetIOASCII::out1D_every      = 1
+CarpetIOASCII::out1D_vars       = "wavetoy::phi IDSWTEsimple::residual"
diff --git a/CarpetDev/CarpetJacobi/param.ccl b/CarpetDev/CarpetJacobi/param.ccl
new file mode 100644
index 000000000..bdeecf41b
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/param.ccl
@@ -0,0 +1,32 @@
+# Parameter definitions for thorn CarpetJacobi
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/param.ccl,v 1.2 2004/01/25 14:57:29 schnetter Exp $
+
+BOOLEAN verbose "Produce log output while running" STEERABLE=always
+{
+} "no"
+
+BOOLEAN veryverbose "Produce much log output while running" STEERABLE=always
+{
+} "no"
+
+CCTK_INT outevery "Produce output every n iterations" STEERABLE=always
+{
+  1:* :: ""
+} 1
+
+CCTK_INT outeveryseconds "Produce output every n seconds" STEERABLE=always
+{
+  0:* :: ""
+} 10
+
+
+
+CCTK_INT solve_minlevel "First level on which to solve" STEERABLE=always
+{
+  0:* :: ""
+} 0
+
+CCTK_INT reflevelpower "Power of the refinement factor" STEERABLE=always
+{
+  0:* :: ""
+} 2
diff --git a/CarpetDev/CarpetJacobi/schedule.ccl b/CarpetDev/CarpetJacobi/schedule.ccl
new file mode 100644
index 000000000..009b40f7a
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/schedule.ccl
@@ -0,0 +1,7 @@
+# Schedule definitions for thorn CarpetJacobi
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/schedule.ccl,v 1.1 2003/09/02 14:35:57 schnetter Exp $
+
+SCHEDULE CarpetJacobi_register AT startup
+{
+  LANG: C
+} "Register the elliptic solver"
diff --git a/CarpetDev/CarpetJacobi/src/Jacobi.cc b/CarpetDev/CarpetJacobi/src/Jacobi.cc
new file mode 100644
index 000000000..c326b66cd
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/src/Jacobi.cc
@@ -0,0 +1,536 @@
+// $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/src/Jacobi.cc,v 1.3 2004/01/25 14:57:29 schnetter Exp $
+
+#include <cassert>
+#include <cmath>
+#include <cstdio>
+#include <ctime>
+#include <iostream>
+#include <vector>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
+#include "TATelliptic.h"
+
+#include "carpet.hh"
+
+#include "th.hh"
+
+
+
+namespace Carpet {
+  // TODO: fix this
+  void CycleTimeLevels (const cGH* cctkGH);
+  void Restrict (const cGH* cctkGH);
+};
+
+
+
+namespace CarpetJacobi {
+  
+  using namespace std;
+  using namespace Carpet;
+  
+  
+  
+  extern "C" {
+    void CarpetJacobi_register ();
+  }
+  
+  int CarpetJacobi_solve (const cGH * const cctkGH,
+                          const int * const var,
+                          const int * const res,
+                          const int nvars,
+                          const int options_table,
+                          const calcfunc calcres,
+                          const calcfunc applybnds,
+                          void * const userdata);
+  
+  
+  
+  void CarpetJacobi_register ()
+  {
+    int const ierr = TATelliptic_RegisterSolver
+      (CarpetJacobi_solve, "CarpetJacobi");
+    assert (!ierr);
+  }
+  
+  
+  
+  int CarpetJacobi_solve (const cGH * const cctkGH,
+                          const int * const var,
+                          const int * const res,
+                          const int nvars,
+                          const int options_table,
+                          const calcfunc calcres,
+                          const calcfunc applybnds,
+                          void * const userdata)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    // Error control
+    int ierr;
+    
+    if (! CCTK_IsThornActive(CCTK_THORNSTRING)) {
+      CCTK_WARN (0, "Thorn " CCTK_THORNSTRING " has not been activated.  It is therefore not possible to call CarpetJacobi_solve.");
+    }
+    
+    // Check arguments
+    assert (cctkGH);
+    
+    assert (nvars > 0);
+    assert (var);
+    assert (res);
+    for (int n=0; n<nvars; ++n) {
+      assert (var[n] >= 0 && var[n] < CCTK_NumVars());
+      assert (CCTK_GroupTypeFromVarI(var[n]) == CCTK_GF);
+      assert (CCTK_GroupDimFromVarI(var[n]) == dim);
+      assert (CCTK_VarTypeI(var[n]) == CCTK_VARIABLE_REAL);
+      assert (CCTK_QueryGroupStorageI(cctkGH,
+                                      CCTK_GroupIndexFromVarI(var[n])));
+    }
+    for (int n=0; n<nvars; ++n) {
+      assert (res[n] >= 0 && res[n] < CCTK_NumVars());
+      assert (CCTK_GroupTypeFromVarI(res[n]) == CCTK_GF);
+      assert (CCTK_GroupDimFromVarI(res[n]) == dim);
+      assert (CCTK_VarTypeI(res[n]) == CCTK_VARIABLE_REAL);
+      assert (CCTK_QueryGroupStorageI(cctkGH,
+                                      CCTK_GroupIndexFromVarI(res[n])));
+    }
+    for (int n=0; n<nvars; ++n) {
+      assert (var[n] != res[n]);
+      for (int nn=0; nn<n; ++nn) {
+        assert (var[nn] != var[n]);
+        assert (var[nn] != res[n]);
+        assert (res[nn] != var[n]);
+        assert (res[nn] != res[n]);
+      }
+    }
+    
+    
+    
+    assert (options_table >= 0);
+    
+    assert (calcres);
+    assert (applybnds);
+    
+    // Get options from options table
+    CCTK_INT maxiters = 10000;
+    CCTK_REAL minerror = 1e-8;
+    CCTK_REAL factor = 1.0e-2;    // slow start
+    CCTK_REAL minfactor = 1.0e-8;
+    CCTK_REAL maxfactor = 1.0;
+    CCTK_REAL acceleration = 1.2; // slow acceleration
+    CCTK_REAL deceleration = 0.1; // emergency brake
+    vector<CCTK_REAL> factors (nvars);
+    for (int n=0; n<nvars; ++n) {
+      factors[n] = 1.0;
+    }
+    
+    ierr = Util_TableGetInt (options_table, &maxiters, "maxiters");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetReal (options_table, &minerror, "minerror");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetReal (options_table, &factor, "factor");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetReal (options_table, &minfactor, "minfactor");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetReal (options_table, &maxfactor, "maxfactor");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetReal (options_table, &acceleration, "acceleration");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetReal (options_table, &deceleration, "deceleration");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetRealArray (options_table,
+                                   nvars, &factors[0], "factors");
+    assert (ierr==nvars || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    
+    assert (maxiters >= 0);
+    assert (minerror >= 0);
+    assert (factor > 0);
+    assert (minfactor > 0 && minfactor <= factor);
+    assert (maxfactor >= factor);
+    assert (acceleration >= 1);
+    assert (deceleration > 0 && deceleration <= 1);
+    for (int n=0; n<nvars; ++n) {
+      assert (factors[n] != 0);
+    }
+    
+    
+    
+    assert (is_global_mode() || is_level_mode());
+    
+    // The level to solve on
+    const int solve_level = is_level_mode() ? reflevel : reflevels-1;
+    
+    if (solve_level < solve_minlevel) {
+      if (verbose || veryverbose) {
+        CCTK_VInfo (CCTK_THORNSTRING,
+                    "Did not solve because the level is too coarse.");
+      }
+      Util_TableSetInt (options_table, 0, "iters");
+      Util_TableSetReal (options_table, -1.0, "error");
+      
+      // Return as if the solving had not converged
+      return -1;
+    }
+    
+#warning "TODO: assert that all levels are at the same time"
+    
+    // Switch to global mode
+    BEGIN_GLOBAL_MODE(cctkGH) {
+      
+      // Reset the initial data time
+      global_time = 0;
+      delta_time = 1;
+      * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time;
+      * const_cast<CCTK_REAL *> (& cctkGH->cctk_delta_time) = delta_time;
+      BEGIN_REFLEVEL_LOOP(cctkGH) {
+        * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time;
+        for (int m=0; m<maps; ++m) {
+          vtt[m]->set_time (reflevel, mglevel, 0);
+          vtt[m]->set_delta
+            (reflevel, mglevel, 1.0 / ipow(reflevelfact, reflevelpower));
+        }
+      } END_REFLEVEL_LOOP;
+      
+      
+      
+      // Init statistics
+      vector<CCTK_REAL> norm_counts (solve_level+1);
+      vector<CCTK_REAL> norm_l2s (solve_level+1);
+      CCTK_REAL norm2 = HUGE_VAL;
+      CCTK_REAL speed = 0.0;
+      CCTK_REAL throttle = 1.0;
+      bool did_hit_min = false;
+      bool did_hit_max = false;
+      const time_t starttime = time(0);
+      time_t nexttime = starttime;
+      
+      if (verbose || veryverbose) {
+        CCTK_VInfo (CCTK_THORNSTRING,
+                    "Solving through adaptive Jacobi iterations");
+      }
+      
+      const int icoarsestep
+        = ipow (mglevelfact * maxreflevelfact, reflevelpower);
+      const int istep
+        = ipow (mglevelfact * maxreflevelfact / ipow(reffact, solve_level),
+                reflevelpower);
+      int iter = istep;
+      for (;; iter+=istep) {
+        
+        global_time = 1.0 * iter / ipow(maxreflevelfact, reflevelpower);
+        * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time;
+//         cout << "CJ global iter " << iter << " time " << global_time << flush << endl;
+        
+        // Calculate residual, smooth, and apply boundary conditions
+        ierr = 0;
+        BEGIN_REFLEVEL_LOOP(cctkGH) {
+          const int do_every
+            = ipow (mglevelfact * (maxreflevelfact/reflevelfact),
+                    reflevelpower);
+          if (reflevel <= solve_level && (iter-istep) % do_every == 0) {
+            
+            // Advance time
+            for (int m=0; m<maps; ++m) {
+              vtt[m]->advance_time (reflevel, mglevel);
+            }
+            * const_cast<CCTK_REAL *> (& cctkGH->cctk_time)
+              = (1.0 * (iter - istep + do_every)
+                 / ipow(maxreflevelfact, reflevelpower));
+            CycleTimeLevels (cctkGH);
+//             cout << "CJ residual iter " << iter << " reflevel " << reflevel << " time " << cctkGH->cctk_time << flush << endl;
+            
+            // Calculate residual
+            BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+              BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
+                if (! ierr) {
+                  ierr = calcres (cctkGH, options_table, userdata);
+                }
+              } END_LOCAL_COMPONENT_LOOP;
+            } END_MAP_LOOP;
+            
+            // Smooth and apply boundary conditions
+            CCTK_REAL levfac = 0;
+            for (int d=0; d<dim; ++d) {
+              levfac += 1 /
+                pow (cctkGH->cctk_delta_space[d] / cctkGH->cctk_levfac[d], 2);
+            }
+            levfac = 1 / levfac;
+            
+            BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+              BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
+                if (! ierr) {
+                  
+                  for (int n=0; n<nvars; ++n) {
+                    CCTK_REAL * restrict varptr
+                      = (CCTK_REAL *)CCTK_VarDataPtrI (cctkGH, 0, var[n]);
+                    assert (varptr);
+                    const CCTK_REAL * restrict resptr
+                      = (const CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, res[n]);
+                    assert (resptr);
+                    assert (resptr != varptr);
+                    const CCTK_REAL fac = factor * factors[n] * levfac;
+                    assert (cctkGH->cctk_dim == 3);
+                    for (int k=cctkGH->cctk_nghostzones[2];
+                         k<cctkGH->cctk_lsh[2]-cctkGH->cctk_nghostzones[2];
+                         ++k) {
+                      for (int j=cctkGH->cctk_nghostzones[1];
+                           j<cctkGH->cctk_lsh[1]-cctkGH->cctk_nghostzones[1];
+                           ++j) {
+                        for (int i=cctkGH->cctk_nghostzones[0];
+                             i<cctkGH->cctk_lsh[0]-cctkGH->cctk_nghostzones[0];
+                             ++i) {
+                          const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k);
+                          varptr[ind] += fac * resptr[ind];
+                        }
+                      }
+                    }
+                  } // for n
+                  
+                  // Apply boundary conditions
+                  ierr = applybnds (cctkGH, options_table, userdata);
+                  
+                }
+              } END_LOCAL_COMPONENT_LOOP;
+            } END_MAP_LOOP;
+            
+          }
+        } END_REFLEVEL_LOOP;
+        
+        if (ierr) {
+          if (verbose || veryverbose) {
+            CCTK_VInfo (CCTK_THORNSTRING,
+                        "Residual calculation or boundary conditions reported error; aborting solver.");
+          }
+          Util_TableSetInt (options_table, iter, "iters");
+          Util_TableDeleteKey (options_table, "error");
+          goto done;
+        }
+        
+        // Restrict and calculate norm
+        ierr = 0;
+        BEGIN_REVERSE_REFLEVEL_LOOP(cctkGH) {
+          const int do_every
+            = ipow (mglevelfact * (maxreflevelfact/reflevelfact),
+                    reflevelpower);
+          if (reflevel <= solve_level && iter % do_every == 0) {
+            
+//             cout << "CJ restrict iter " << iter << " reflevel " << reflevel << " time " << cctkGH->cctk_time << flush << endl;
+            
+            if (! ierr) {
+              if (reflevel < solve_level) {
+                Restrict (cctkGH);
+              }
+            }
+            
+            BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+              BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
+                if (! ierr) {
+                  ierr = applybnds (cctkGH, options_table, userdata);
+                }
+              } END_LOCAL_COMPONENT_LOOP;
+            } END_MAP_LOOP;
+            
+            // Initialise norm
+            CCTK_REAL norm_count = 0;
+            CCTK_REAL norm_l2 = 0;
+            
+            BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+              BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
+                if (! ierr) {
+                  
+                  // Calculate norm
+                  for (int n=0; n<nvars; ++n) {
+                    const CCTK_REAL * restrict resptr
+                      = (const CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, res[n]);
+                    assert (resptr);
+                    const CCTK_REAL fac = factors[n];
+                    assert (cctkGH->cctk_dim == 3);
+                    for (int k=cctkGH->cctk_nghostzones[2];
+                         k<cctkGH->cctk_lsh[2]-cctkGH->cctk_nghostzones[2];
+                         ++k) {
+                      for (int j=cctkGH->cctk_nghostzones[1];
+                           j<cctkGH->cctk_lsh[1]-cctkGH->cctk_nghostzones[1];
+                           ++j) {
+                        for (int i=cctkGH->cctk_nghostzones[0];
+                             i<cctkGH->cctk_lsh[0]-cctkGH->cctk_nghostzones[0];
+                             ++i) {
+                          const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k);
+                          ++norm_count;
+                          // TODO: scale the norm by the resolution?
+                          norm_l2 += pow(fac * resptr[ind], 2);
+                        }
+                      }
+                    }
+                  } // for n
+                  
+                }
+              } END_LOCAL_COMPONENT_LOOP;
+            } END_MAP_LOOP;
+            
+            const int sum_handle = CCTK_ReductionArrayHandle ("sum");
+            assert (sum_handle >= 0);
+            CCTK_REAL reduce_in[2], reduce_out[2];
+            reduce_in[0] = norm_count;
+            reduce_in[1] = norm_l2;
+            ierr = CCTK_ReduceLocArrayToArray1D
+              (cctkGH, -1, sum_handle, reduce_in, reduce_out, 2,
+               CCTK_VARIABLE_REAL);
+            norm_counts.at(reflevel) = reduce_out[0];
+            norm_l2s.at(reflevel) = reduce_out[1];
+            
+#warning "TODO"
+#if 0
+            CCTK_OutputVarAs (cctkGH, "WaveToyFO::phi", "phi-ell");
+            CCTK_OutputVarAs (cctkGH, "IDSWTEsimpleFO::residual", "residual-ell");
+#endif
+            
+          }
+        } END_REVERSE_REFLEVEL_LOOP;
+        
+        if (ierr) {
+          if (verbose || veryverbose) {
+            CCTK_VInfo (CCTK_THORNSTRING,
+                        "Residual calculation or boundary conditions reported error; aborting solver.");
+          }
+          Util_TableSetInt (options_table, iter, "iters");
+          Util_TableDeleteKey (options_table, "error");
+          goto done;
+        }
+        
+        // Save old norm
+        const CCTK_REAL old_norm2 = norm2;
+        
+        // Calculate new norm
+        CCTK_REAL norm_count = 0;
+        CCTK_REAL norm_l2 = 0;
+        for (int rl=0; rl<=solve_level; ++rl) {
+          norm_count += norm_counts[rl];
+          norm_l2 += norm_l2s[rl];
+        }
+        norm2 = sqrt(norm_l2 / norm_count);
+        
+        // Calculate convergence speed
+        speed = old_norm2 / norm2;
+        
+        // Log output
+        if (verbose || veryverbose) {
+          const time_t currenttime = time(0);
+          if ((iter % icoarsestep == 0 && iter >= maxiters)
+#if HAVE_ISNAN
+              || isnan(norm2)
+#endif
+              || norm2 <= minerror
+              || (iter % outevery == 0 && currenttime >= nexttime)) {
+            if (verbose || veryverbose) {
+              if (did_hit_min) {
+                CCTK_VInfo (CCTK_THORNSTRING,
+                            "Convergence factor became too small: artificially kept at minfactor (may lead to instability)");
+              }
+              if (did_hit_max) {
+                CCTK_VInfo (CCTK_THORNSTRING,
+                            "Convergence factor became too large: artificially kept at maxfactor (may lead to slower convergence)");
+              }
+              did_hit_min = false;
+              did_hit_max = false;
+            }
+            if (veryverbose) {
+              CCTK_VInfo (CCTK_THORNSTRING,
+                          "Iteration %d, time %g, factor %g, throttle %g, residual %g, speed %g",
+                          iter, double(currenttime-starttime),
+                          double(factor), double(throttle), double(norm2),
+                          double(speed));
+            } else {
+              CCTK_VInfo (CCTK_THORNSTRING,
+                          "Iteration %d, time %g, residual %g",
+                          iter, double(currenttime-starttime), double(norm2));
+            }
+            if (outeveryseconds > 0) {
+              while (nexttime <= currenttime) nexttime += outeveryseconds;
+            }
+          }
+        }
+        
+#if HAVE_ISNAN
+        // Exit if things go bad
+        if (isnan(norm2)) {
+          if (verbose || veryverbose) {
+            CCTK_VInfo (CCTK_THORNSTRING,
+                        "Encountered NaN.  Aborting solver.");
+          }
+          break;
+        }
+#endif
+        
+        // Return when desired accuracy is reached
+        if (norm2 <= minerror) {
+          if (verbose || veryverbose) {
+            CCTK_VInfo (CCTK_THORNSTRING, "Done solving.");
+          }
+          Util_TableSetInt (options_table, iter, "iters");
+          Util_TableSetReal (options_table, norm2, "error");
+          ierr = 0;
+          goto done;
+        }
+        
+        // Exit if the maximum number of iterations has been reached
+        if (iter % icoarsestep == 0 && iter >= maxiters) {
+          break;
+        }
+        
+        // Adjust speed
+        if (speed > 1.0) {
+          throttle = acceleration;
+        } else {
+          throttle = deceleration;
+        }
+        factor *= throttle;
+        if (factor < minfactor) {
+          did_hit_min = true;
+          factor = minfactor;
+        }
+        if (factor > maxfactor) {
+          did_hit_max = true;
+          factor = maxfactor;
+        }
+        
+      } // for iter
+      
+      // Did not solve
+      if (verbose || veryverbose) {
+        CCTK_VInfo (CCTK_THORNSTRING, "Did not converge.");
+      }
+      Util_TableSetInt (options_table, iter, "iters");
+      Util_TableSetReal (options_table, norm2, "error");
+      ierr = -1;
+    
+    done: ;
+      
+      
+      
+      // Reset the initial time
+#warning "TODO: reset the initial time a bit more intelligently"
+      global_time = 0;
+      delta_time = 1;
+      * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time;
+      * const_cast<CCTK_REAL *> (& cctkGH->cctk_delta_time) = delta_time;
+      BEGIN_REFLEVEL_LOOP(cctkGH) {
+        * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time;
+        for (int m=0; m<maps; ++m) {
+          vtt[m]->set_time (reflevel, mglevel, 0);
+          vtt[m]->set_delta (reflevel, mglevel, 1.0 / reflevelfact);
+        }
+      } END_REFLEVEL_LOOP;
+      
+    } END_GLOBAL_MODE;
+    
+    return ierr;
+  }
+
+} // namespace CarpetJacobi
diff --git a/CarpetDev/CarpetJacobi/src/make.code.defn b/CarpetDev/CarpetJacobi/src/make.code.defn
new file mode 100644
index 000000000..efd1ad4c4
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/src/make.code.defn
@@ -0,0 +1,9 @@
+# Main make.code.defn file for thorn CarpetJacobi
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/src/make.code.defn,v 1.1 2003/09/02 14:35:58 schnetter Exp $
+
+# Source files in this directory
+SRCS = Jacobi.cc
+
+# Subdirectories containing source files
+SUBDIRS = 
+