Changes for PETSc 2.1.0

git-svn-id: http://svn.cactuscode.org/arrangements/CactusElliptic/EllPETSc/trunk@61 1d96b42b-98df-4a6a-9d84-1b24288d4588

7 files changed, 274 insertions, 181 deletions

diff --git a/src/Startup.c b/src/Startup.c
index 84c8b9c..24306df 100644
--- a/src/Startup.c
+++ b/src/Startup.c
@@ -9,38 +9,25 @@ static const char *rcsid = "$Header$";
 
 CCTK_FILEVERSION(CactusElliptic_EllPETSc_Startup_c)
 
+void petsc_confmetric(cGH *GH, int *MetricPsiI, int *FieldI, int *MI,
+		      int *NI, int *AbsTol, int *RelTol);
+
+void petsc_metric(cGH *GH, int *MetricI, int *FieldI, int *MI,
+		  int *NI, int *AbsTol, int *RelTol);
+
+void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, 
+		int *AbsTol, int *RelTol);
+
 /* Registration of the petsc solvers with the Elliptic solver registry.
    This routine registers petsc_confmetric under the name "petsc" for 
    the class of elliptic equations "LinConfMetric" */
    
 int EllPETSc_Register(cGH *GH) 
 {
-  void petsc_confmetric(cGH *GH, int *MetricPsiI, int *FieldI, int *MI,
-			int *NI, int *AbsTol, int *RelTol);
-  
-  void petsc_metric(cGH *GH, int *MetricI, int *FieldI, int *MI,
-		    int *NI, int *AbsTol, int *RelTol);
-
-  void petsc_flat(cGH *GH, int *FieldIndex, int *MIndex, int *NIndex, 
-			 int *AbsTol, int *RelTol);
   DECLARE_CCTK_PARAMETERS 
   
-  if (CCTK_Equals(elliptic_verbose,"yes"))
-  {
-      printf("PETSc: Registering petsc for Ell_LinConfMetric...\n");
-  }
   Ell_RegisterSolver(petsc_confmetric,"petsc","Ell_LinConfMetric");
-
-  if (CCTK_Equals(elliptic_verbose,"yes"))
-  {
-    printf("PETSc: Registering petsc for Ell_LinMetric...\n");
-  }
   Ell_RegisterSolver(petsc_metric,"petsc","Ell_LinMetric");
-
-  if (CCTK_Equals(elliptic_verbose,"yes"))
-  {
-    printf("PETSc: Registering petsc for Ell_LinFlat...\n");
-  }
   Ell_RegisterSolver(petsc_flat,"petsc","Ell_LinFlat");
 
   return 0;
diff --git a/src/ellpetsc.h b/src/ellpetsc.h
new file mode 100644
index 0000000..a1dceb2
--- /dev/null
+++ b/src/ellpetsc.h
@@ -0,0 +1,21 @@
+/*@@
+  @header ellpetsc.h
+  @author Gabrielle Allen
+  @date 17th September 2001
+  @desc
+  Include file for thorn EllPETSc
+  @enddesc
+  @version $Header$
+@@*/
+
+#ifndef _ELLPETSC_H_
+#define _ELLPETSC_H_ 1
+
+#define XDM 0
+#define XDP 1
+#define YDM 2
+#define YDP 3
+#define ZDM 4
+#define ZDP 5
+
+#endif /* _ELLPETSC_H_ */
diff --git a/src/make.configuration.defn b/src/make.configuration.defn
index 4fc63b8..7ea39c3 100644
--- a/src/make.configuration.defn
+++ b/src/make.configuration.defn
@@ -51,7 +51,7 @@ endif
 
 ### PETSC on linux
 ifeq ($(PETSC_ARCH),linux)
-     PLATFORM_LIBS = lapack blas g2c mpich
+     PLATFORM_LIBS = flapack fblas g2c mpich
 endif
 
 ### Otherwise
diff --git a/src/petsc_confmetric.c b/src/petsc_confmetric.c
index 96f1dbb..d85307a 100644
--- a/src/petsc_confmetric.c
+++ b/src/petsc_confmetric.c
@@ -15,8 +15,8 @@
 
 #include "assert.h"
 #include "mpi.h"
-#include "sles.h"
-
+#include "petscsles.h"
+#include "ellpetsc.h"
 #include "CactusPUGH/PUGH/src/include/pugh.h"
 
 static const char *rcsid = "$Header$";
@@ -354,16 +354,16 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
 			 19,PETSC_NULL,           /* Diagonal  */
 			 19,PETSC_NULL,           /* Off diagonal */
 			 &A[0]);                  /* The output  */
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
   ierr = VecCreateMPI(pughGH->PUGH_COMM_WORLD,(endpoint-startpoint),rank,&soln);
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
   ierr = VecCreateMPI(pughGH->PUGH_COMM_WORLD,(endpoint-startpoint),rank,&b);
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
   /*$}$*/ 
   /*$else {$*/
   /*  ierr = MatRestoreValuesMemory(A);  */
   /*$CCTK_WARN(0,"no reuse");
-    CHKERRA(ierr);
+    CHKERRQ(ierr);
   }$*/
 
 
@@ -702,7 +702,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
 		  else {
                     J = wsp[DATINDEX(pughGH,i+l,j+m,k+n)];
                     ierr = MatSetValues(A[0],1,&I,1,&J,&(a(l,m,n)),INSERT_VALUES); 
-		    CHKERRA(ierr);
+		    CHKERRQ(ierr);
                   }
                 }
               }
@@ -711,7 +711,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
 	  /* printf("RHS(%d %d %d)+N: %f",i,j,k,rhsval);*/
 
           ierr   = VecSetValues(b,1,&I,&rhsval,INSERT_VALUES);
-	  CHKERRA(ierr);
+	  CHKERRQ(ierr);
           
         }
       }
@@ -719,18 +719,18 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
   }
   
   if (verbose) CCTK_INFO ("Assembling the vectors");
-  ierr = MatAssemblyBegin(A[0],MAT_FINAL_ASSEMBLY); CHKERRA(ierr);
-  ierr = VecAssemblyBegin(soln);                    CHKERRA(ierr);
-  ierr = VecAssemblyBegin(b);                       CHKERRA(ierr)
-  ierr = MatAssemblyEnd(A[0],MAT_FINAL_ASSEMBLY);   CHKERRA(ierr);
-  ierr = VecAssemblyEnd(soln);                      CHKERRA(ierr);
-  ierr = VecAssemblyEnd(b);                         CHKERRA(ierr);
-  ierr = MatSetOption(A[0],MAT_NO_NEW_NONZERO_LOCATIONS); CHKERRA(ierr);
+  ierr = MatAssemblyBegin(A[0],MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
+  ierr = VecAssemblyBegin(soln);                    CHKERRQ(ierr);
+  ierr = VecAssemblyBegin(b);                       CHKERRQ(ierr)
+  ierr = MatAssemblyEnd(A[0],MAT_FINAL_ASSEMBLY);   CHKERRQ(ierr);
+  ierr = VecAssemblyEnd(soln);                      CHKERRQ(ierr);
+  ierr = VecAssemblyEnd(b);                         CHKERRQ(ierr);
+  ierr = MatSetOption(A[0],MAT_NO_NEW_NONZERO_LOCATIONS); CHKERRQ(ierr);
 
   OptionsSetValue("-ksp_monitor","");
   
   CCTK_INFO("CREATING SLES");
-  ierr = SLESCreate(pughGH->PUGH_COMM_WORLD,&sles);  CHKERRA(ierr);
+  ierr = SLESCreate(pughGH->PUGH_COMM_WORLD,&sles);  CHKERRQ(ierr);
   
   
   /* 
@@ -743,15 +743,15 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
      means, we'll leave it for now.  */
 
   CCTK_INFO("CREATING SLES OPERATOR");
-  ierr = SLESSetOperators(sles,A[0],A[0],DIFFERENT_NONZERO_PATTERN); CHKERRA(ierr);
+  ierr = SLESSetOperators(sles,A[0],A[0],DIFFERENT_NONZERO_PATTERN); CHKERRQ(ierr);
 
   /* Set linear solver defaults for this problem. Later this
      should be replaced/modified with appropriate parsing from the 
      parameter parser. For now it is not. These defaults are
      reasonable, I hope.  */
   CCTK_INFO("SLESGet KSP/PC");
-  ierr = SLESGetKSP(sles,&ksp); CHKERRA(ierr);
-  ierr = SLESGetPC(sles,&pc);   CHKERRA(ierr);
+  ierr = SLESGetKSP(sles,&ksp); CHKERRQ(ierr);
+  ierr = SLESGetPC(sles,&pc);   CHKERRQ(ierr);
   
  /* Get the PC Type */
   if (CCTK_Equals(petsc_PC_type,"PCJACOBI")) 
@@ -772,7 +772,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
     CCTK_WARN(1,"Don't understand petsc_PC_type. Using PCNONE\n");
     ierr = PCSetType(pc,PCNONE);
   }
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
 
   
   /* Now the same thing for the KSP Type */
@@ -800,7 +800,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
     CCTK_WARN (1,"I don't understand petsc_KSP_type. Using BiCGSTAB\n");
     ierr = KSPSetType(ksp,KSPBCGS);
   }
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
 
 
 
@@ -822,7 +822,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
   
   ierr = KSPSetTolerances(ksp,rtol,atol,PETSC_DEFAULT,
                           PETSC_DEFAULT);
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
 
 
     /* We are warned in the manual that 
@@ -849,7 +849,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
    */
   /*$ierr = KSPGMRESSetOrthogonalization(ksp, 
                 KSPGMRESUnmodifiedGramSchmidtOrthogonalization); 
-  CHKERRA(ierr);$*/
+  CHKERRQ(ierr);$*/
 
   /* We also learn that  ...
 
@@ -861,7 +861,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
 
    */
   CCTK_INFO("KSPSetInitialGuess\n");
-  ierr = KSPSetInitialGuessNonzero(ksp); CHKERRA(ierr);
+  ierr = KSPSetInitialGuessNonzero(ksp); CHKERRQ(ierr);
   
   /* 
     Set runtime options. Since we don't use PETSC runtime options
@@ -869,7 +869,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
     we parse things out. But that may not be such a good idea.
     So lets put it here.
   */
-  ierr = SLESSetFromOptions(sles);   CHKERRA(ierr);
+  ierr = SLESSetFromOptions(sles);   CHKERRQ(ierr);
     
 
   /* OK so finally we are able to ... */
@@ -878,7 +878,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
                       Solve the linear system
      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
   CCTK_INFO("SLES solve...");
-  ierr = SLESSolve(sles,b,soln,&its); CHKERRA(ierr); 
+  ierr = SLESSolve(sles,b,soln,&its); CHKERRQ(ierr); 
   
   /* Here we can form a "res = Ax  - b" and find its norm to get
      an idea of how well we've converged. Put this in later
@@ -925,14 +925,14 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
   PUGH_SyncGroup(GH,"ellpetsc::petscworkspace");
   
   /* And finally free up the matrix memory   FIXME 
-  ierr = MatReleaseValuesMemory(A); CHKERRA(ierr); */
+  ierr = MatReleaseValuesMemory(A); CHKERRQ(ierr); */
 
   /* This code is not used anymore */
   if (verbose) CCTK_INFO("Destroying Matrices");
-  ierr = SLESDestroy(sles); CHKERRA(ierr);
-  ierr = VecDestroy(soln);  CHKERRA(ierr);
-  ierr = VecDestroy(b);     CHKERRA(ierr);  
-  ierr = MatDestroy(A[0]);  CHKERRA(ierr);
+  ierr = SLESDestroy(sles); CHKERRQ(ierr);
+  ierr = VecDestroy(soln);  CHKERRQ(ierr);
+  ierr = VecDestroy(b);     CHKERRQ(ierr);  
+  ierr = MatDestroy(A[0]);  CHKERRQ(ierr);
   free(A);
     
   /*$PetscFinalize();$*/
diff --git a/src/petsc_confmetric_solver.c b/src/petsc_confmetric_solver.c
index 0219ed3..59441d3 100644
--- a/src/petsc_confmetric_solver.c
+++ b/src/petsc_confmetric_solver.c
@@ -15,8 +15,8 @@
 
 #include "assert.h"
 #include "mpi.h"
-#include "sles.h"
-
+#include "petscsles.h"
+#include "ellpetsc.h"
 #include "CactusPUGH/PUGH/src/include/pugh.h"
 
 static const char *rcsid = "$Header$";
@@ -200,7 +200,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
 
     /* Set the PETSc communicator to set of PUGH and
        initialzie PETSc */
-    ierr = PetscSetCommWorld(pughGH->PUGH_COMM_WORLD); CHKERRA(ierr);
+    ierr = PetscSetCommWorld(pughGH->PUGH_COMM_WORLD); CHKERRQ(ierr);
     PetscInitialize(&argc,&argv,NULL,NULL); 
 
     CCTK_INFO("PETSc initialized");
@@ -371,11 +371,11 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
 			 19,PETSC_NULL,           /* Diagonal  */
 			 19,PETSC_NULL,           /* Off diagonal */
 			 &A[0]);                  /* The output  */
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
   ierr = VecCreateMPI(pughGH->PUGH_COMM_WORLD,(endpoint-startpoint),rank,&soln);
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
   ierr = VecCreateMPI(pughGH->PUGH_COMM_WORLD,(endpoint-startpoint),rank,&b);
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
  
   /* Compare the PETSc layout to Cactus, this better be a match */
   ierr = VecGetOwnershipRange(soln,&pvstart,&pvend);
@@ -819,7 +819,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
 		  else {
                     J = wsp[DATINDEX(pEx,i+l,j+m,k+n)];
                     ierr = MatSetValues(A[0],1,&I,1,&J,&(a(l,m,n)),INSERT_VALUES); 
-		    CHKERRA(ierr);
+		    CHKERRQ(ierr);
                   }
                 }
               }
@@ -828,7 +828,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
 	    rhsval = -rhsval - Nlin[ijk] / ac;
 
           ierr   = VecSetValues(b,1,&I,&rhsval,INSERT_VALUES);
-	  CHKERRA(ierr);
+	  CHKERRQ(ierr);
           
         }
       }
@@ -836,13 +836,13 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
   }
   
   if (verbose) CCTK_INFO ("Assembling the vectors");
-  ierr = MatAssemblyBegin(A[0],MAT_FINAL_ASSEMBLY); CHKERRA(ierr);
-  ierr = VecAssemblyBegin(soln);                    CHKERRA(ierr);
-  ierr = VecAssemblyBegin(b);                       CHKERRA(ierr)
-  ierr = MatAssemblyEnd(A[0],MAT_FINAL_ASSEMBLY);   CHKERRA(ierr);
-  ierr = VecAssemblyEnd(soln);                      CHKERRA(ierr);
-  ierr = VecAssemblyEnd(b);                         CHKERRA(ierr);
-  ierr = MatSetOption(A[0],MAT_NO_NEW_NONZERO_LOCATIONS); CHKERRA(ierr);
+  ierr = MatAssemblyBegin(A[0],MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
+  ierr = VecAssemblyBegin(soln);                    CHKERRQ(ierr);
+  ierr = VecAssemblyBegin(b);                       CHKERRQ(ierr);
+  ierr = MatAssemblyEnd(A[0],MAT_FINAL_ASSEMBLY);   CHKERRQ(ierr);
+  ierr = VecAssemblyEnd(soln);                      CHKERRQ(ierr);
+  ierr = VecAssemblyEnd(b);                         CHKERRQ(ierr);
+  ierr = MatSetOption(A[0],MAT_NO_NEW_NONZERO_LOCATIONS); CHKERRQ(ierr);
 
   if (nabla_form == 3) {
     if (verbose)
@@ -855,10 +855,16 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
     free(uzz3);
   }
   if (trips==0) 
-    OptionsSetValue("-ksp_monitor","");
-  
-  if (verbose) CCTK_INFO("CREATING SLES");
-  ierr = SLESCreate(pughGH->PUGH_COMM_WORLD,&sles);  CHKERRA(ierr);
+  {
+    PetscOptionsSetValue("-ksp_monitor","");
+  }
+
+  if (verbose) 
+  {
+    CCTK_INFO("CREATING SLES");
+  }
+
+  ierr = SLESCreate(pughGH->PUGH_COMM_WORLD,&sles);  CHKERRQ(ierr);
   
   
   /* 
@@ -871,15 +877,15 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
      means, we'll leave it for now.  */
 
   if (verbose) CCTK_INFO("CREATING SLES OPERATOR");
-  ierr = SLESSetOperators(sles,A[0],A[0],DIFFERENT_NONZERO_PATTERN); CHKERRA(ierr);
+  ierr = SLESSetOperators(sles,A[0],A[0],DIFFERENT_NONZERO_PATTERN); CHKERRQ(ierr);
 
   /* Set linear solver defaults for this problem. Later this
      should be replaced/modified with appropriate parsing from the 
      parameter parser. For now it is not. These defaults are
      reasonable, I hope.  */
   if (verbose) CCTK_INFO("SLESGet KSP/PC");
-  ierr = SLESGetKSP(sles,&ksp); CHKERRA(ierr);
-  ierr = SLESGetPC(sles,&pc);   CHKERRA(ierr);
+  ierr = SLESGetKSP(sles,&ksp); CHKERRQ(ierr);
+  ierr = SLESGetPC(sles,&pc);   CHKERRQ(ierr);
   
  /* Get the PC Type */
   if (CCTK_Equals(petsc_PC_type,"PCJACOBI")) 
@@ -900,7 +906,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
     CCTK_WARN(1,"Don't understand petsc_PC_type. Using PCNONE\n");
     ierr = PCSetType(pc,PCNONE);
   }
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
 
   
   /* Now the same thing for the KSP Type */
@@ -928,7 +934,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
     CCTK_WARN (1,"I don't understand petsc_KSP_type. Using BiCGSTAB\n");
     ierr = KSPSetType(ksp,KSPBCGS);
   }
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
 
 
 
@@ -950,7 +956,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
   
   ierr = KSPSetTolerances(ksp,rtol,atol,PETSC_DEFAULT,
                           PETSC_DEFAULT);
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
 
 
     /* We are warned in the manual that 
@@ -977,7 +983,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
    */
   /*$ierr = KSPGMRESSetOrthogonalization(ksp, 
                 KSPGMRESUnmodifiedGramSchmidtOrthogonalization); 
-  CHKERRA(ierr);$*/
+  CHKERRQ(ierr);$*/
 
   /* We also learn that  ...
 
@@ -989,7 +995,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
 
    */
   if (verbose) CCTK_INFO("KSPSetInitialGuess\n");
-  ierr = KSPSetInitialGuessNonzero(ksp); CHKERRA(ierr);
+  ierr = KSPSetInitialGuessNonzero(ksp); CHKERRQ(ierr);
   
   /* 
     Set runtime options. Since we don't use PETSC runtime options
@@ -997,7 +1003,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
     we parse things out. But that may not be such a good idea.
     So lets put it here.
   */
-  ierr = SLESSetFromOptions(sles);   CHKERRA(ierr);
+  ierr = SLESSetFromOptions(sles);   CHKERRQ(ierr);
     
 
   /* OK so finally we are able to ... */
@@ -1006,7 +1012,7 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
                       Solve the linear system
      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
   if (verbose) CCTK_INFO("SLES solve...");
-  ierr = SLESSolve(sles,b,soln,&its); CHKERRA(ierr); 
+  ierr = SLESSolve(sles,b,soln,&its); CHKERRQ(ierr); 
   
   /* Here we can form a "res = Ax  - b" and find its norm to get
      an idea of how well we've converged. Put this in later
@@ -1047,14 +1053,14 @@ void petsc_confmetric_solver(cGH *GH, int *MetricPsiI, int MetricPsiISize,
   PUGH_SyncGroup(GH,"ellpetsc::petscworkspace");
   
   /* And finally free up the matrix memory   FIXME 
-  ierr = MatReleaseValuesMemory(A); CHKERRA(ierr); */
+  ierr = MatReleaseValuesMemory(A); CHKERRQ(ierr); */
 
   /* This code is not used anymore */
   if (verbose) CCTK_INFO("Destroying Matrices");
-  ierr = SLESDestroy(sles); CHKERRA(ierr);
-  ierr = VecDestroy(soln);  CHKERRA(ierr);
-  ierr = VecDestroy(b);     CHKERRA(ierr);  
-  ierr = MatDestroy(A[0]);  CHKERRA(ierr);
+  ierr = SLESDestroy(sles); CHKERRQ(ierr);
+  ierr = VecDestroy(soln);  CHKERRQ(ierr);
+  ierr = VecDestroy(b);     CHKERRQ(ierr);  
+  ierr = MatDestroy(A[0]);  CHKERRQ(ierr);
   free(A);
     
   /*$PetscFinalize();$*/
diff --git a/src/petsc_flat_solver.c b/src/petsc_flat_solver.c
index 74e0216..0a8ec9e 100644
--- a/src/petsc_flat_solver.c
+++ b/src/petsc_flat_solver.c
@@ -15,8 +15,8 @@
 #include "cctk_Parameters.h"
 
 #include "mpi.h"
-#include "sles.h"
-
+#include "petscsles.h"
+#include "ellpetsc.h"
 #include "CactusPUGH/PUGH/src/include/pugh.h"
 
 static const char *rcsid = "$Header$";
@@ -48,7 +48,8 @@ static SLES    sles;         /* linear solver context */
 
 
 void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex, 
-		int *AbsTol, int *RelTol) {
+		int *AbsTol, int *RelTol) 
+{
 
   DECLARE_CCTK_PARAMETERS    /* CCTK passed parameters */
 
@@ -105,10 +106,6 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   /* flags to signal if storage is ON or OFF (used/not-used) for M/N*/
   int Mstorage=0, Nstorage=0;
 
-
-
-
-
   octant       = CCTK_Equals(domain,"octant");
   verbose      = CCTK_Equals(petsc_verbose,"yes")||
                  CCTK_Equals(petsc_verbose,"debug");
@@ -122,14 +119,18 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
 
   /* Get the link to pugh Extension */
   pughGH = (pGH*)GH->extensions[CCTK_GHExtensionHandle("PUGH")];
-  if (!pughGH) CCTK_WARN(0,"ETERNAL ERROR: Cannot find PUGH Extension Handle\n");
+  if (!pughGH) 
+  {
+    CCTK_WARN(0,"ETERNAL ERROR: Cannot find PUGH Extension Handle");
+  }
 
   /* Get the extras extension for 3D grid functions */
   pEx = pughGH->GFExtras[2];
   pCon = pughGH->Connectivity[2];
 
   /* Things to do on first iteration */
-  if (trips==0) {
+  if (trips==0) 
+  {
     int argc;
     char **argv; 
 
@@ -139,7 +140,7 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
 #endif
 
     argc = CCTK_CommandLine(&argv);
-    ierr = PetscSetCommWorld(pughGH->PUGH_COMM_WORLD); CHKERRA(ierr);
+    ierr = PetscSetCommWorld(pughGH->PUGH_COMM_WORLD); CHKERRQ(ierr);
     PetscInitialize(&argc,&argv,NULL,NULL); 
 
     CCTK_INFO("PETSc initialized");
@@ -157,11 +158,13 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   /* if we have a negative index, this GF is not needed, 
      there for don't even look for it. when index positive,
      set flag Mstorage=1, dito for N */
-  if (MIndex>=0)  { 
+  if (MIndex>=0)  
+  { 
     Mlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,MIndex);
     Mstorage = 1;
   }
-  if (NIndex>=0) {
+  if (NIndex>=0) 
+  {
     Nlin = (CCTK_REAL*) CCTK_VarDataPtrI(GH,0,NIndex);
     Nstorage = 1;
   }
@@ -171,7 +174,10 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
 
   /* initialize the linear index lookup table with -1; 
      this will indicate a boudary later */
-  for (i=0;i<pEx->npoints;i++) wsp[i] = -1.0;
+  for (i=0;i<pEx->npoints;i++) 
+  {
+    wsp[i] = -1.0;
+  }
 
   /* Get myproc from the CCTK, get the gridspacing */
   myproc = PUGH_MyProc(GH);
@@ -199,7 +205,8 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
      (we use PUGH's "rn" variables (number of gp on proc[#] in 
      direction [i]) and subtract the ghostzones if nec.) */
 
-  for (i=0;i<myproc;i++) {
+  for (i=0;i<myproc;i++) 
+  {
     
     nxs=pEx->rnsize[i][0];
     nxs=((pCon->neighbours[i][XDM]<0) && !(octant) ? nxs-1 : 
@@ -241,7 +248,10 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
      do a sync that means each processor knows the indiced
      in the ghost zones (eg, on its neighbors) in a FD stencil of 1 */
   retcode = PUGH_SyncGroup(GH,"ellpetsc::petscworkspace");
-  if (retcode<0) CCTK_WARN(1,"Synchronization failed\n");
+  if (retcode<0) 
+  {
+    CCTK_WARN(1,"Synchronization failed");
+  }
   
   /* So woohoo. Now for each point in our ijk space, we have
      information about our row in the matrix (workspace->data[ijk])
@@ -253,7 +263,8 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   nzs = pEx->nsize[2]-2;
 
 #ifdef DEBUG
-  if (debug) {
+  if (debug) 
+  {
     printf("nxyzs %d %d %d \n",nxs,nys,nzs);
     printf("lnxyz: %d %d %d \n",pEx->lnsize[0],pEx->lnsize[1],pEx->lnsize[2]);
   }
@@ -291,23 +302,24 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
 			 19,PETSC_NULL,           /* Diagonal  */
 			 19,PETSC_NULL,           /* Off diagonal */
 			 &A[0]);                  /* The output  */
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
   ierr = VecCreateMPI(pughGH->PUGH_COMM_WORLD,(endpoint-startpoint),rank,&soln);
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
   ierr = VecCreateMPI(pughGH->PUGH_COMM_WORLD,(endpoint-startpoint),rank,&b);
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
  
 
   /* Compare the PETSc layout to Cactus, this better be a match */
   ierr = VecGetOwnershipRange(soln,&pvstart,&pvend);
   ierr = MatGetOwnershipRange(A[0],&pstart,&pend);
-  if (debug) { 
+  if (debug) 
+  { 
     printf("CAC M-Layout: %d %d \n",startpoint, endpoint);
     printf("PET M-Layout: %d %d \n",pstart, pend);
     printf("PET V-Layout: %d %d \n",pvstart,pvend);
   }
   if (pstart != startpoint && pend != endpoint) 
-    CCTK_WARN(1,"WARNING: PETSC and data layouts differ! (why??)\n");
+    CCTK_WARN(1,"WARNING: PETSC and data layouts differ! (why??)");
          
 
   /* Initialize the stencil array */
@@ -321,11 +333,15 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   a( 0, 0,-1) = 1.0;  /* ab */
 
 
-  for (k=kmin;k<kmax;k++) {
-    for (j=jmin;j<jmax;j++) {
-      for (i=imin;i<imax;i++) {
+  for (k=kmin;k<kmax;k++) 
+  {
+    for (j=jmin;j<jmax;j++) 
+    {
+      for (i=imin;i<imax;i++) 
+      {
 
-        if (wsp[DATINDEX(pEx,i,j,k)] >= 0) {
+        if (wsp[DATINDEX(pEx,i,j,k)] >= 0) 
+        {
 
           /* Set up indices */
           int ijk;              /* The lineax index for the array */
@@ -342,16 +358,27 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
           I = wsp[ijk];
        
           /* M phi */
-	  if (Mstorage) a(0,0,0) = a(0,0,0) + Mlin[ijk];
+	  if (Mstorage) 
+	  {
+	    a(0,0,0) = a(0,0,0) + Mlin[ijk];
+	  }
 
-          /* Great now set the values of the matrix. Here we force dirichlet. */
+          /* Now set the values of the matrix. Here we force dirichlet. */
 	  VecSetValues(soln,1,&I,&(ell_field[ijk]),INSERT_VALUES);
           
 #ifdef DEBUG
-	  if (ijk==1918) {
+	  if (ijk==1918) 
+	  {
 	    for (l=-1;l<=1;l++)
+	    {
 	      for (m=-1;m<=1;m++)
-		for (n=-1;n<=1;n++) printf(" (%d %d %d): %f \n",l,m,n,a(l,m,n));
+	      {
+		for (n=-1;n<=1;n++) 
+		{
+		  printf(" (%d %d %d): %f \n",l,m,n,a(l,m,n));
+		}
+	      }
+	    }
 	    printf("\n");
 	  }
 #endif
@@ -376,7 +403,8 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
                       }
       
           /* renormalize */
-          if (matnormalize) {
+          if (matnormalize) 
+	  {
             ac = a(0,0,0);
             for (J=0;J<27;J++) a[J] = a[J] / ac;
           } 
@@ -390,27 +418,37 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
              19-way banded matrix. */
 
           for (l=-1;l<=1;l++)
+	  {
             for (m=-1;m<=1;m++)
-              for (n=-1;n<=1;n++) {
-                if (l*m*n == 0) {       /* This is the 19 point ... if none are
-                                         * zero, then we have no stencil here.
-                                         */
-                  if (wsp[DATINDEX(pEx,i+l,j+m,k+n)] < 0) {
+	    {
+              for (n=-1;n<=1;n++) 
+	      {
+                if (l*m*n == 0) 
+		{       /* This is the 19 point ... if none are
+			 * zero, then we have no stencil here.
+			 */
+                  if (wsp[DATINDEX(pEx,i+l,j+m,k+n)] < 0) 
+		  {
                     /* This is a boundary. */
                     rhsval += a(l,m,n) * ell_field[DATINDEX(pEx,i+l,j+m,k+n)];
                   } 
-		  else {
+		  else 
+		  {
                     J = wsp[DATINDEX(pEx,i+l,j+m,k+n)];
                     ierr = MatSetValues(A[0],1,&I,1,&J,&(a(l,m,n)),INSERT_VALUES); 
-		    CHKERRA(ierr);
+		    CHKERRQ(ierr);
                   }
                 }
               }
-	  
-	  if (Nstorage) rhsval = -rhsval - Nlin[ijk] / ac;
+	    }
+	  }
+	  if (Nstorage) 
+	  {
+	    rhsval = -rhsval - Nlin[ijk] / ac;
+	  }
 
           ierr   = VecSetValues(b,1,&I,&rhsval,INSERT_VALUES);
-	  CHKERRA(ierr);
+	  CHKERRQ(ierr);
           
 	  /* reset the central value*/
 	  a( 0, 0, 0) = -6.0;
@@ -420,20 +458,29 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
     }   /* for j */
   }   /* for k */
   
-  if (verbose) CCTK_INFO ("Assembling the vectors");
-  ierr = MatAssemblyBegin(A[0],MAT_FINAL_ASSEMBLY); CHKERRA(ierr);
-  ierr = VecAssemblyBegin(soln);                    CHKERRA(ierr);
-  ierr = VecAssemblyBegin(b);                       CHKERRA(ierr)
-  ierr = MatAssemblyEnd(A[0],MAT_FINAL_ASSEMBLY);   CHKERRA(ierr);
-  ierr = VecAssemblyEnd(soln);                      CHKERRA(ierr);
-  ierr = VecAssemblyEnd(b);                         CHKERRA(ierr);
-  ierr = MatSetOption(A[0],MAT_NO_NEW_NONZERO_LOCATIONS); CHKERRA(ierr);
+  if (verbose) 
+  {
+    CCTK_INFO ("Assembling the vectors");
+  }
+  ierr = MatAssemblyBegin(A[0],MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
+  ierr = VecAssemblyBegin(soln);                    CHKERRQ(ierr);
+  ierr = VecAssemblyBegin(b);                       CHKERRQ(ierr);
+  ierr = MatAssemblyEnd(A[0],MAT_FINAL_ASSEMBLY);   CHKERRQ(ierr);
+  ierr = VecAssemblyEnd(soln);                      CHKERRQ(ierr);
+  ierr = VecAssemblyEnd(b);                         CHKERRQ(ierr);
+  ierr = MatSetOption(A[0],MAT_NO_NEW_NONZERO_LOCATIONS); CHKERRQ(ierr);
   
   if (trips==0) 
-    OptionsSetValue("-ksp_monitor","");
-  
-  if (verbose) CCTK_INFO("CREATING SLES");
-  ierr = SLESCreate(pughGH->PUGH_COMM_WORLD,&sles);  CHKERRA(ierr);
+  {
+    PetscOptionsSetValue("-ksp_monitor","");
+  }
+
+  if (verbose) 
+  {
+    CCTK_INFO("CREATING SLES");
+  }
+
+  ierr = SLESCreate(pughGH->PUGH_COMM_WORLD,&sles);  CHKERRQ(ierr);
   
   
   /* 
@@ -445,16 +492,22 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
      using A to precondition itself. Since I don't know what this
      means, we'll leave it for now.  */
 
-  if (verbose) CCTK_INFO("CREATING SLES OPERATOR");
-  ierr = SLESSetOperators(sles,A[0],A[0],DIFFERENT_NONZERO_PATTERN); CHKERRA(ierr);
+  if (verbose) 
+  {
+    CCTK_INFO("CREATING SLES OPERATOR");
+  }
+  ierr = SLESSetOperators(sles,A[0],A[0],DIFFERENT_NONZERO_PATTERN); CHKERRQ(ierr);
 
   /* Set linear solver defaults for this problem. Later this
      should be replaced/modified with appropriate parsing from the 
      parameter parser. For now it is not. These defaults are
      reasonable, I hope.  */
-  if (verbose) CCTK_INFO("SLESGet KSP/PC");
-  ierr = SLESGetKSP(sles,&ksp); CHKERRA(ierr);
-  ierr = SLESGetPC(sles,&pc);   CHKERRA(ierr);
+  if (verbose) 
+  {
+    CCTK_INFO("SLESGet KSP/PC");
+  }
+  ierr = SLESGetKSP(sles,&ksp); CHKERRQ(ierr);
+  ierr = SLESGetPC(sles,&pc);   CHKERRQ(ierr);
   
  /* Get the PC Type */
   if (CCTK_Equals(petsc_PC_type,"PCJACOBI")) 
@@ -471,11 +524,12 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
     ierr = PCSetType(pc,PCLU);
   else if (CCTK_Equals(petsc_PC_type,"PCNONE"))
     ierr = PCSetType(pc,PCNONE);
-  else {
-    CCTK_WARN(1,"Don't understand petsc_PC_type. Using PCNONE\n");
+  else 
+  {
+    CCTK_WARN(1,"Don't understand petsc_PC_type. Using PCNONE");
     ierr = PCSetType(pc,PCNONE);
   }
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
 
   
   /* Now the same thing for the KSP Type */
@@ -499,36 +553,47 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
     ierr = KSPSetType(ksp,KSPCR);
   else if (CCTK_Equals(petsc_KSP_type,"KSPLSQR"))
     ierr = KSPSetType(ksp,KSPLSQR);
-  else {
-    CCTK_WARN (1,"I don't understand petsc_KSP_type. Using BiCGSTAB\n");
+  else 
+  {
+    CCTK_WARN (1,"I don't understand petsc_KSP_type. Using BiCGSTAB");
     ierr = KSPSetType(ksp,KSPBCGS);
   }
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
 
 
 
  /* Set up tolerances */
   
-  if (PetscTolStyle == ELLCONV_ABSOLUTE) {
+  if (PetscTolStyle == ELLCONV_ABSOLUTE) 
+  {
     rtol = 1.0e-15; 
     atol = tolerance;
-  } else if (PetscTolStyle == ELLCONV_RELATIVE) {
+  } 
+  else if (PetscTolStyle == ELLCONV_RELATIVE) 
+  {
     rtol = tolerance; 
     atol = 1.0e-15;
-  } else if (PetscTolStyle == ELLCONV_EITHER) {
+  } 
+  else if (PetscTolStyle == ELLCONV_EITHER) 
+  {
     rtol = tolerance; 
     atol = tolerance;
-  } else {
+  } 
+  else 
+  {
     printf("PETSC Solver: PetscTolStyle set incorrectly [%d]\n",
             PetscTolStyle);
   }
   
   ierr = KSPSetTolerances(ksp, rtol, atol, PETSC_DEFAULT, PETSC_DEFAULT);
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
 
-  if (verbose) CCTK_INFO("KSPSetInitialGuess\n");
+  if (verbose) 
+  {
+    CCTK_INFO("KSPSetInitialGuess");
+  }
   ierr = KSPSetInitialGuessNonzero(ksp); 
-  CHKERRA(ierr);
+  CHKERRQ(ierr);
   
   /* 
     Set runtime options. Since we don't use PETSC runtime options
@@ -536,7 +601,7 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
     we parse things out. But that may not be such a good idea.
     So lets put it here.
   */
-  ierr = SLESSetFromOptions(sles);   CHKERRA(ierr);
+  ierr = SLESSetFromOptions(sles);   CHKERRQ(ierr);
     
 
   /* OK so finally we are able to ... */
@@ -544,8 +609,12 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
                       Solve the linear system
      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-  if (verbose) CCTK_INFO("SLES solve...");
-  ierr = SLESSolve(sles,b,soln,&its); CHKERRA(ierr); 
+  if (verbose) 
+  {
+    CCTK_INFO("SLES solve...");
+  }
+
+  ierr = SLESSolve(sles,b,soln,&its); CHKERRQ(ierr); 
   
   /* Here we can form a "res = Ax  - b" and find its norm to get
      an idea of how well we've converged. Put this in later
@@ -559,20 +628,22 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
      on (eg, the VecCreateMPI has made the solution local) */
 
   VecGetArray(soln,&values); 
-  for (k=kmin;k<kmax;k++) {
-    for (j=jmin;j<jmax;j++) {
-      for (i=imin;i<imax;i++) {
+  for (k=kmin;k<kmax;k++) 
+  {
+    for (j=jmin;j<jmax;j++) 
+    {
+      for (i=imin;i<imax;i++) 
+      {
         /* Set up indices */
         int ijk;                /* The data point for the array */
         int I;                  /* The col in the matrix */
 
 	/* these guys are easy */
         ijk = DATINDEX(pEx,i,j,k);
-        if (wsp[ijk] >= 0) {
-           
+        if (wsp[ijk] >= 0) 
+        {
           /* Now this one. "Fortran-order" the matrix. But remember
-             we have stripped off the ghost zones. Hence ig-1 and nxs...
-          */
+             we have stripped off the ghost zones. Hence ig-1 and nxs */
           I = wsp[ijk]-startpoint;
           ell_field[ijk] = values[I];
         }
@@ -585,19 +656,25 @@ void petsc_flat(cGH *GH, int FieldIndex, int MIndex, int NIndex,
   PUGH_SyncGroup(GH,"ellpetsc::petscworkspace");
   
   /* And finally free up the matrix memory */  
-  /*$ierr = MatReleaseValuesMemory(A); CHKERRA(ierr);$*/ 
+  /*$ierr = MatReleaseValuesMemory(A); CHKERRQ(ierr);$*/ 
 
   /* This code is not used anymore */
-  if (verbose) CCTK_INFO("Destroying Matrices");
-  ierr = SLESDestroy(sles); CHKERRA(ierr);
-  ierr = VecDestroy(soln);  CHKERRA(ierr);
-  ierr = VecDestroy(b);     CHKERRA(ierr);  
-  ierr = MatDestroy(A[0]);  CHKERRA(ierr);
+  if (verbose) 
+  {
+    CCTK_INFO("Destroying Matrices");
+  }
+  ierr = SLESDestroy(sles); CHKERRQ(ierr);
+  ierr = VecDestroy(soln);  CHKERRQ(ierr);
+  ierr = VecDestroy(b);     CHKERRQ(ierr);  
+  ierr = MatDestroy(A[0]);  CHKERRQ(ierr);
   free(A);
     
   /*$PetscFinalize();$*/
 
-  if (verbose) CCTK_INFO("LEAVING ELLPETSC");
+  if (verbose) 
+  {
+    CCTK_INFO("LEAVING ELLPETSC");
+  }
   
 }
 
diff --git a/src/petsc_wrapper.c b/src/petsc_wrapper.c
index 1f25b3e..74d3512 100644
--- a/src/petsc_wrapper.c
+++ b/src/petsc_wrapper.c
@@ -7,17 +7,17 @@ static const char *rcsid = "$Header$";
 
 CCTK_FILEVERSION(CactusElliptic_EllPETSc_petsc_wrapper_c)
 
-/* The wrapper functions for the core PETSc solver, that performs the  actual solve. 
-   One wrapper function for each class, because these functions are registered 
-   with a fixed set up arguments. Still, they all call petsc_confmetric_solver
-   as the core solver. */
+/* The wrapper functions for the core PETSc solver, that performs the 
+   actual solve. One wrapper function for each class, because these 
+   functions are registered with a fixed set up arguments. Still, they 
+   all call petsc_confmetric_solver as the core solver. */
 
-/* Wrapper function for the class of ellitpic equations that needs a metric */
+/* Wrapper function for the class of elliptic equations that needs a metric */
 void petsc_metric(cGH *GH, int *MetricI, int FieldIndex, 
 		  int MIndex, int NIndex, 
 		  CCTK_REAL *AbsTol, CCTK_REAL *RelTol) {
 
-  /* Arrays for the M/N size info */
+/* Arrays for the M/N size info */
   
   void petsc_confmetric_solver(cGH *GH, int *MetricI, int MetricISize, 
 			       int FieldIndex, int MIndex, int NIndex, 
@@ -25,7 +25,9 @@ void petsc_metric(cGH *GH, int *MetricI, int FieldIndex,
   
   
   if (GH->cctk_dim>3) 
+  {
    CCTK_WARN(0,"This elliptic solver implementation does not do dimension>3!");
+  }
 
   petsc_confmetric_solver(GH, MetricI, 6, FieldIndex, 
 			  MIndex, NIndex,