Commits

Matt Knepley committed 4f423c0 Merge

Comments (0)

Files changed (14)

config/linux-gnu-ia64

+#! /bin/csh
+set DATAFILESPATH=/home/balay/datafiles
+

config/linux-gnu-ia64-intel

+#! /bin/csh
+set DATAFILESPATH=/home/balay/datafiles
+

python/PETSc/packages/blopex.py

 class Configure(PETSc.package.Package):
   def __init__(self, framework):
     PETSc.package.Package.__init__(self, framework)
-    self.download   = ['ftp://ftp.mcs.anl.gov/pub/petsc/externalpackages/blopex_abstract_Nov_2005.tar.gz']
+    self.download   = ['ftp://ftp.mcs.anl.gov/pub/petsc/externalpackages/blopex_abstract_Aug_2006.tar.gz']
     self.functions  = ['lobpcg_solve']
     self.includes   = ['interpreter.h']
     self.liblist    = [['libBLOPEX.a']]

src/contrib/blopex/driver/driver.c

-/* Program usage:  mpirun -np <procs> driver [-help] [all PETSc options] */ 
-
+/* This code was developed by Merico Argentati, Andrew Knyazev, Ilya Lashuk and Evgueni Ovtchinnikov */
 
+/* Program usage:  mpirun -np <procs> driver [-help] [all PETSc options] */ 
 
 static char help[] = "Test driver for 'abstract lobpcg' in PETSC\n\
 Usage: mpirun -np <procs> driver [-help] [all PETSc options]\n\
 #include "petscksp.h"
 #include "petscda.h"
 #include <assert.h>
-#include "lobpcg.h" 
-#include "src/contrib/blopex/petsc-interface/petsc-interface.h" 
+#include "lobpcg.h"
+#include "src/contrib/blopex/petsc-interface/petsc-interface.h"
 #include "interpreter.h"
 #include "multivector.h"
 
 {
       PetscErrorCode     ierr;
       
-      ierr = KSPSolve(((aux_data_struct*)data)->ksp, x, y);
+      ierr = KSPSolve(((aux_data_struct*)data)->ksp, (Vec)x, (Vec)y);
       assert(!ierr);
 }
 

src/contrib/blopex/driver_fiedler/driver_fiedler.c

+/* This code was developed by Merico Argentati, Andrew Knyazev, Ilya Lashuk and Evgueni Ovtchinnikov */
+
+static char help[] = "'Fiedler' test driver for 'abstract lobpcg' in PETSC\n\
+Usage: mpirun -np <procs> driver_fiedler [-help] [all PETSc options]\n\
+Special options:\n\
+-matrix <filename>      (mandatory) specify file with 'stiffness' matrix \
+(in petsc format)\n\
+-mass_matrix <filename> (optional) 'mass' matrix for generalized eigenproblem\n\
+-n_eigs <integer>      Number of eigenvalues to calculate\n\
+-tol <real number>     absolute tolerance for residuals\n\
+-full_out              Produce more output\n\
+-seed <integer>        seed for random number generator\n\
+-itr <integer>         Maximal number of iterations\n\
+-output_file <string>  Filename to write calculated eigenvectors.\n\
+-shift <real number>   Apply shift to 'stiffness' matrix\n\
+Example:\n\
+mpirun -np 2 ./driver_fiedler -matrix my_matrix.bin -n_eigs 3 -tol 1e-6 -itr 20\n";
+
+#include "petscksp.h"
+#include <assert.h>
+#include "lobpcg.h" 
+#include "src/contrib/blopex/petsc-interface/petsc-interface.h" 
+#include "interpreter.h"
+#include "multivector.h"
+#include "temp_multivector.h"
+
+typedef struct
+{
+  KSP                      ksp;
+  Mat                      A;
+  Mat                      B;
+  mv_InterfaceInterpreter  ii;
+} aux_data_struct;
+
+void Precond_FnSingleVector(void * data, void * x, void * y)
+{
+      PetscErrorCode     ierr;
+      
+      ierr = KSPSolve(((aux_data_struct*)data)->ksp, (Vec)x, (Vec)y);
+      assert(!ierr);
+}
+
+void Precond_FnMultiVector(void * data, void * x, void * y)
+{
+   ((aux_data_struct*)data)->ii.Eval(Precond_FnSingleVector, data, x, y);
+}
+
+void OperatorASingleVector (void * data, void * x, void * y)
+{
+   PetscErrorCode     ierr;
+   
+   ierr = MatMult(((aux_data_struct*)data)->A, (Vec)x, (Vec)y);
+   assert(!ierr);
+}
+
+void OperatorAMultiVector(void * data, void * x, void * y)
+{
+   ((aux_data_struct*)data)->ii.Eval(OperatorASingleVector, data, x, y);
+}
+
+
+void OperatorBSingleVector (void * data, void * x, void * y)
+{
+   PetscErrorCode     ierr;
+   
+   ierr = MatMult(((aux_data_struct*)data)->B, (Vec)x, (Vec)y);
+   assert(!ierr);
+}
+
+void OperatorBMultiVector(void * data, void * x, void * y)
+{
+   ((aux_data_struct*)data)->ii.Eval(OperatorBSingleVector, data, x, y);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc,char **args)
+{
+   Vec            u;
+   Vec            tmp_vec;
+   Mat            A;
+   Mat            B;
+
+   PetscErrorCode ierr;
+  
+   mv_MultiVectorPtr          eigenvectors;   
+   mv_MultiVectorPtr          constraints;
+   mv_TempMultiVector*        raw_constraints;
+   mv_TempMultiVector*        raw_eigenvectors;
+   
+   double *                   eigs;
+   double *                   eigs_hist;
+   double *                   resid;
+   double *                   resid_hist;
+   int                        iterations;
+   PetscMPIInt                rank;
+   int                        n_eigs = 1;
+   int                        seed = 1;
+   int                        i,j;
+   PetscLogDouble             t1,t2,elapsed_time;
+   double                     tol=1e-08;
+   PetscTruth                 full_output=PETSC_FALSE;
+   KSP                        ksp;
+   lobpcg_Tolerance           lobpcg_tol;
+   int                        maxIt = 100;
+   mv_InterfaceInterpreter    ii;
+   lobpcg_BLASLAPACKFunctions blap_fn;
+   aux_data_struct            aux_data;
+   PetscViewer                fd;               /* viewer */
+   PetscTruth                 matrix_present;
+   char                       filename[PETSC_MAX_PATH_LEN];
+   char                       mass_filename[PETSC_MAX_PATH_LEN];
+   PetscTruth                 mass_matrix_present;
+   PetscReal                  shift;
+   PetscTruth                 shift_present;
+   char                       output_filename[PETSC_MAX_PATH_LEN];
+   PetscTruth                 output_filename_present;
+   char                       tmp_str[PETSC_MAX_PATH_LEN];
+   
+   PetscInitialize(&argc,&args,(char *)0,help);
+   
+   /* read command-line parameters */
+   ierr = PetscOptionsGetInt(PETSC_NULL,"-n_eigs",&n_eigs,PETSC_NULL);CHKERRQ(ierr);
+   ierr = PetscOptionsGetReal(PETSC_NULL,"-tol", &tol,PETSC_NULL); CHKERRQ(ierr);
+   ierr = PetscOptionsGetString(PETSC_NULL,"-matrix",filename,PETSC_MAX_PATH_LEN-1,
+           &matrix_present);
+   CHKERRQ(ierr);
+   if (!matrix_present) 
+    SETERRQ(1,"Must indicate binary file to read matrix from with the "
+            "'-matrix' option");
+   ierr = PetscOptionsGetString(PETSC_NULL,"-mass_matrix",mass_filename,
+           PETSC_MAX_PATH_LEN-1,&mass_matrix_present);
+   CHKERRQ(ierr);
+   ierr = PetscOptionsHasName(PETSC_NULL,"-full_out",&full_output); CHKERRQ(ierr);
+   ierr = PetscOptionsGetInt(PETSC_NULL,"-seed",&seed,PETSC_NULL);CHKERRQ(ierr);
+   if (seed<1)
+    seed=1;
+   ierr = PetscOptionsGetInt(PETSC_NULL,"-itr",&maxIt,PETSC_NULL);CHKERRQ(ierr);
+   ierr = PetscOptionsGetReal(PETSC_NULL,"-shift",&shift,&shift_present);
+   ierr = PetscOptionsGetString(PETSC_NULL,"-output_file",output_filename,
+            PETSC_MAX_PATH_LEN-1, &output_filename_present);
+           
+           
+   /* load matrices */
+   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd);
+   CHKERRQ(ierr);
+   ierr = MatLoad(fd,MATMPIAIJ,&A);CHKERRQ(ierr);
+   ierr = PetscViewerDestroy(fd);CHKERRQ(ierr);
+   
+   if (mass_matrix_present)
+   {
+       ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,mass_filename,FILE_MODE_READ,&fd);
+       CHKERRQ(ierr);
+       ierr = MatLoad(fd,MATMPIAIJ,&B);CHKERRQ(ierr);
+       ierr = PetscViewerDestroy(fd);CHKERRQ(ierr);
+   }
+     
+   /* apply shift to stiffness matrix if asked to do so */
+   if (shift_present)
+   {
+      if (mass_matrix_present)
+      {
+        ierr = MatAXPY(A,shift,B,SUBSET_NONZERO_PATTERN);
+        CHKERRQ(ierr);
+      }
+      else
+      {
+        ierr = MatShift(A,shift);
+        CHKERRQ(ierr);
+      }
+   }
+        
+   
+   /* 
+    Create parallel vectors.
+     - We form 1 vector from scratch and then duplicate as needed.
+   */
+
+   MatGetVecs(A,&u,NULL);
+
+   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
+               Create the linear solver and set various options
+    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+ /* Here we START measuring time for preconditioner setup */
+   ierr = PetscGetTime(&t1);CHKERRQ(ierr);
+
+   /* 
+      Create linear solver context
+   */
+   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
+
+   /* 
+      Set operators. Here the matrix that defines the linear system
+      also serves as the preconditioning matrix.
+   */
+   ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
+
+   /* 
+     Set runtime options, e.g.,
+         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
+     These options will override those specified above as long as
+     KSPSetFromOptions() is called _after_ any other customization
+     routines.
+   */
+   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+
+   /* probably this call actually builds the preconditioner */
+   ierr = KSPSetUp(ksp);CHKERRQ(ierr);
+
+   /* Here we STOP measuring time for preconditioner setup */
+   ierr = PetscGetTime(&t2);CHKERRQ(ierr);
+   elapsed_time=t2-t1;
+
+   PetscPrintf(PETSC_COMM_WORLD,"Preconditioner setup, seconds: %f\n",elapsed_time);
+
+   /* request memory for eig-vals */
+   ierr = PetscMalloc(sizeof(double)*n_eigs,&eigs); CHKERRQ(ierr);
+
+   /* request memory for eig-vals history */
+   ierr = PetscMalloc(sizeof(double)*n_eigs*(maxIt+1),&eigs_hist); CHKERRQ(ierr);
+
+   /* request memory for resid. norms */
+   ierr = PetscMalloc(sizeof(double)*n_eigs,&resid); CHKERRQ(ierr);
+
+   /* request memory for resid. norms hist. */
+   ierr = PetscMalloc(sizeof(double)*n_eigs*(maxIt+1),&resid_hist); CHKERRQ(ierr);
+
+   LOBPCG_InitRandomContext();
+
+   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
+
+   PETSCSetupInterpreter( &ii );
+   eigenvectors = mv_MultiVectorCreateFromSampleVector(&ii, n_eigs,u);
+
+   /* building constraints (constant vector */
+   constraints= mv_MultiVectorCreateFromSampleVector(&ii, 1,u);
+   raw_constraints = (mv_TempMultiVector*)mv_MultiVectorGetData (constraints);
+   tmp_vec = (Vec)(raw_constraints->vector)[0];
+   ierr = VecSet(tmp_vec,1.0); CHKERRQ(ierr); 
+
+   for (i=0; i<seed; i++) /* this cycle is to imitate changing random seed */
+      mv_MultiVectorSetRandom (eigenvectors, 1234);
+
+   lobpcg_tol.absolute = tol;
+   lobpcg_tol.relative = 1e-50;
+
+   blap_fn.dpotrf = PETSC_dpotrf_interface;
+   blap_fn.dsygv = PETSC_dsygv_interface;
+
+   aux_data.A = A;
+   aux_data.B = B;
+   aux_data.ksp = ksp;
+   aux_data.ii = ii;
+
+   /* Here we START measuring time for solution process */
+   ierr = PetscGetTime(&t1);CHKERRQ(ierr);
+
+   lobpcg_solve
+   ( 
+      eigenvectors,
+      &aux_data,
+      OperatorAMultiVector,
+      mass_matrix_present?&aux_data:NULL,
+      mass_matrix_present?OperatorBMultiVector:NULL,
+      &aux_data,
+      Precond_FnMultiVector,
+      constraints,
+      blap_fn,
+      lobpcg_tol,
+      maxIt,
+      !rank,
+      &iterations,
+	    eigs,
+      eigs_hist, 
+      n_eigs,                  	      
+      resid, 
+      resid_hist,
+      n_eigs                                     
+   );
+
+   /* Here we STOP measuring time for solution process */
+   ierr = PetscGetTime(&t2);CHKERRQ(ierr);
+   elapsed_time=t2-t1;
+
+   PetscPrintf(PETSC_COMM_WORLD,"Solution process, seconds: %e\n",elapsed_time);
+   
+   /* shift eigenvalues back */
+   for (i=0; i<n_eigs; i++)
+      eigs[i]-=shift;
+      
+   PetscPrintf(PETSC_COMM_WORLD,"Final eigenvalues:\n");     
+   for (i=0;i<n_eigs;i++)
+	 {
+		 ierr = PetscPrintf(PETSC_COMM_WORLD,"%e\n",eigs[i]);
+		 CHKERRQ(ierr);
+	 }
+   
+   if (full_output)
+   {
+     PetscPrintf(PETSC_COMM_WORLD,"Output from LOBPCG, eigenvalues history:\n");
+     for (j=0; j<iterations+1; j++)
+       for (i=0;i<n_eigs;i++)
+       {
+         ierr = PetscPrintf(PETSC_COMM_WORLD,"%e\n",*(eigs_hist+j*n_eigs+i));
+         CHKERRQ(ierr);
+	     }
+     PetscPrintf(PETSC_COMM_WORLD,"Output from LOBPCG, residual norms:\n");
+     for (i=0;i<n_eigs;i++)
+	   {
+		   ierr = PetscPrintf(PETSC_COMM_WORLD,"%e\n",resid[i]);
+		   CHKERRQ(ierr);
+	   }
+
+     PetscPrintf(PETSC_COMM_WORLD,"Output from LOBPCG, residual norms history:\n");
+     for (j=0; j<iterations+1; j++)
+       for (i=0;i<n_eigs;i++)
+       {
+         ierr = PetscPrintf(PETSC_COMM_WORLD,"%e\n",*(resid_hist+j*n_eigs+i));
+         CHKERRQ(ierr);
+       }
+   }	
+
+   /* write eigenvectors to disk, if told to do so */
+   
+   if (output_filename_present)
+   {
+      raw_eigenvectors = (mv_TempMultiVector*)mv_MultiVectorGetData (eigenvectors);
+      tmp_vec = (Vec)(raw_constraints->vector)[0];
+      for ( i = 0; i < n_eigs; i++ ) 
+      {
+        sprintf( tmp_str, "%s_%d.petsc", output_filename, i );
+        PetscViewerBinaryOpen(PETSC_COMM_WORLD, tmp_str, FILE_MODE_WRITE, &fd);
+        /* PetscViewerSetFormat(fd,PETSC_VIEWER_ASCII_MATLAB); */
+        ierr = VecView((Vec)(raw_eigenvectors->vector)[i],fd); CHKERRQ(ierr); 
+        ierr = PetscViewerDestroy(fd);CHKERRQ(ierr);
+      }
+  
+   }
+   
+   /*
+      Free work space.  All PETSc objects should be destroyed when they
+      are no longer needed.
+   */
+   ierr = VecDestroy(u);CHKERRQ(ierr);
+   ierr = MatDestroy(A);CHKERRQ(ierr);
+   if (mass_matrix_present)
+     ierr = MatDestroy(B);CHKERRQ(ierr);
+   ierr = KSPDestroy(ksp);CHKERRQ(ierr);
+
+   LOBPCG_DestroyRandomContext();
+   mv_MultiVectorDestroy(eigenvectors);
+   mv_MultiVectorDestroy(constraints);
+
+   /* free memory used for eig-vals */
+   ierr = PetscFree(eigs); 
+   CHKERRQ(ierr);
+   ierr = PetscFree(eigs_hist); 
+   CHKERRQ(ierr);
+   ierr = PetscFree(resid); 
+   CHKERRQ(ierr);
+   ierr = PetscFree(resid_hist); 
+   CHKERRQ(ierr);
+
+   ierr = PetscFinalize();CHKERRQ(ierr);
+   return 0;
+}

src/contrib/blopex/driver_fiedler/makefile

+
+CFLAGS           = ${BLOPEX_INCLUDE} 
+FFLAGS           =
+CPPFLAGS         =
+FPPFLAGS         =
+LOCDIR           = src/contrib/blopex/driver_fiedler
+EXAMPLESC        = driver.c
+EXAMPLESF        = 
+MANSEC           = KSP
+
+
+include ${PETSC_DIR}/bmake/common/base
+
+driver_fiedler: driver_fiedler.o  chkopts
+	-${CLINKER} -o driver_fiedler driver_fiedler.o  ${PETSC_KSP_LIB}
+	${RM} driver_fiedler.o
+
+

src/contrib/blopex/driver_fiedler/readme

+The file "driver_fiedler.c" contains a driver for spectral segmentation problems. The driver solves 
+partial generalized eigenvalue problem with constant vector as a constraint.
+
+Usage: mpirun -np <procs> driver_fiedler [-help] [driver options] [all PETSc options]
+Driver options:
+-matrix <filename>      (mandatory) Specify file with 'stiffness' matrix. Matrix should be stored in petsc 
+                        binary format.
+
+-mass_matrix <filename> (optional) 'mass' matrix for generalized eigenproblem. Assumed identity if omitted. 
+                        If present, and if shift (see below) is present, then nonzero pattern of 'mass'
+			matrix must be a subset of nonzero pattern of 'stiffness' matrix.
+
+-n_eigs <integer>       Number of (smallest) eigenvalues to calculate.
+
+-tol <real number>     absolute tolerance for residuals
+
+-full_out              Produce more output
+
+-seed <integer>        seed for random number generator
+
+-itr <integer>         Maximal number of iterations
+
+-output_file <string>  Filename (prefix) to write calculated eigenvectors. For example, if the value of this 
+                       parameter is 'vectors', and 3 eigenvectors are calculated, then 3 files will be created:
+		       'vectors_0.petsc', 'vectors_1.petsc', and 'vectors_2.petsc'. If this parameter is 
+		       omitted, no files are created.
+
+-shift <real number>   Apply (add) shift to 'stiffness' matrix. If 'mass' matrix is present, then 'mass' matrix
+                       multiplied by shift is added to 'stiffness' matrix. The shift is subsequently subtracted 
+		       from 'final' eigenvalues. This option may be useful if for shifted matrix some 
+		       particular preconditioning technique is applicable.
+Example:
+mpirun -np 2 driver -matrix my_matrix.bin -n_eigs 3 -tol 1e-6 -itr 20

src/contrib/blopex/petsc-interface/makefile

 LIBBASE  = libpetscksp
 DIRS     =
 MANSEC   = 
-LOCDIR   = src/contrib/blopex/petsc-interface
+LOCDIR   = src/contrib/blopex/petsc-interface/
 
 include ${PETSC_DIR}/bmake/common/base
 include ${PETSC_DIR}/bmake/common/test

src/contrib/blopex/petsc-interface/petsc-interface.c

+/* This code was developed by Merico Argentati, Andrew Knyazev, Ilya Lashuk and Evgueni Ovtchinnikov */
+
 #include "petscsys.h"
 #include "petscvec.h"
 #include "petscmat.h"
 int
 LOBPCG_InitRandomContext(void)
 {
-  PetscErrorCode ierr;
+	PetscErrorCode ierr;
   /* PetscScalar rnd_bound = 1.0; */
   
-  ierr = PetscRandomCreate(PETSC_COMM_WORLD,&LOBPCG_RandomContext);CHKERRQ (ierr);
-  ierr = PetscRandomSetFromOptions(LOBPCG_RandomContext);CHKERRQ(ierr);
+  ierr = PetscRandomCreate(PETSC_COMM_WORLD,&LOBPCG_RandomContext);CHKERRQ(ierr);
+  ierr = PetscRandomSetFromOptions(LOBPCG_RandomContext);CHKERRQ(ierr); 
   
-  ierr = PetscRandomSetInterval(LOBPCG_RandomContext,(PetscScalar)-1.0,(PetscScalar)1.0);CHKERRQ(ierr);
-  return 0;
+  ierr = PetscRandomSetInterval(LOBPCG_RandomContext,(PetscScalar)-1.0,(PetscScalar)1.0);
+	CHKERRQ(ierr);
+	return 0;
 }
 
 int 

src/contrib/blopex/petsc-interface/petsc-interface.h

+/* This code was developed by Merico Argentati, Andrew Knyazev, Ilya Lashuk and Evgueni Ovtchinnikov */
+
 #ifndef PETSC_INTERFACE_HEADER
 #define PETSC_INTERFACE_HEADER
 
 #include "interpreter.h"
 
-#ifdef __cplusplus
-extern "C" {
-#endif
-
 int PETSC_dpotrf_interface (char *uplo, int *n, double *a, int * lda, int *info);
 
 int PETSC_dsygv_interface (int *itype, char *jobz, char *uplo, int *
 int
 PETSCSetupInterpreter( mv_InterfaceInterpreter *ii );
 
-#ifdef __cplusplus
-}
-#endif
-
 #endif /* PETSC_INTERFACE_HEADER */

src/dm/mesh/examples/tests/xsifter0.cxx

 
 #undef __FUNCT__
 #define __FUNCT__ "BasicBaseTest"
-PetscErrorCode BasicBaseTest(const ALE_X::Obj<sifter_type>& sifter, Options options)
+PetscErrorCode BasicBaseTest(const ALE::Obj<sifter_type>& sifter, ALE_X::Test::Options options)
 {
-  ALE_X::Obj<sifter_type::BaseSequence> base = sifter.base();
+  ALE::Obj<sifter_type::BaseSequence> base = sifter->base();
 
   PetscFunctionBegin;
+  ALE::LogStage stage = ALE::LogStageRegister("Base Test");
   ALE::LogStagePush(stage);
   std::cout << "Basic base:" << std::endl;
   sifter_type::BaseSequence::iterator begin, end, itor;
-  begin = base.begin();
-  end   = base.end();
+  begin = base->begin();
+  end   = base->end();
   itor = begin;
   std::cout << *itor;
   for(; itor != end; ++itor) {
 #define __FUNCT__ "main"
 int main(int argc, char *argv[])
 {
-  Options        options;
+  ALE_X::Test::Options        options;
   PetscErrorCode ierr;
 
   PetscFunctionBegin;

src/dm/mesh/sieve/XSifter.hh

 
 #include <iostream>
 
-// ALE extensions
 
 #ifndef  included_ALE_hh
 #include <ALE.hh>
       typedef lex1<key_type, key_order_type>                                   order1_type;
       typedef lex2<key_type, rec_type, key_order_type, xxx_order_type>         order2_type;
     private:
+      order1_type        _order1;
+      order2_type        _order2;
+      key_extractor_type _kex;
     public:
       bool operator()(const rec_type& rec1, const rec_type& rec2) { 
-        static order2_type        _order2;
-        static key_extractor_type _kex;
-        return _order2(_kex(rec1),rec1,_kex(rec2),rec2);
+        return this->_order2(this->_kex(rec1),rec1,this->_kex(rec2),rec2);
       };
       template <typename CompatibleKey_>
       bool operator()(const rec_type& rec1, const CompatibleKey_& key) {
         // We want key to be less than any (key, ...)
-        return !_compare1(key,_kex(rec1));
+        return !this->_order1(key,this->_kex(rec1));
       };
       template <typename CompatibleKey_, typename CompatibleXXXKey_>
       bool operator()(const rec_type& rec1, const ALE::pair<CompatibleKey_, CompatibleXXXKey_>& keyPair) {
         // We want (key,xxxkey) to be less than any (key, xxxkey, ...)
-        return !_compare2(keyPair.first,keyPair.second,_kex(rec1),rec1);
+        return !this->_order2(keyPair.first,keyPair.second,this->_kex(rec1),rec1);
       };
     };// class RecKeyXXXOrder
 
       typedef typename value_extractor_type::result_type       value_type;
       typedef typename index_type::size_type                   size_type;
       typedef typename index_type::iterator                    itor_type;
+      typedef typename index_type::reverse_iterator            ritor_type;
       //
       class iterator {
       public:
         /* value_type defined in the containing StridedIndexSequence */
       protected:
         // Parent sequence
-        sequence_type&  _sequence;
+        sequence_type  *_sequence;
         // Underlying iterator & segment boundary
         itor_type       _itor;
         itor_type       _segBndry;
         inner_key_extractor_type _ikex;
         value_extractor_type     _ex;
       public:
-        iterator(sequence_type& sequence, const itor_type& itor, const itor_type& segBndry) : 
+        iterator() : _sequence(NULL) {};
+        iterator(sequence_type *sequence, const itor_type& itor, const itor_type& segBndry) : 
           _sequence(sequence), _itor(itor), _segBndry(segBndry) {};
-        iterator(const iterator& iter):_sequence(iter._sequence), _itor(iter._itor),_segBndry(iter.segBndry){};
+        iterator(const iterator& iter):_sequence(iter._sequence), _itor(iter._itor),_segBndry(iter._segBndry){};
         virtual ~iterator() {};
         virtual bool              operator==(const iterator& iter) const {return this->_itor == iter._itor;};
         virtual bool              operator!=(const iterator& iter) const {return this->_itor != iter._itor;};
         // FIX: operator*() should return a const reference, but it won't compile that way, because _ex() returns const value_type
         virtual const value_type  operator*() const {return _ex(*(this->_itor));};
         virtual iterator   operator++() {
-          this->_sequence.next(this->_itor, this->_segBndry, inner_strided_flag);
+          this->_sequence->next(this->_itor, this->_segBndry, inner_strided_flag);
           return *this;
         };
         virtual iterator   operator++(int n) {iterator tmp(*this); ++(*this); return tmp;};
       };// class iterator
     protected:
-      index_type&     _index;
+      index_type      *_index;
+      //
       outer_key_type  _ohigh, _olow;
       bool            _have_olow, _have_ohigh;
       outer_key_type  _ihigh, _ilow;
       //
       // Basic interface
       //
+      StridedIndexSequence() : _index(NULL) {};
       StridedIndexSequence(const StridedIndexSequence& seq) : _index(seq._index), _olow(seq._olow), _ohigh(seq._ohigh), _have_olow(seq._have_olow), _have_ohigh(seq._have_ohigh), _ilow(seq._ilow), _ihigh(seq._ihigh), _have_ilow(seq._have_ilow), _have_ihigh(seq._have_ihigh)
       {};
-      StridedIndexSequence(index_type& index)  :  _index(index) {
+      StridedIndexSequence(index_type *index)  :  _index(index) {
         this->_have_olow = false; this->_have_ohigh = false;
         this->_have_ilow = false; this->_have_ihigh = false;
       };
       StridedIndexSequence& operator=(const StridedIndexSequence& seq) {
         copy(seq,*this); return *this;
       };
-      void reset(index_type& index) {
+      void reset(index_type *index) {
         this->_index      = index;
         this->_have_olow  = false;
         this->_have_ohigh = false;
       iterator begin() {
         static itor_type itor;
         // Determine the lower outer limit iterator
-        if(this->_have_olow()) {
+        if(this->_have_olow) {
           // ASSUMPTION: index ordering operator can compare against outer_keys
-          this->_itor = this->_index.lower_bound(this->_olow);
+          itor = this->_index->lower_bound(this->_olow);
         }
         else {
-          this->_itor = this->_index.begin();
+          itor = this->_index->begin();
         }
         // Now determine the inner lower limit and set the iterator to that limit within the first segment
         if(this->_have_ilow) {
           // ASSUMPTION: index ordering operator can compare against (outer_key, inner_key) pairs
-          itor = this->_index.lower_bound(ALE::pair<outer_key_type, inner_key_type>(this->_okex(*itor),this->_ilow));
+          itor = this->_index->lower_bound(ALE::pair<outer_key_type, inner_key_type>(this->_okex(*itor),this->_ilow));
         }
         else {
           // the itor is already in the right place: nothing to do
         }  
         // ASSUMPTION: index ordering operator can compare against (outer_key, inner_key) pairs
         static itor_type segBndry;
-        segBndry = this->_index.upper_bound(ALE::pair<outer_key_type, inner_key_type>(this->_okex(*itor),_ikex(*itor)));
-        return iterator(*this, itor, segBndry);
+        segBndry = this->_index->upper_bound(ALE::pair<outer_key_type, inner_key_type>(this->_okex(*itor),this->_ikex(*itor)));
+        return iterator(this, itor, segBndry);
       }; // begin()
       //
       void next(itor_type& itor, itor_type& segBndry, bool inner_strided = false) {
           // ASSUMPTION: index ordering operator can compare against (outer_key, inner_key) pairs
           olow = this->_okex(*itor);
           ilow = this->_ikex(*itor);
-          segBndry = this->_index.upper_bound(ALE::pair<outer_key_type, inner_key_type>(olow,ilow));
+          segBndry = this->_index->upper_bound(ALE::pair<outer_key_type, inner_key_type>(olow,ilow));
         }// inner strided
         // Otherwise, we iterate *within* a segment until its end is reached; then the following segment is started.
         else {
             olow = this->_okex(*itor);
             // Compute the lower boundary of the new segment
             // ASSUMPTION: index ordering operator can compare against outer_keys
-            itor = this->_index.upper_bound(olow);
+            itor = this->_index->upper_bound(olow);
             // Extract the new outer key
             olow = this->_okex(*itor);
             // Now determine the inner lower limit and set the iterator to that limit within the new segment
             if(this->_have_ilow) {
               ilow = this->_ilow;
               // ASSUMPTION: index ordering operator can compare against (outer_key, inner_key) pairs
-              itor = this->_index.lower_bound(ALE::pair<outer_key_type, inner_key_type>(olow,ilow));
+              itor = this->_index->lower_bound(ALE::pair<outer_key_type, inner_key_type>(olow,ilow));
             }
             else {
               // the itor is already in the right place; need to extract the ilow key
             }
             // Finally, compute the new segment's boundary
             // ASSUMPTION: index ordering operator can compare against (outer_key, inner_key) pairs
-            segBndry = this->_index.upper_bound(ALE::pair<outer_key_type, inner_key_type>(olow,ilow));
+            segBndry = this->_index->upper_bound(ALE::pair<outer_key_type, inner_key_type>(olow,ilow));
           }
         }// inner not strided
       };// next()
       //
       iterator end() {
-        itor_type itor;
+        static itor_type itor;
+        static ritor_type ritor;
         // Determine the upper outer limit
-        outer_key_type ohigh;
-        if(!this->_have_ohigh) {
-          itor = this->_index.rbegin();
-          ohigh = this->_okex(*itor);
+        static outer_key_type ohigh;
+        if(this->_have_ohigh) {
+          ohigh = this->_ohigh;
+        }
+        else {
+          ritor = this->_index->rbegin();
+          ohigh = this->_okex(*ritor);
         }
         // Determine the inner outer limit
-        inner_key_type ihigh;
+        static inner_key_type ihigh;
         if(this->_have_ihigh) {
           ihigh = this->_ihigh;
           // ASSUMPTION: index ordering operator can compare against (outer_key, inner_key) pairs
-          itor = this->_index.upper_bound(ALE::pair<outer_key_type, inner_key_type>(ohigh,ihigh));
+          itor = this->_index->upper_bound(ALE::pair<outer_key_type, inner_key_type>(ohigh,ihigh));
         }
         else {
           // ASSUMPTION: index ordering operator can compare against outer_keys
-          itor = this->_index.upper_bound(ohigh);
+          itor = this->_index->upper_bound(ohigh);
         }
         // use segBndry == itor
-        return iterator(*this, itor, itor); 
+        return iterator(this, itor, itor); 
       };// end()
       //
+      virtual bool contains(const outer_key_type& ok, const inner_key_type& ik) {
+        // FIX: This has to be implemented correctly, using the index ordering operator.
+        //return (this->_index->find(ALE::pair<outer_key_type,inner_key_type>(ok,ik)) != this->_index->end());
+        return true;
+      };
+      //
       void setInnerLow(const inner_key_type& ilow) {
         this->_ilow = ilow; this->_have_ilow = true;
       };
   //
   // Sifter definition
   template<typename Arrow_, 
-           typename ArrowSupportOrder_= SifterDef::ColorTargetOrder<Arrow_>, 
-           typename ArrowConeOrder_   = SifterDef::ColorSourceOrder<Arrow_>, 
+           typename ArrowSupportOrder_= SifterDef::TargetColorOrder<Arrow_>, 
+           typename ArrowConeOrder_   = SifterDef::SourceColorOrder<Arrow_>, 
            typename Predicate_ = int, typename PredicateOrder_ = std::less<Predicate_> >
   struct Sifter { // struct Sifter
     //
     public:
       // Predicate stored alongside the arrow data
       predicate_type _predicate;
+      // Basic interface
+      Rec(const arrow_type& a) : arrow_type(a) {};
+      Rec(const arrow_type& a, const predicate_type& p) : arrow_type(a), _predicate(p) {};
+      // Extended interface
       predicate_type predicate() const{return this->_predicate;};
       source_type    source() const {return this->arrow_type::source();};
       target_type    target() const {return this->arrow_type::target();};
-    };
+      // Modifier objects
+      struct predicateChanger {
+        predicateChanger(const predicate_type& newPredicate) : _newPredicate(newPredicate) {};
+        void operator()(Rec& r) { r._predicate = this->_newPredicate;}
+      private:
+        const predicate_type _newPredicate;
+      };
+    };// struct Rec
     //
     typedef Rec                              rec_type;
     //
       public SifterDef::StridedIndexSequence<Index_, OuterKeyExtractor_, InnerKeyExtractor_, ValueExtractor_, inner_strided_flag> {
       // ArrowSequence extends StridedIndexSequence with extra iterator methods.
     public:
-      typedef SifterDef::StridedIndexSequence<Index_, OuterKeyExtractor_, InnerKeyExtractor_, ValueExtractor_> super;
-      typedef Sifter                                                                                           container_type;
-      typedef typename super::index_type                                                                       index_type;
-      typedef typename super::outer_key_type                                                                   outer_key_type;
-      typedef typename super::inner_key_type                                                                   inner_key_type;
+      typedef SifterDef::StridedIndexSequence<Index_, OuterKeyExtractor_, InnerKeyExtractor_, ValueExtractor_, inner_strided_flag> super;
+      typedef Sifter                                                                                                               container_type;
+      typedef typename super::index_type                                                                                           index_type;
+      typedef typename super::outer_key_type                                                                                       outer_key_type;
+      typedef typename super::inner_key_type                                                                                       inner_key_type;
       
       // Need to extend the inherited iterators to be able to extract arrow color
       class iterator : public super::iterator {
       public:
+        iterator() : super::iterator() {};
         iterator(const typename super::iterator& super_iter) : super::iterator(super_iter) {};
         virtual const source_type& source() const {return this->_itor->_source;};
         virtual const color_type&  color()  const {return this->_itor->_color;};
         virtual const arrow_type&  arrow()  const {return *(this->_itor);};
       };
     protected:
-      container_type _container;
+      container_type *_container;
     public:
       //
       // Basic ArrowSequence interface
       //
+      ArrowSequence() : super(), _container(NULL) {};
       ArrowSequence(const ArrowSequence& seq) : super(seq), _container(seq._container) {};
-      ArrowSequence(const container_type& container, index_type& index) : super(index), _container(container) {};
+      ArrowSequence(container_type *container, index_type *index) : super(index), _container(container) {};
       virtual ~ArrowSequence() {};
       void copy(const ArrowSequence& seq, ArrowSequence& cseq) {
         super::copy(seq,cseq);
       ArrowSequence& operator=(const ArrowSequence& seq) {
         copy(seq,*this); return *this;
       };
-      void reset(const container_type& container, index_type& index) {
+      void reset(container_type *container, index_type *index) {
         this->super::reset(index);
         this->_container = container;
       };
       //
       // Extended ArrowSequence interface
       //
-      
       virtual iterator begin() {
-        return super::begin();
+        return this->super::begin();
       };
       //
       virtual iterator end() {
-        return super::end();
+        return this->super::end();
       };
       //
       template<typename ostream_type>
         os << " ]" << std::endl;
       };
       void addArrow(const arrow_type& a) {
-        this->_container.addArrow(a);
+        this->_container->addArrow(a);
       };
       //
-      virtual bool contains(const outer_key_type& ok, const inner_key_type& ik) {
-        return (this->_index.find(ALE::pair<outer_key_type,inner_key_type>(ok,ik)) != this->_index.end());
-      };
     };// class ArrowSequence    
       
     //
     //
     typedef ArrowSequence<typename ::boost::multi_index::index<rec_set_type, UpwardTag>::type,
                           ::boost::multi_index::const_mem_fun<rec_type, predicate_type, &rec_type::predicate>,
-                          ::boost::multi_index::identity<rec_type>,
+                          ::boost::multi_index::const_mem_fun<rec_type, target_type, &rec_type::target>,
                           ::boost::multi_index::const_mem_fun<rec_type, target_type, &rec_type::target>, 
                           true>                                                       
     BaseSequence;
 
     typedef ArrowSequence<typename ::boost::multi_index::index<rec_set_type, UpwardTag>::type,
                           ::boost::multi_index::const_mem_fun<rec_type, predicate_type, &rec_type::predicate>,
-                          ::boost::multi_index::identity<rec_type>,
+                          ::boost::multi_index::const_mem_fun<rec_type, target_type, &rec_type::target>,
                           ::boost::multi_index::const_mem_fun<rec_type, source_type, &rec_type::source> >     
     ConeSequence;
     //
     // Extended interface
     //
     void addArrow(const arrow_type& a) {
-      this->_rec_set.insert(a);
+      this->_rec_set.insert(rec_type(a));
     };
     void cone(const target_type& t, ConeSequence& seq) {
-      seq.reset(*this, ::boost::multi_index::get<UpwardTag>(this->_rec_set));
+      seq.reset(this, &::boost::multi_index::get<UpwardTag>(this->_rec_set));
       seq.setInnerLimits(t,t);
     };
     ConeSequence& cone(const target_type& t) {
       return cseq;
     };
     void base(BaseSequence& seq) {
-      seq.reset(*this, ::boost::multi_index::get<UpwardTag>(this->_rec_set));
+      seq.reset(this, &::boost::multi_index::get<UpwardTag>(this->_rec_set));
     };
     BaseSequence& base() {
       static BaseSequence bseq;

src/docs/tex/manual/manual.tex

 {\rm This manual is intended for use with PETSc 2.3.1}}
 {95/11 - Revision 2.3.1}{February, 3, 2006}
 
-\newpage
-
-\hbox{ }
+\cleardoublepage
 
 \vspace{1in}
 \date{\today}
 \addcontentsline{toc}{chapter}{Abstract}
 \input{abstract.tex}
 
+\cleardoublepage
+
 \input{gettinginfo.tex}
 
 \medskip \medskip
 
-%
+\cleardoublepage
+
 % Acknowledgements for users manual
-\newpage \hbox{ } \newpage
 \input{acknowl.tex}
 
 % Blank page makes double sided printout look bettter.
-%\newpage \hbox{ } \newpage
+
+\cleardoublepage
 
 \tableofcontents
 
 % --------------------------------------------------------------------
 %                            PART 1
 % --------------------------------------------------------------------
+\cleardoublepage
 \part{Introduction to PETSc}
 \label{part_intro}
+\cleardoublepage
 \chapter{Getting Started}
 \input{part1tmp.tex}
 
 % --------------------------------------------------------------------
 %                            PART 2
 % --------------------------------------------------------------------
+\cleardoublepage
 \part{Programming with PETSc}
 \label{part_usage}
 \input{part2tmp.tex}
 %------------------------------------------------------------------
 
 
+\cleardoublepage
 \bibliographystyle{plain}
 \addtocounter{chapter}{1}
 \addcontentsline{toc}{chapter}{Bibliography}

src/docs/tex/manual/part2.tex

 %    \tt mode from running of the page)
 %
 
+\cleardoublepage
 \chapter{Vectors and Distributing Parallel Data} 
 \label{chapter_vectors}
 \sindex{vectors}
 
 %-------------------------------------------------------------
 %-------------------------------------------------------------
+\cleardoublepage
 \chapter{Matrices}
 \label{chapter_matrices}
 \sindex{matrices}
 
 
 % ------------------------------------------------------------------
+\cleardoublepage
 \chapter{KSP: Linear Equations Solvers} \sindex{linear system solvers}
 \label{ch_ksp}
 
 
 
 % ---------------------------------------------------------------
+\cleardoublepage
 \chapter{SNES: Nonlinear Solvers}
 \sindex{nonlinear equation solvers}
 \label{chapter_snes}
 supported for the AIJ matrix format.
 
 % ---------------------------------------------------------------
+\cleardoublepage
 \chapter{TS: Scalable ODE Solvers}
 \label{chapter_ts}
 \sindex{ODE solvers}
 http://www.parallab.uib.no/projects/molecul/matrix/
 
 
+\cleardoublepage
 \chapter{High Level Support for Multigrid with DMMG}
 \label{chapter_dmmg}
 
 option \trl{-dmmg_grid_sequence} and \trl{-snes_view} for most runs.
 
 
+\cleardoublepage
 \chapter{Using ADIC and ADIFOR with PETSc}
 
 Automatic differentiation is an incredible technique to generate code
 \end{tabbing}
 
 %------------------------------------------------------------------
+\cleardoublepage
 \chapter{Using Matlab with PETSc}
 \label{ch_matlab}\sindex{Matlab}
 
 with \trl{PETSC_MATLAB_ENGINE_(MPI_Comm)}.
 
 
+\cleardoublepage
 \chapter{Using ESI with PETSc}
 \label{ch_esi} \findex{ESI} \findex{Equation Solver Interface}
 The \href{http://z.ca.sandia.gov/esi}{Equation Solver Interface} 
 
 
 %------------------------------------------------------------------
+\cleardoublepage
 \chapter{PETSc for Fortran Users}
 \label{ch_fortran}
 
 % --------------------------------------------------------------------
 %                            PART 3
 % --------------------------------------------------------------------
+\cleardoublepage
 \part{Additional Information}
 \label{part_usefulstuff}
 
 %---------------------------------------------------------------------
+\cleardoublepage
 \chapter{Profiling} 
 \label{ch_profiling} \sindex{profiling}
 
 \findex{PreLoadEnd()} \findex{PreLoadStage()} \findex{-preload}
 
 %---------------------------------------------------------------------
+\cleardoublepage
 \chapter{Hints for Performance Tuning} 
 \label{ch_performance} \sindex{performance tuning}
 
 \end{itemize}
 
 %---------------------------------------------------------------------
+\cleardoublepage
 \chapter{Other PETSc Features}
 
 \section{PETSc on a process subset}
 
 
 %----------------------------------------------------------------------
+\cleardoublepage
 \chapter{Makefiles}
 \label{ch_makefiles}
 
 
 %------------------------------------------------------------------
 
+\cleardoublepage
 \chapter{Unimportant and Advanced Features of Matrices and Solvers}
 \label{ch_advanced}
 
   PCDestroy(PC pc);
 \end{tabbing}
 
+\cleardoublepage
 \addtocounter{chapter}{1}
 \addcontentsline{toc}{chapter}{Index}
 \begin{theindex}