From ff3b59df665d781db8826999ac93e2a9e8af9571 Mon Sep 17 00:00:00 2001
From: "denis.davydov" <denis.davydov@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Wed, 21 May 2014 19:58:56 +0000
Subject: [PATCH] updated the test for matrix_scalar_product and
 matrix_norm_square

git-svn-id: https://svn.dealii.org/branches/branch_petscscalar_complex@32962 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../parallel_sparse_matrix_01.cc              | 97 ++++++++-----------
 .../parallel_sparse_matrix_01.mpirun=3.output | 12 +++
 .../parallel_sparse_matrix_01.output          |  6 --
 3 files changed, 51 insertions(+), 64 deletions(-)
 create mode 100644 tests/petsc_complex/parallel_sparse_matrix_01.mpirun=3.output
 delete mode 100644 tests/petsc_complex/parallel_sparse_matrix_01.output

diff --git a/tests/petsc_complex/parallel_sparse_matrix_01.cc b/tests/petsc_complex/parallel_sparse_matrix_01.cc
index 16f73b2cae..0483a943c8 100644
--- a/tests/petsc_complex/parallel_sparse_matrix_01.cc
+++ b/tests/petsc_complex/parallel_sparse_matrix_01.cc
@@ -49,43 +49,22 @@
 
 using namespace dealii;
 
-static const unsigned int dim = 2;
-
-void test ()
+template<int dim>
+void test (const unsigned int poly_degree = 1)
 {
 
   MPI_Comm mpi_communicator (MPI_COMM_WORLD);
    const unsigned int n_mpi_processes (Utilities::MPI::n_mpi_processes(mpi_communicator));
    const unsigned int this_mpi_process (Utilities::MPI::this_mpi_process(mpi_communicator));
-   
-   ConditionalOStream pcout (std::cout,
-           (Utilities::MPI::this_mpi_process(mpi_communicator)
-            == 0));
   
   parallel::distributed::Triangulation<dim> tria(mpi_communicator,
                    typename Triangulation<dim>::MeshSmoothing
                    (Triangulation<dim>::smoothing_on_refinement |
                     Triangulation<dim>::smoothing_on_coarsening));
 
-  if (false)
-   {
-    Triangulation<dim>   triangulation1;
-    Triangulation<dim>   triangulation2;
-    Triangulation<dim>   triangulation12;
-    GridGenerator::hyper_cube (triangulation1, -1,0); //create a square [-1,0]^d domain
-    GridGenerator::hyper_cube (triangulation2, -1,0); //create a square [-1,0]^d domain
-    Point<dim> shift_vector;
-    shift_vector[0] = 1.0;
-    GridTools::shift(shift_vector,triangulation2);
-    GridGenerator::merge_triangulations (triangulation1, triangulation2, tria);
-   }  
-  else
-  {
-     GridGenerator::hyper_cube (tria, -1,0);
-	 //tria.refine_global (1);
-  } 
+  GridGenerator::hyper_cube (tria, -1,0);
+  tria.refine_global (3);
   
-  const unsigned int poly_degree = 1;
   FE_Q<dim> fe(poly_degree);
       
   DoFHandler<dim> dof_handler(tria);
@@ -95,7 +74,7 @@ void test ()
   IndexSet locally_relevant_dofs;
   DoFTools::extract_locally_relevant_dofs (dof_handler,locally_relevant_dofs);
   
-  PETScWrappers::MPI::Vector vector, conjugate;
+  PETScWrappers::MPI::Vector vector;
   PETScWrappers::MPI::SparseMatrix mass_matrix;
   
 
@@ -116,8 +95,6 @@ void test ()
   //assemble mass matrix:
   mass_matrix = 0.0;
 	  {
-		  //a separate quadrature formula
-		  //enough for mass and kinetic matrices assembly
 		  QGauss<dim> quadrature_formula(poly_degree+1);
 		  FEValues<dim> fe_values (fe, quadrature_formula,
 				                   update_values | update_gradients | update_JxW_values);
@@ -158,44 +135,48 @@ void test ()
 			 mass_matrix.compress (VectorOperation::add);
 	  }			   
   
-  vector = std::complex<double>(1.0,-1.0);
-  
-  //deallog<<"original vector: "<<std::endl;
-  //vector.print(std::cout);
+  for (unsigned int i = 0; i < locally_owned_dofs.n_elements(); i++) {
+       double re = 0,
+              im = 0;
+      if ( i % 2 ) {
+         re = 1.0*i;
+         im = 1.0*(this_mpi_process+1);
+      } else if (i % 3) {
+      	 re = 0.0;
+      	 im = -1.0*(this_mpi_process+1);
+      } else {
+         re = 3.0*i;
+         im = 0.0;
+      }
+      const PetscScalar val = re+im*PETSC_i;
+      vector   ( locally_owned_dofs.nth_index_in_set(i)) = val;
+  }
+  vector.compress(VectorOperation::insert);	
+  constraints.distribute(vector);
   
   deallog<<"symmetric: "<<mass_matrix.is_symmetric()<<std::endl;
   deallog<<"hermitian: "<<mass_matrix.is_hermitian()<<std::endl;
   
-  //mass_matrix.print(std::cout);
-  
   PETScWrappers::MPI::Vector tmp(vector);
   mass_matrix.vmult (tmp, vector);
   
-  //pcout<<"vmult: "<<std::endl;
-  //tmp.print(std::cout);
-  //correct output ^^^ is:
-  // 0.25000 - 0.25000i
-  // 0.25000 - 0.25000i
-  // 0.25000 - 0.25000i
-  // 0.25000 - 0.25000i
-  
-  conjugate = vector;
-  conjugate.conjugate();
-  //pcout<<"conjugate vector: "<<std::endl;
-  //conjugate.print(std::cout);
-  
-  
-  deallog<<"conj*vmult:"<< conjugate*tmp<<std::endl;
-  
-  deallog<<"vector*vmult:"<< vector*tmp<<std::endl;
+  const std::complex<double> norm1 = vector*tmp;
+  deallog<<"(conj(vector),M vector): "<<std::endl;
+  deallog<<"real part:      "<<PetscRealPart(norm1)<<std::endl;
+  deallog<<"imaginary part: "<<PetscImaginaryPart(norm1)<<std::endl;
   
-  deallog<<"matrix_scalar_product: ";
-  const std::complex<double> norm =
+  const std::complex<double> norm2 =
 	  mass_matrix.matrix_scalar_product(vector, vector);
-  deallog<<norm<<std::endl;
-  //correct output is:
-  // (2,0)
-}//*/
+  deallog<<"matrix_scalar_product(vec,vec): "<<std::endl;
+  deallog<<"real part:      "<<PetscRealPart(norm2)<<std::endl;
+  deallog<<"imaginary part: "<<PetscImaginaryPart(norm2)<<std::endl;
+    
+  const std::complex<double> norm3 =
+	  mass_matrix.matrix_norm_square(vector);
+  deallog<<"matrix_norm_square(vec): "<<std::endl;
+  deallog<<"real part:      "<<PetscRealPart(norm3)<<std::endl;
+  deallog<<"imaginary part: "<<PetscImaginaryPart(norm3)<<std::endl;
+}
 
 
 int main (int argc, char *argv[])
@@ -207,6 +188,6 @@ int main (int argc, char *argv[])
   
    Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
    
-   test ();
+   test<2>();
    
 }
\ No newline at end of file
diff --git a/tests/petsc_complex/parallel_sparse_matrix_01.mpirun=3.output b/tests/petsc_complex/parallel_sparse_matrix_01.mpirun=3.output
new file mode 100644
index 0000000000..a75f75beb6
--- /dev/null
+++ b/tests/petsc_complex/parallel_sparse_matrix_01.mpirun=3.output
@@ -0,0 +1,12 @@
+
+DEAL::symmetric: 1
+DEAL::hermitian: 1
+DEAL::(conj(vector),M vector): 
+DEAL::real part:      337.563
+DEAL::imaginary part: 0
+DEAL::matrix_scalar_product(vec,vec): 
+DEAL::real part:      337.563
+DEAL::imaginary part: 0
+DEAL::matrix_norm_square(vec): 
+DEAL::real part:      337.563
+DEAL::imaginary part: 0
diff --git a/tests/petsc_complex/parallel_sparse_matrix_01.output b/tests/petsc_complex/parallel_sparse_matrix_01.output
deleted file mode 100644
index 52906fdae5..0000000000
--- a/tests/petsc_complex/parallel_sparse_matrix_01.output
+++ /dev/null
@@ -1,6 +0,0 @@
-
-DEAL::symmetric: 1
-DEAL::hermitian: 1
-DEAL::conj*vmult:(2.00000,0.00000)
-DEAL::vector*vmult:(0.00000,-2.00000)
-DEAL::matrix_scalar_product: (2.00000,0.00000)
-- 
2.39.5