]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Make sure that we can compile lac/ with PETSc when the latter is using complex scalars.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 21 Sep 2009 03:47:43 +0000 (03:47 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 21 Sep 2009 03:47:43 +0000 (03:47 +0000)
git-svn-id: https://svn.dealii.org/trunk@19482 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/constraint_matrix.templates.h
deal.II/lac/include/lac/petsc_matrix_base.h
deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h
deal.II/lac/include/lac/petsc_solver.h
deal.II/lac/include/lac/petsc_vector_base.h
deal.II/lac/include/lac/vector.templates.h
deal.II/lac/source/petsc_matrix_base.cc
deal.II/lac/source/petsc_parallel_sparse_matrix.cc
deal.II/lac/source/petsc_solver.cc
deal.II/lac/source/petsc_vector_base.cc

index d022682537a0dd5f4d70a41c9cc1be04d32794ba..3fd91ec336be8b18a2352525866819a6d3775f70 100644 (file)
@@ -72,7 +72,7 @@ ConstraintMatrix::condense (const VectorType &uncondensed,
   Assert (condensed.size()+n_constraints() == uncondensed.size(),
          ExcDimensionMismatch(condensed.size()+n_constraints(),
                               uncondensed.size()));
-  
+
                                   // store for each line of the
                                   // vector its new line number after
                                   // compression. If the shift is -1,
@@ -85,7 +85,7 @@ ConstraintMatrix::condense (const VectorType &uncondensed,
   unsigned int                                shift           = 0;
   unsigned int n_rows = uncondensed.size();
 
-  if (next_constraint == lines.end()) 
+  if (next_constraint == lines.end())
                                     // if no constraint is to be handled
     for (unsigned int row=0; row!=n_rows; ++row)
       new_line.push_back (row);
@@ -127,7 +127,7 @@ ConstraintMatrix::condense (const VectorType &uncondensed,
     else
                                       // line must be distributed
       {
-       for (unsigned int q=0; q!=next_constraint->entries.size(); ++q) 
+       for (unsigned int q=0; q!=next_constraint->entries.size(); ++q)
          condensed(new_line[next_constraint->entries[q].first])
            +=
            uncondensed(row) * next_constraint->entries[q].second;
@@ -148,9 +148,10 @@ ConstraintMatrix::condense (VectorType &vec) const
   std::vector<ConstraintLine>::const_iterator constraint_line = lines.begin();
   for (; constraint_line!=lines.end(); ++constraint_line)
     {
-      for (unsigned int q=0; q!=constraint_line->entries.size(); ++q) 
+      for (unsigned int q=0; q!=constraint_line->entries.size(); ++q)
         vec(constraint_line->entries[q].first)
-          += (vec(constraint_line->line) *
+          += (static_cast<typename VectorType::value_type>
+             (vec(constraint_line->line)) *
               constraint_line->entries[q].second);
       vec(constraint_line->line) = 0.;
 
@@ -175,11 +176,11 @@ ConstraintMatrix::condense (const SparseMatrix<number> &uncondensed,
                                   // check whether we work on real vectors
                                   // or we just used a dummy when calling
                                   // the other function above.
-  const bool use_vectors = (uncondensed_vector.size() == 0 && 
+  const bool use_vectors = (uncondensed_vector.size() == 0 &&
                            condensed_vector.size() == 0) ? false : true;
 
   const SparsityPattern &uncondensed_struct = uncondensed.get_sparsity_pattern ();
-  
+
   Assert (sorted == true, ExcMatrixNotClosed());
   Assert (uncondensed_struct.is_compressed() == true, ExcMatrixNotClosed());
   Assert (condensed.get_sparsity_pattern().is_compressed() == true, ExcMatrixNotClosed());
@@ -210,7 +211,7 @@ ConstraintMatrix::condense (const SparseMatrix<number> &uncondensed,
   unsigned int                                shift           = 0;
   const unsigned int n_rows = uncondensed_struct.n_rows();
 
-  if (next_constraint == lines.end()) 
+  if (next_constraint == lines.end())
                                     // if no constraint is to be handled
     for (unsigned int row=0; row!=n_rows; ++row)
       new_line.push_back (row);
@@ -256,7 +257,7 @@ ConstraintMatrix::condense (const SparseMatrix<number> &uncondensed,
          if (new_line[uncondensed_struct.get_column_numbers()[j]] != -1)
            condensed.add (new_line[row], new_line[uncondensed_struct.get_column_numbers()[j]],
                           uncondensed.global_entry(j));
-         else 
+         else
            {
                                             // let c point to the
                                             // constraint of this column
@@ -277,12 +278,12 @@ ConstraintMatrix::condense (const SparseMatrix<number> &uncondensed,
                                   // row of the inhomogeneous constraint in
                                   // the matrix with Gauss elimination
              if (use_vectors == true)
-               condensed_vector(new_line[row]) -= uncondensed.global_entry(j) * 
+               condensed_vector(new_line[row]) -= uncondensed.global_entry(j) *
                                                   c->inhomogeneity;
            }
 
        if (use_vectors == true)
-         condensed_vector(new_line[row]) += uncondensed_vector(row);     
+         condensed_vector(new_line[row]) += uncondensed_vector(row);
       }
     else
                                       // line must be distributed
@@ -292,12 +293,12 @@ ConstraintMatrix::condense (const SparseMatrix<number> &uncondensed,
                                           // for each column: distribute
          if (new_line[uncondensed_struct.get_column_numbers()[j]] != -1)
                                             // column is not constrained
-           for (unsigned int q=0; q!=next_constraint->entries.size(); ++q) 
+           for (unsigned int q=0; q!=next_constraint->entries.size(); ++q)
              condensed.add (new_line[next_constraint->entries[q].first],
                             new_line[uncondensed_struct.get_column_numbers()[j]],
                             uncondensed.global_entry(j) *
                             next_constraint->entries[q].second);
-       
+
          else
                                             // not only this line but
                                             // also this col is constrained
@@ -307,7 +308,7 @@ ConstraintMatrix::condense (const SparseMatrix<number> &uncondensed,
              std::vector<ConstraintLine>::const_iterator c = lines.begin();
              while (c->line != uncondensed_struct.get_column_numbers()[j])
                ++c;
-             
+
              for (unsigned int p=0; p!=c->entries.size(); ++p)
                for (unsigned int q=0; q!=next_constraint->entries.size(); ++q)
                  condensed.add (new_line[next_constraint->entries[q].first],
@@ -319,14 +320,14 @@ ConstraintMatrix::condense (const SparseMatrix<number> &uncondensed,
              if (use_vectors == true)
                for (unsigned int q=0; q!=next_constraint->entries.size(); ++q)
                  condensed_vector (new_line[next_constraint->entries[q].first])
-                   -= uncondensed.global_entry(j) * 
+                   -= uncondensed.global_entry(j) *
                       next_constraint->entries[q].second *
                       c->inhomogeneity;
            };
 
                                   // condense the vector
        if (use_vectors == true)
-         for (unsigned int q=0; q!=next_constraint->entries.size(); ++q) 
+         for (unsigned int q=0; q!=next_constraint->entries.size(); ++q)
            condensed_vector(new_line[next_constraint->entries[q].first])
              +=
              uncondensed_vector(row) * next_constraint->entries[q].second;
@@ -355,7 +356,7 @@ ConstraintMatrix::condense (SparseMatrix<number> &uncondensed,
          ExcNotQuadratic());
   if (use_vectors == true)
     {
-      Assert (vec.size() == sparsity.n_rows(), 
+      Assert (vec.size() == sparsity.n_rows(),
              ExcDimensionMismatch(vec.size(), sparsity.n_rows()));
     }
 
@@ -363,7 +364,7 @@ ConstraintMatrix::condense (SparseMatrix<number> &uncondensed,
   for (unsigned int i=0; i<uncondensed.m(); ++i)
     average_diagonal += std::fabs (uncondensed.diag_element(i));
   average_diagonal /= uncondensed.m();
-  
+
                                   // store for each index whether it must be
                                   // distributed or not. If entry is
                                   // invalid_unsigned_int, no distribution is
@@ -372,7 +373,7 @@ ConstraintMatrix::condense (SparseMatrix<number> &uncondensed,
                                   // handles this index
   std::vector<unsigned int> distribute (sparsity.n_rows(),
                                         numbers::invalid_unsigned_int);
-  
+
   for (unsigned int c=0; c<lines.size(); ++c)
     distribute[lines[c].line] = c;
 
@@ -387,7 +388,7 @@ ConstraintMatrix::condense (SparseMatrix<number> &uncondensed,
                entry != uncondensed.end(row); ++entry)
             {
               const unsigned int column = entry->column();
-              
+
                                                // end of row reached?
                                                // this should not
                                                // happen, since we only
@@ -395,7 +396,7 @@ ConstraintMatrix::condense (SparseMatrix<number> &uncondensed,
                                                // matrices!
               Assert (column != SparsityPattern::invalid_entry,
                       ExcMatrixNotClosed());
-           
+
               if (distribute[column] != numbers::invalid_unsigned_int)
                                                  // distribute entry at
                                                  // regular row @p row
@@ -417,7 +418,7 @@ ConstraintMatrix::condense (SparseMatrix<number> &uncondensed,
                                   // row of the inhomogeneous constraint in
                                   // the matrix with Gauss elimination
                  if (use_vectors == true)
-                   vec(row) -= 
+                   vec(row) -=
                      entry->value() * lines[distribute[column]].inhomogeneity;
 
                                                    // set old value to zero
@@ -451,7 +452,7 @@ ConstraintMatrix::condense (SparseMatrix<number> &uncondensed,
                                                  // old entry to zero
                 {
                   for (unsigned int q=0;
-                       q!=lines[distribute[row]].entries.size(); ++q) 
+                       q!=lines[distribute[row]].entries.size(); ++q)
                     uncondensed.add (lines[distribute[row]].entries[q].first,
                                      column,
                                      entry->value() *
@@ -479,7 +480,7 @@ ConstraintMatrix::condense (SparseMatrix<number> &uncondensed,
                                         lines[distribute[column]].entries[q].second);
 
                      if (use_vectors == true)
-                       vec(lines[distribute[row]].entries[p].first) -= 
+                       vec(lines[distribute[row]].entries[p].first) -=
                          entry->value() * lines[distribute[row]].entries[p].second *
                          lines[distribute[column]].inhomogeneity;
                    }
@@ -493,7 +494,7 @@ ConstraintMatrix::condense (SparseMatrix<number> &uncondensed,
                                   // take care of vector
          if (use_vectors == true)
            {
-             for (unsigned int q=0; q!=lines[distribute[row]].entries.size(); ++q) 
+             for (unsigned int q=0; q!=lines[distribute[row]].entries.size(); ++q)
                vec(lines[distribute[row]].entries[q].first)
                  += (vec(row) * lines[distribute[row]].entries[q].second);
 
@@ -516,7 +517,7 @@ ConstraintMatrix::condense (BlockSparseMatrix<number> &uncondensed,
   const bool use_vectors = vec.n_blocks() == 0 ? false : true;
 
   const unsigned int blocks = uncondensed.n_block_rows();
-  
+
   const BlockSparsityPattern &
     sparsity = uncondensed.get_sparsity_pattern ();
 
@@ -533,7 +534,7 @@ ConstraintMatrix::condense (BlockSparseMatrix<number> &uncondensed,
 
   if (use_vectors == true)
     {
-      Assert (vec.size() == sparsity.n_rows(), 
+      Assert (vec.size() == sparsity.n_rows(),
              ExcDimensionMismatch(vec.size(), sparsity.n_rows()));
       Assert (vec.n_blocks() == sparsity.n_block_rows(),
              ExcDimensionMismatch(vec.n_blocks(), sparsity.n_block_rows()));
@@ -547,7 +548,7 @@ ConstraintMatrix::condense (BlockSparseMatrix<number> &uncondensed,
 
   const BlockIndices &
     index_mapping = sparsity.get_column_indices();
-  
+
                                   // store for each index whether it must be
                                   // distributed or not. If entry is
                                   // numbers::invalid_unsigned_int,
@@ -557,7 +558,7 @@ ConstraintMatrix::condense (BlockSparseMatrix<number> &uncondensed,
                                   // index
   std::vector<unsigned int> distribute (sparsity.n_rows(),
                                         numbers::invalid_unsigned_int);
-  
+
   for (unsigned int c=0; c<lines.size(); ++c)
     distribute[lines[c].line] = c;
 
@@ -569,7 +570,7 @@ ConstraintMatrix::condense (BlockSparseMatrix<number> &uncondensed,
       const std::pair<unsigned int,unsigned int>
        block_index = index_mapping.global_to_local(row);
       const unsigned int block_row = block_index.first;
-      
+
       if (distribute[row] == numbers::invalid_unsigned_int)
                                         // regular line. loop over
                                         // all columns and see
@@ -592,7 +593,7 @@ ConstraintMatrix::condense (BlockSparseMatrix<number> &uncondensed,
                 {
                   const unsigned int global_col
                     = index_mapping.local_to_global(block_col,entry->column());
-                   
+
                   if (distribute[global_col] != numbers::invalid_unsigned_int)
                                                      // distribute entry at
                                                      // regular row @p row
@@ -601,7 +602,7 @@ ConstraintMatrix::condense (BlockSparseMatrix<number> &uncondensed,
                                                      // entry to zero
                     {
                       const double old_value = entry->value ();
-                       
+
                       for (unsigned int q=0;
                            q!=lines[distribute[global_col]].entries.size(); ++q)
                         uncondensed.add (row,
@@ -615,7 +616,7 @@ ConstraintMatrix::condense (BlockSparseMatrix<number> &uncondensed,
                                   // row of the inhomogeneous constraint in
                                   // the matrix with Gauss elimination
                      if (use_vectors == true)
-                       vec(row) -= entry->value() * 
+                       vec(row) -= entry->value() *
                                    lines[distribute[global_col]].inhomogeneity;
 
                       entry->value() = 0.;
@@ -639,7 +640,7 @@ ConstraintMatrix::condense (BlockSparseMatrix<number> &uncondensed,
                 {
                   const unsigned int global_col
                     = index_mapping.local_to_global (block_col, entry->column());
-                   
+
                   if (distribute[global_col] ==
                       numbers::invalid_unsigned_int)
                                                      // distribute
@@ -653,9 +654,9 @@ ConstraintMatrix::condense (BlockSparseMatrix<number> &uncondensed,
                                                      // zero
                     {
                       const double old_value = entry->value();
-                         
+
                       for (unsigned int q=0;
-                           q!=lines[distribute[row]].entries.size(); ++q) 
+                           q!=lines[distribute[row]].entries.size(); ++q)
                         uncondensed.add (lines[distribute[row]].entries[q].first,
                                          global_col,
                                          old_value *
@@ -673,7 +674,7 @@ ConstraintMatrix::condense (BlockSparseMatrix<number> &uncondensed,
                                                      // otherwise
                     {
                       const double old_value = entry->value ();
-                         
+
                       for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
                        {
                          for (unsigned int q=0; q!=lines[distribute[global_col]].entries.size(); ++q)
@@ -684,7 +685,7 @@ ConstraintMatrix::condense (BlockSparseMatrix<number> &uncondensed,
                                             lines[distribute[global_col]].entries[q].second);
 
                          if (use_vectors == true)
-                           vec(lines[distribute[row]].entries[p].first) -= 
+                           vec(lines[distribute[row]].entries[p].first) -=
                              old_value * lines[distribute[row]].entries[p].second *
                              lines[distribute[global_col]].inhomogeneity;
                        }
@@ -697,7 +698,7 @@ ConstraintMatrix::condense (BlockSparseMatrix<number> &uncondensed,
                                           // take care of vector
          if (use_vectors == true)
            {
-             for (unsigned int q=0; q!=lines[distribute[row]].entries.size(); ++q) 
+             for (unsigned int q=0; q!=lines[distribute[row]].entries.size(); ++q)
                vec(lines[distribute[row]].entries[q].first)
                  += (vec(row) * lines[distribute[row]].entries[q].second);
 
@@ -749,7 +750,7 @@ distribute_local_to_global (const Vector<double>            &local_vector,
                                   // vector that tells about whether a
                                   // certain constraint exists. then, we
                                   // simply copy over the data.
-      const std::vector<ConstraintLine>::const_iterator position = 
+      const std::vector<ConstraintLine>::const_iterator position =
        find_constraint(local_dof_indices[i]);
 
       if (position==lines.end())
@@ -764,7 +765,7 @@ distribute_local_to_global (const Vector<double>            &local_vector,
          if (val != 0)
            for (unsigned int j=0; j<n_local_dofs; ++j)
              {
-               const std::vector<ConstraintLine>::const_iterator 
+               const std::vector<ConstraintLine>::const_iterator
                  position_j = find_constraint(local_dof_indices[j]);
 
                if (position_j == lines.end())
@@ -821,7 +822,7 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                                   // feature further down.
   Vector<double> dummy(0);
   distribute_local_to_global (local_matrix, dummy, local_dof_indices,
-                             global_matrix, dummy, 
+                             global_matrix, dummy,
                              internal::bool2type<IsBlockMatrix<MatrixType>::value>());
 }
 
@@ -841,7 +842,7 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                                   // actual implementation follows further
                                   // down.
   distribute_local_to_global (local_matrix, local_vector, local_dof_indices,
-                             global_matrix, global_vector, 
+                             global_matrix, global_vector,
                              internal::bool2type<IsBlockMatrix<MatrixType>::value>());
 }
 
@@ -888,7 +889,7 @@ ConstraintMatrix::distribute (const VectorType &condensed,
   unsigned int                                shift           = 0;
   unsigned int n_rows = uncondensed.size();
 
-  if (next_constraint == lines.end()) 
+  if (next_constraint == lines.end())
                                     // if no constraint is to be handled
     for (unsigned int row=0; row!=n_rows; ++row)
       old_line.push_back (row);
@@ -921,7 +922,7 @@ ConstraintMatrix::distribute (const VectorType &condensed,
                                   // only evaluated so often as there are
                                   // entries in new_line[*] which tells us
                                   // which constraints exist
-  for (unsigned int line=0; line<uncondensed.size(); ++line) 
+  for (unsigned int line=0; line<uncondensed.size(); ++line)
     if (old_line[line] != -1)
                                       // line was not condensed away
       uncondensed(line) = condensed(old_line[line]);
@@ -949,14 +950,16 @@ ConstraintMatrix::distribute (VectorType &vec) const
   Assert (sorted == true, ExcMatrixNotClosed());
 
   std::vector<ConstraintLine>::const_iterator next_constraint = lines.begin();
-  for (; next_constraint != lines.end(); ++next_constraint) 
+  for (; next_constraint != lines.end(); ++next_constraint)
     {
                                       // fill entry in line
                                       // next_constraint.line by adding the
                                       // different contributions
-      double new_value = next_constraint->inhomogeneity;
+      typename VectorType::value_type
+       new_value = next_constraint->inhomogeneity;
       for (unsigned int i=0; i<next_constraint->entries.size(); ++i)
-       new_value += (vec(next_constraint->entries[i].first) *
+       new_value += (static_cast<typename VectorType::value_type>
+                     (vec(next_constraint->entries[i].first)) *
                       next_constraint->entries[i].second);
       vec(next_constraint->line) = new_value;
     }
@@ -972,7 +975,7 @@ namespace internals
                                   // which global entries (global_row) are
                                   // given rise by local entries
                                   // (local_row) or some constraints.
-  struct distributing 
+  struct distributing
   {
     distributing (const unsigned int global_row = deal_II_numbers::invalid_unsigned_int,
                  const unsigned int local_row = deal_II_numbers::invalid_unsigned_int);
@@ -1070,7 +1073,7 @@ namespace internals
       }
 
     if (&*pos1->constraints == 0)
-      pos1->constraints = 
+      pos1->constraints =
        new std::vector<std::pair<unsigned int,double> > (1,constraint);
     else
       pos1->constraints->push_back (constraint);
@@ -1113,7 +1116,7 @@ namespace internals
            j2 = j-istep;
            temp = my_indices[i].global_row;
            templ = my_indices[i].local_row;
-           if (my_indices[j2].global_row > temp) 
+           if (my_indices[j2].global_row > temp)
              {
                while ((j >= istep) && (my_indices[j2].global_row > temp))
                  {
@@ -1154,7 +1157,7 @@ namespace internals
     block_starts[0] = 0;
     for (unsigned int i=1;i<num_blocks;++i)
       {
-       row_iterator first_block = 
+       row_iterator first_block =
          std::lower_bound (col_indices,
                            row_indices.end(),
                            block_object.get_row_indices().block_start(i));
@@ -1188,7 +1191,7 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                                   // check whether we work on real vectors
                                   // or we just used a dummy when calling
                                   // the other function above.
-  const bool use_vectors = (local_vector.size() == 0 && 
+  const bool use_vectors = (local_vector.size() == 0 &&
                            global_vector.size() == 0) ? false : true;
 
   Assert (local_matrix.n() == local_dof_indices.size(),
@@ -1284,7 +1287,7 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
        {
          have_indirect_rows = true;
          internals::insert_index(my_indices, position->entries[q].first,
-                                 std::make_pair<unsigned int,double> 
+                                 std::make_pair<unsigned int,double>
                                  (local_row, position->entries[q].second));
        }
 
@@ -1313,7 +1316,7 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                                   // hanging nodes in 3d). however, in the
                                   // line below, we do actually do
                                   // something with this dof
-      const typename MatrixType::value_type new_diagonal 
+      const typename MatrixType::value_type new_diagonal
        = (std::fabs(local_matrix(local_row,local_row)) != 0 ?
           std::fabs(local_matrix(local_row,local_row)) : average_diagonal);
       global_matrix.add(global_row, global_row, new_diagonal);
@@ -1329,7 +1332,7 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
 
   typedef std::vector<std::pair<unsigned int,double> > constraint_format;
 
-                                  // now do the actual job. 
+                                  // now do the actual job.
   for (unsigned int i=0; i<n_actual_dofs; ++i)
     {
       const unsigned int row = my_indices[i].global_row;
@@ -1343,7 +1346,7 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                                   // all on this set of dofs (saves a lot
                                   // of checks). the only check we actually
                                   // need to perform is whether the matrix
-                                  // element is zero. 
+                                  // element is zero.
       if (have_indirect_rows == false)
        {
          Assert(loc_row < n_local_dofs, ExcInternalError());
@@ -1416,7 +1419,7 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                    {
                      constraint_format &constraint_j = *my_indices[j].constraints;
                      for (unsigned int p=0; p<constraint_j.size(); ++p)
-                       col_val += matrix_ptr[constraint_j[p].first] 
+                       col_val += matrix_ptr[constraint_j[p].first]
                                   *
                                   constraint_j[p].second;
                    }
@@ -1440,13 +1443,13 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                    {
                      double add_this = loc_col != deal_II_numbers::invalid_unsigned_int ?
                        local_matrix(constraint_i[q].first, loc_col) : 0;
-  
+
                      if (my_indices[j].constraints != 0)
                        {
                          constraint_format &constraint_j = *my_indices[j].constraints;
                          for (unsigned int p=0; p<constraint_j.size(); ++p)
                            add_this += local_matrix(constraint_i[q].first,
-                                                    constraint_j[p].first) 
+                                                    constraint_j[p].first)
                                        *
                                        constraint_j[p].second;
                        }
@@ -1454,7 +1457,7 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                    }
                }
 
-                  
+
 
                                   // if we got some nontrivial value,
                                   // append it to the array of values.
@@ -1488,7 +1491,7 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
 
              if (my_indices[i].constraints != 0)
                {
-                 std::vector<std::pair<unsigned int,double> > &constraint_i = 
+                 std::vector<std::pair<unsigned int,double> > &constraint_i =
                    *my_indices[i].constraints;
 
                  for (unsigned int q=0; q<constraint_i.size(); ++q)
@@ -1536,7 +1539,7 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                                   // the other function for additional
                                   // comments.
 
-  const bool use_vectors = (local_vector.size() == 0 && 
+  const bool use_vectors = (local_vector.size() == 0 &&
                            global_vector.size() == 0) ? false : true;
 
   Assert (local_matrix.n() == local_dof_indices.size(),
@@ -1600,11 +1603,11 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
        {
          have_indirect_rows = true;
          internals::insert_index(my_indices, position->entries[q].first,
-                                 std::make_pair<unsigned int,double> 
+                                 std::make_pair<unsigned int,double>
                                  (local_row, position->entries[q].second));
        }
 
-      const typename MatrixType::value_type new_diagonal 
+      const typename MatrixType::value_type new_diagonal
        = (std::fabs(local_matrix(local_row,local_row)) != 0 ?
           std::fabs(local_matrix(local_row,local_row)) : average_diagonal);
       global_matrix.add(global_row, global_row, new_diagonal);
@@ -1684,12 +1687,12 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                                   // constraints
                          if (my_indices[j].constraints != 0)
                            {
-                             constraint_format &constraint_j = 
+                             constraint_format &constraint_j =
                                *my_indices[j].constraints;
 
                              for (unsigned int p=0; p<constraint_j.size(); ++p)
                                col_val += local_matrix(loc_row,
-                                                       constraint_j[p].first) 
+                                                       constraint_j[p].first)
                                           *
                                           constraint_j[p].second;
                            }
@@ -1704,18 +1707,18 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
 
                          for (unsigned int q=0; q<constraint_i.size(); ++q)
                            {
-                             double add_this = 
+                             double add_this =
                                loc_col != deal_II_numbers::invalid_unsigned_int ?
                                local_matrix(constraint_i[q].first, loc_col) : 0;
-  
+
                              if (my_indices[j].constraints != 0)
                                {
-                                 constraint_format &constraint_j = 
+                                 constraint_format &constraint_j =
                                    *my_indices[j].constraints;
 
                                  for (unsigned int p=0; p<constraint_j.size(); ++p)
                                    add_this += local_matrix(constraint_i[q].first,
-                                                            constraint_j[p].first) 
+                                                            constraint_j[p].first)
                                                *
                                                constraint_j[p].second;
                                }
@@ -1726,7 +1729,7 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                      if (col_val != 0)
                        {
                          *col_ptr++ = localized_indices[j];
-                         *val_ptr++ = 
+                         *val_ptr++ =
                            static_cast<typename MatrixType::value_type>(col_val);
                        }
                    }
@@ -1762,7 +1765,7 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
 
              if (my_indices[i].constraints != 0)
                {
-                 std::vector<std::pair<unsigned int,double> > &constraint_i = 
+                 std::vector<std::pair<unsigned int,double> > &constraint_i =
                    *my_indices[i].constraints;
 
                  for (unsigned int q=0; q<constraint_i.size(); ++q)
@@ -1776,7 +1779,7 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                    }
                }
              if (val != 0)
-               global_vector(my_indices[i].global_row) += 
+               global_vector(my_indices[i].global_row) +=
                  static_cast<typename VectorType::value_type>(val);
            }
        }
@@ -1854,7 +1857,7 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
                }
              else
                {
-                 std::vector<unsigned int>::iterator it = 
+                 std::vector<unsigned int>::iterator it =
                    std::lower_bound(actual_dof_indices.begin(),
                                     actual_dof_indices.end(),
                                     new_index);
@@ -1943,7 +1946,7 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
        {
          have_indirect_rows = true;
          internals::insert_index(my_indices, position->entries[q].first,
-                                 std::make_pair<unsigned int,double> 
+                                 std::make_pair<unsigned int,double>
                                  (local_row, position->entries[q].second));
        }
 
@@ -2032,7 +2035,7 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
                                   // constraints
                  if (my_indices[j].constraints != 0)
                    {
-                     std::vector<std::pair<unsigned int,double> > &constraint_j = 
+                     std::vector<std::pair<unsigned int,double> > &constraint_j =
                        *my_indices[j].constraints;
 
                      for (unsigned int p=0; p<constraint_j.size(); ++p)
@@ -2047,7 +2050,7 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
                                   // given column.
              if (my_indices[i].constraints != 0)
                {
-                 std::vector<std::pair<unsigned int,double> > &constraint_i = 
+                 std::vector<std::pair<unsigned int,double> > &constraint_i =
                    *my_indices[i].constraints;
                  for (unsigned int q=0; q<constraint_i.size(); ++q)
                    {
@@ -2057,10 +2060,10 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
                          if (dof_mask[constraint_i[q].first][loc_col] == true)
                            goto add_this_index;
                        }
-  
+
                      if (my_indices[j].constraints != 0)
                        {
-                         std::vector<std::pair<unsigned int,double> > &constraint_j = 
+                         std::vector<std::pair<unsigned int,double> > &constraint_j =
                            *my_indices[j].constraints;
 
                          for (unsigned int p=0; p<constraint_j.size(); ++p)
@@ -2171,7 +2174,7 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
                }
              else
                {
-                 std::vector<unsigned int>::iterator it = 
+                 std::vector<unsigned int>::iterator it =
                    std::lower_bound(actual_dof_indices.begin(),
                                     actual_dof_indices.end(),
                                     new_index);
@@ -2220,8 +2223,8 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
                  const unsigned int next_block_col = block_starts[block_col+1];
                  sparsity_pattern.block(block,block_col).
                    add_entries(row,
-                               index_it, 
-                               actual_dof_indices.begin() + next_block_col, 
+                               index_it,
+                               actual_dof_indices.begin() + next_block_col,
                                true);
                  index_it = actual_dof_indices.begin() + next_block_col;
                }
@@ -2288,7 +2291,7 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
        {
          have_indirect_rows = true;
          internals::insert_index(my_indices, position->entries[q].first,
-                                 std::make_pair<unsigned int,double> 
+                                 std::make_pair<unsigned int,double>
                                  (local_row, position->entries[q].second));
        }
 
@@ -2384,7 +2387,7 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
                                   // constraints
                          if (my_indices[j].constraints != 0)
                            {
-                             std::vector<std::pair<unsigned int,double> > 
+                             std::vector<std::pair<unsigned int,double> >
                                &constraint_j = *my_indices[j].constraints;
 
                              for (unsigned int p=0; p<constraint_j.size(); ++p)
@@ -2399,7 +2402,7 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
                                   // given column.
                      if (my_indices[i].constraints != 0)
                        {
-                         std::vector<std::pair<unsigned int,double> > 
+                         std::vector<std::pair<unsigned int,double> >
                            &constraint_i = *my_indices[i].constraints;
                          for (unsigned int q=0; q<constraint_i.size(); ++q)
                            {
@@ -2410,10 +2413,10 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
                                  if (dof_mask[constraint_i[q].first][loc_col] == true)
                                    goto add_this_index;
                                }
-  
+
                              if (my_indices[j].constraints != 0)
                                {
-                                 std::vector<std::pair<unsigned int,double> > 
+                                 std::vector<std::pair<unsigned int,double> >
                                    &constraint_j = *my_indices[j].constraints;
 
                                  for (unsigned int p=0; p<constraint_j.size(); ++p)
@@ -2435,8 +2438,8 @@ add_entries_local_to_global (const std::vector<unsigned int> &local_dof_indices,
                                   // that accumulated under the given
                                   // process into the global matrix row and
                                   // into the vector
-             sparsity_pattern.block(block, block_col).add_entries(row, 
-                                                                  cols.begin(), 
+             sparsity_pattern.block(block, block_col).add_entries(row,
+                                                                  cols.begin(),
                                                                   col_ptr,
                                                                   true);
            }
index 592b81777cef58bb32ec8ede53865c3106aac062..17121547a0aeb7a8b3402c48b3ea19c7a9c75327 100644 (file)
@@ -323,7 +323,7 @@ namespace PETScWrappers
                                         * previously used.
                                         */
       MatrixBase &
-      operator = (const double d);
+      operator = (const value_type d);
                                        /**
                                         * Release all memory and return
                                         * to a state just like after
index 38128e8440517c0ca8eee47bc9eb9064015c053c..c58591aea41cd2a616db6b98b90c4fac0c550cca 100644 (file)
@@ -34,7 +34,7 @@ namespace PETScWrappers
 {
   namespace MPI
   {
-    
+
 /**
  * Implementation of a parallel sparse matrix class based on PETSC, with rows
  * of the matrix distributed across an MPI network. All the functionality is
@@ -104,7 +104,7 @@ namespace PETScWrappers
  * which. Thus, it is not sufficient to just count the number of degrees of
  * freedom that belong to a particular process, but you have to have the
  * numbers for all processes available at all processes.
- * 
+ *
  * @ingroup PETScWrappers
  * @ingroup Matrix1
  * @author Wolfgang Bangerth, 2004
@@ -147,7 +147,7 @@ namespace PETScWrappers
                                           * empty matrix.
                                           */
         SparseMatrix ();
-      
+
                                          /**
                                           * Create a sparse matrix of
                                           * dimensions @p m times @p n, with
@@ -164,7 +164,7 @@ namespace PETScWrappers
                                           * local_row and @p local_columns
                                           * parameters, see the class
                                           * documentation.
-                                          * 
+                                          *
                                           * The @p is_symmetric flag determines
                                           * whether we should tell PETSc that
                                           * the matrix is going to be symmetric
@@ -197,7 +197,7 @@ namespace PETScWrappers
                                           * local_row and @p local_columns
                                           * parameters, see the class
                                           * documentation.
-                                          * 
+                                          *
                                           * Just as for the other
                                           * constructors: PETSc is able to
                                           * cope with the situation that more
@@ -226,7 +226,7 @@ namespace PETScWrappers
                       const unsigned int               local_columns,
                       const std::vector<unsigned int> &row_lengths,
                       const bool                       is_symmetric = false);
-        
+
                                          /**
                                           * Initialize using the given
                                           * sparsity pattern with
@@ -238,7 +238,7 @@ namespace PETScWrappers
                                           * local_columns_per_process
                                           * parameters, see the class
                                           * documentation.
-                                          * 
+                                          *
                                           * Note that PETSc can be very slow
                                           * if you do not provide it with a
                                           * good estimate of the lengths of
@@ -291,7 +291,7 @@ namespace PETScWrappers
                                           * but keep the sparsity pattern
                                           * previously used.
                                           */
-        SparseMatrix & operator = (const double d);
+        SparseMatrix & operator = (const value_type d);
 
                                          /**
                                           * Throw away the present matrix and
@@ -366,14 +366,14 @@ namespace PETScWrappers
                      const std::vector<unsigned int> &local_columns_per_process,
                      const unsigned int               this_process,
                      const bool                       preset_nonzero_locations = true);
-        
+
                                          /**
                                           * Return a reference to the MPI
                                           * communicator object in use with
                                           * this matrix.
                                           */
         virtual const MPI_Comm & get_mpi_communicator () const;
-        
+
                                     /** @addtogroup Exceptions
                                      * @{ */
                                         /**
@@ -383,7 +383,7 @@ namespace PETScWrappers
                        int, int,
                        << "The number of local rows " << arg1
                        << " must be larger than the total number of rows " << arg2);
-                                        //@}        
+                                        //@}
       private:
 
                                          /**
@@ -431,7 +431,7 @@ namespace PETScWrappers
                                        *  To allow calling protected
                                        *  prepare_add() and
                                        *  prepare_set().
-                                       */  
+                                       */
        friend class BlockMatrixBase<SparseMatrix>;
     };
 
@@ -440,7 +440,7 @@ namespace PETScWrappers
 // -------- template and inline functions ----------
 
     inline
-    const MPI_Comm &    
+    const MPI_Comm &
     SparseMatrix::get_mpi_communicator () const
     {
       return communicator;
index 380652cfe42e1f8134a8dd98a3f00ac4c630db6e..45deb8ea3be54ca23bdde7c19bf5fa58f77e4c30 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004, 2005, 2006, 2007 by the deal.II authors
+//    Copyright (C) 2004, 2005, 2006, 2007, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -32,8 +32,8 @@ namespace PETScWrappers
   class MatrixBase;
   class VectorBase;
   class PreconditionerBase;
-  
-  
+
+
 /**
  * Base class for solver classes using the PETSc solvers. Since solvers in
  * PETSc are selected based on flags passed to a generic solver object,
@@ -62,7 +62,7 @@ namespace PETScWrappers
  * fact that the default argument to the solver classes, @p PETSC_COMM_SELF,
  * is the appropriate argument for the sequential case (which is why it is the
  * default argument), so this error only shows up in parallel mode.
- * 
+ *
  * @ingroup PETScWrappers
  * @author Wolfgang Bangerth, 2004
  */
@@ -87,7 +87,7 @@ namespace PETScWrappers
                                         */
       SolverBase (SolverControl  &cn,
                   const MPI_Comm &mpi_communicator);
-      
+
                                        /**
                                         * Destructor.
                                         */
@@ -165,7 +165,7 @@ namespace PETScWrappers
       int
       convergence_test (KSP                 ksp,
                         const int           iteration,
-                        const PetscScalar   residual_norm,
+                        const PetscReal     residual_norm,
                         KSPConvergedReason *reason,
                         void               *solver_control);
 
@@ -204,7 +204,7 @@ namespace PETScWrappers
                                             * Destructor
                                             */
           ~SolverData ();
-          
+
 #if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR < 2)
                                            /**
                                             * A PETSc solver object.
@@ -245,7 +245,7 @@ namespace PETScWrappers
                                         * pipe additional data to the
                                         * solver.
                                         */
-      struct AdditionalData 
+      struct AdditionalData
       {
                                            /**
                                             * Constructor. By default,
@@ -253,13 +253,13 @@ namespace PETScWrappers
                                             * to one.
                                             */
           AdditionalData (const double omega = 1);
-       
+
                                            /**
                                             * Relaxation parameter.
                                             */
           double omega;
       };
-       
+
                                        /**
                                         * Constructor. In contrast to
                                         * deal.II's own solvers, there is no
@@ -327,7 +327,7 @@ namespace PETScWrappers
                                         */
       struct AdditionalData
       {};
-       
+
                                        /**
                                         * Constructor. In contrast to
                                         * deal.II's own solvers, there is no
@@ -395,7 +395,7 @@ namespace PETScWrappers
                                         */
       struct AdditionalData
       {};
-       
+
                                        /**
                                         * Constructor. In contrast to
                                         * deal.II's own solvers, there is no
@@ -463,7 +463,7 @@ namespace PETScWrappers
                                         */
       struct AdditionalData
       {};
-       
+
                                        /**
                                         * Constructor. In contrast to
                                         * deal.II's own solvers, there is no
@@ -538,7 +538,7 @@ namespace PETScWrappers
                                             */
           AdditionalData (const unsigned int restart_parameter = 30,
                          const bool right_preconditioning = false);
-       
+
                                            /**
                                             * Maximum number of
                                             * tmp vectors.
@@ -551,7 +551,7 @@ namespace PETScWrappers
                                            */
          bool right_preconditioning;
       };
-       
+
                                        /**
                                         * Constructor. In contrast to
                                         * deal.II's own solvers, there is no
@@ -619,7 +619,7 @@ namespace PETScWrappers
                                         */
       struct AdditionalData
       {};
-       
+
                                        /**
                                         * Constructor. In contrast to
                                         * deal.II's own solvers, there is no
@@ -685,7 +685,7 @@ namespace PETScWrappers
                                         */
       struct AdditionalData
       {};
-       
+
                                        /**
                                         * Constructor. In contrast to
                                         * deal.II's own solvers, there is no
@@ -753,7 +753,7 @@ namespace PETScWrappers
                                         */
       struct AdditionalData
       {};
-       
+
                                        /**
                                         * Constructor. In contrast to
                                         * deal.II's own solvers, there is no
@@ -826,7 +826,7 @@ namespace PETScWrappers
                                         */
       struct AdditionalData
       {};
-       
+
                                        /**
                                         * Constructor. In contrast to
                                         * deal.II's own solvers, there is no
@@ -894,7 +894,7 @@ namespace PETScWrappers
                                         */
       struct AdditionalData
       {};
-       
+
                                        /**
                                         * Constructor. In contrast to
                                         * deal.II's own solvers, there is no
@@ -962,7 +962,7 @@ namespace PETScWrappers
                                         */
       struct AdditionalData
       {};
-       
+
                                        /**
                                         * Constructor. In contrast to
                                         * deal.II's own solvers, there is no
@@ -1010,14 +1010,14 @@ namespace PETScWrappers
                                         */
       virtual void set_solver_type (KSP &ksp) const;
   };
-  
+
 
 /**
  * An implementation of the solver interface using the PETSc PREONLY
  * solver. Actually this is NOT a real solution algorithm. Its only
  * purpose is to provide a solver object, when the preconditioner
  * should be used as real solver. It is very useful in conjunction with
- * the complete LU decomposition preconditioner <tt> PreconditionLU </tt>, 
+ * the complete LU decomposition preconditioner <tt> PreconditionLU </tt>,
  * which in conjunction with this solver class becomes a direct solver.
  *
  * @ingroup PETScWrappers
@@ -1033,7 +1033,7 @@ namespace PETScWrappers
                                         */
       struct AdditionalData
       {};
-       
+
                                        /**
                                         * Constructor. In contrast to
                                         * deal.II's own solvers, there is no
@@ -1065,7 +1065,7 @@ namespace PETScWrappers
       SolverPreOnly (SolverControl        &cn,
                     const MPI_Comm       &mpi_communicator = PETSC_COMM_SELF,
                     const AdditionalData &data = AdditionalData());
-      
+
     protected:
                                        /**
                                         * Store a copy of the flags for this
index 5cf061e7cd822154c7d0b521305bcca30a911e32..b2f920801fe63da77ed99e232cbc362764bfddcc 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004, 2005, 2006, 2007 by the deal.II authors
+//    Copyright (C) 2004, 2005, 2006, 2007, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -47,7 +47,7 @@ namespace PETScWrappers
                                   /**
                                    * @cond internal
                                    */
-  
+
 /**
  * A namespace for internal implementation details of the PETScWrapper
  * members.
@@ -84,7 +84,7 @@ namespace PETScWrappers
                                           */
         VectorReference (const VectorBase  &vector,
                          const unsigned int index);
-        
+
       public:
                                          /**
                                           * This looks like a copy operator,
@@ -107,7 +107,7 @@ namespace PETScWrappers
                                           * the assignment operator const.
                                           */
         const VectorReference & operator = (const VectorReference &r) const;
-        
+
                                          /**
                                           * Set the referenced element of the
                                           * vector to <tt>s</tt>.
@@ -162,7 +162,7 @@ namespace PETScWrappers
                         << " of a distributed vector, but only elements "
                         << arg2 << " through " << arg3
                         << " are stored locally and can be accessed.");
-        
+
       private:
                                          /**
                                           * Point to the vector we are
@@ -187,7 +187,7 @@ namespace PETScWrappers
                                    /**
                                     * @endcond
                                     */
-  
+
 
 /**
  * Base class for all vector classes that are implemented on top of the PETSc
@@ -224,7 +224,7 @@ namespace PETScWrappers
                                         * class.
                                         */
       typedef PetscScalar               value_type;
-      typedef PetscScalar               real_type;
+      typedef PetscReal                 real_type;
       typedef size_t                    size_type;
       typedef internal::VectorReference reference;
       typedef const internal::VectorReference const_reference;
@@ -235,7 +235,7 @@ namespace PETScWrappers
                                         * to initialize the data.
                                         */
       VectorBase ();
-      
+
                                        /**
                                         * Copy constructor. Sets the dimension
                                         * to that of the given vector, and
@@ -279,8 +279,8 @@ namespace PETScWrappers
                                         * other values is deprecated and may
                                         * be disallowed in the future.
                                         */
-      VectorBase & operator = (const PetscScalar s);    
-      
+      VectorBase & operator = (const PetscScalar s);
+
                                        /**
                                         * Test for equality. This function
                                         * assumes that the present vector and
@@ -290,7 +290,7 @@ namespace PETScWrappers
                                         * sense anyway.
                                         */
       bool operator == (const VectorBase &v) const;
-    
+
                                        /**
                                         * Test for inequality. This function
                                         * assumes that the present vector and
@@ -383,7 +383,7 @@ namespace PETScWrappers
                                         */
       void add (const std::vector<unsigned int> &indices,
                 const std::vector<PetscScalar>  &values);
-      
+
                                        /**
                                         * Return the scalar product of two
                                         * vectors. The vectors must have the
@@ -421,7 +421,7 @@ namespace PETScWrappers
                                         * powers of the absolute values
                                         * of the elements.
                                         */
-      real_type lp_norm (const PetscScalar p) const;
+      real_type lp_norm (const real_type p) const;
 
                                        /**
                                         * Maximum absolute value of the
@@ -461,7 +461,7 @@ namespace PETScWrappers
                                        * The model for this function
                                        * should be similer to add ().
                                         */
-      VectorBase & mult (); 
+      VectorBase & mult ();
 
                                        /**
                                         * Return whether the vector contains
@@ -482,13 +482,13 @@ namespace PETScWrappers
                                         * really all positive (or zero).
                                         */
       bool is_non_negative () const;
-      
+
                                        /**
                                         * Multiply the entire vector by a
                                         * fixed factor.
                                         */
       VectorBase & operator *= (const PetscScalar factor);
-    
+
                                        /**
                                         * Divide the entire vector by a
                                         * fixed factor.
@@ -513,26 +513,26 @@ namespace PETScWrappers
                                         * scalar and not a vector.
                                         */
       void add (const PetscScalar s);
-    
+
                                        /**
                                         * Simple vector addition, equal to the
                                         * <tt>operator +=</tt>.
                                         */
       void add (const VectorBase &V);
-    
+
                                        /**
                                         * Simple addition of a multiple of a
                                         * vector, i.e. <tt>*this += a*V</tt>.
                                         */
       void add (const PetscScalar a, const VectorBase &V);
-    
+
                                        /**
                                         * Multiple addition of scaled vectors,
                                         * i.e. <tt>*this += a*V+b*W</tt>.
                                         */
       void add (const PetscScalar a, const VectorBase &V,
                 const PetscScalar b, const VectorBase &W);
-    
+
                                        /**
                                         * Scaling and simple vector addition,
                                         * i.e.
@@ -540,7 +540,7 @@ namespace PETScWrappers
                                         */
       void sadd (const PetscScalar s,
                  const VectorBase     &V);
-    
+
                                        /**
                                         * Scaling and simple addition, i.e.
                                         * <tt>*this = s*(*this)+a*V</tt>.
@@ -548,7 +548,7 @@ namespace PETScWrappers
       void sadd (const PetscScalar s,
                  const PetscScalar a,
                  const VectorBase     &V);
-    
+
                                        /**
                                         * Scaling and multiple addition.
                                         */
@@ -557,7 +557,7 @@ namespace PETScWrappers
                  const VectorBase     &V,
                  const PetscScalar b,
                  const VectorBase     &W);
-    
+
                                        /**
                                         * Scaling and multiple addition.
                                         * <tt>*this = s*(*this)+a*V + b*W + c*X</tt>.
@@ -566,10 +566,10 @@ namespace PETScWrappers
                  const PetscScalar a,
                  const VectorBase     &V,
                  const PetscScalar b,
-                 const VectorBase     &W, 
+                 const VectorBase     &W,
                  const PetscScalar c,
                  const VectorBase     &X);
-    
+
                                        /**
                                         * Scale each element of this
                                         * vector by the corresponding
@@ -580,12 +580,12 @@ namespace PETScWrappers
                                         * diagonal scaling matrix.
                                         */
       void scale (const VectorBase &scaling_factors);
-    
+
                                        /**
                                         * Assignment <tt>*this = a*V</tt>.
                                         */
       void equ (const PetscScalar a, const VectorBase &V);
-    
+
                                        /**
                                         * Assignment <tt>*this = a*V + b*W</tt>.
                                         */
@@ -654,7 +654,7 @@ namespace PETScWrappers
                                         * to standard functions.
                                         */
       void swap (VectorBase &v);
-      
+
                                        /**
                                         * Conversion operator to gain access
                                         * to the underlying PETSc type. If you
@@ -674,7 +674,7 @@ namespace PETScWrappers
                                        * for this class).
                                        */
       unsigned int memory_consumption () const;
-      
+
     protected:
                                        /**
                                         * A generic vector object in
@@ -741,7 +741,7 @@ namespace PETScWrappers
 
 
 
-// ------------------- inline and template functions --------------  
+// ------------------- inline and template functions --------------
 
 /**
  * Global function @p swap which overloads the default implementation
@@ -784,7 +784,7 @@ namespace PETScWrappers
     }
 
 
-    
+
     inline
     const VectorReference &
     VectorReference::operator = (const PetscScalar &value) const
@@ -798,15 +798,15 @@ namespace PETScWrappers
           ierr = VecAssemblyEnd (vector);
           AssertThrow (ierr == 0, ExcPETScError(ierr));
         }
-      
+
       const signed int petsc_i = index;
-      
+
       const int ierr
         = VecSetValues (vector, 1, &petsc_i, &value, INSERT_VALUES);
       AssertThrow (ierr == 0, ExcPETScError(ierr));
 
       vector.last_action = VectorBase::LastAction::insert;
-      
+
       return *this;
     }
 
@@ -837,14 +837,14 @@ namespace PETScWrappers
                                        // zero
       if (value == PetscScalar())
         return *this;
-      
+
                                        // use the PETSc function to add something
       const signed int petsc_i = index;
       const int ierr
         = VecSetValues (vector, 1, &petsc_i, &value, ADD_VALUES);
       AssertThrow (ierr == 0, ExcPETScError(ierr));
 
-      
+
       return *this;
     }
 
@@ -865,7 +865,7 @@ namespace PETScWrappers
 
           vector.last_action = VectorBase::LastAction::add;
         }
-      
+
                                        // we have to do above actions in any
                                        // case to be consistent with the MPI
                                        // communication model (see the
@@ -882,7 +882,7 @@ namespace PETScWrappers
       const int ierr
         = VecSetValues (vector, 1, &petsc_i, &subtractand, ADD_VALUES);
       AssertThrow (ierr == 0, ExcPETScError(ierr));
-      
+
       return *this;
     }
 
@@ -903,7 +903,7 @@ namespace PETScWrappers
 
           vector.last_action = VectorBase::LastAction::insert;
         }
-      
+
                                        // we have to do above actions in any
                                        // case to be consistent with the MPI
                                        // communication model (see the
@@ -918,11 +918,11 @@ namespace PETScWrappers
 
       const PetscScalar new_value
         = static_cast<PetscScalar>(*this) * value;
-      
+
       const int ierr
         = VecSetValues (vector, 1, &petsc_i, &new_value, INSERT_VALUES);
       AssertThrow (ierr == 0, ExcPETScError(ierr));
-      
+
       return *this;
     }
 
@@ -943,7 +943,7 @@ namespace PETScWrappers
 
           vector.last_action = VectorBase::LastAction::insert;
         }
-      
+
                                        // we have to do above actions in any
                                        // case to be consistent with the MPI
                                        // communication model (see the
@@ -958,15 +958,15 @@ namespace PETScWrappers
 
       const PetscScalar new_value
         = static_cast<PetscScalar>(*this) / value;
-      
+
       const int ierr
         = VecSetValues (vector, 1, &petsc_i, &new_value, INSERT_VALUES);
       AssertThrow (ierr == 0, ExcPETScError(ierr));
-      
+
       return *this;
     }
   }
-  
+
 
 
   inline
@@ -977,7 +977,7 @@ namespace PETScWrappers
     const int ierr = VecGetOwnershipRange (static_cast<const Vec &>(vector),
                                           &begin, &end);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
-    
+
     return ((index >= static_cast<unsigned int>(begin)) &&
             (index < static_cast<unsigned int>(end)));
   }
@@ -999,7 +999,7 @@ namespace PETScWrappers
   {
     return static_cast<PetscScalar>(internal::VectorReference (*this, index));
   }
-  
+
 
 #endif // DOXYGEN
 }
index c1d9f089ee5c0c6a2e11c0544d29f7e560ced8d9..9e817471eece4f37010f7c1c7a44b167eaaefd81 100644 (file)
@@ -49,7 +49,7 @@ namespace internal
   {
     Assert (false,
            ExcMessage ("Complex numbers do not have an ordering."));
-    
+
     return false;
   }
 
@@ -63,7 +63,7 @@ namespace internal
     else
       std::printf (" %5.2f", double(t));
   }
-  
+
 
 
   template <typename T>
@@ -75,10 +75,41 @@ namespace internal
     else
       std::printf (" %5.2f+%5.2fi",
                   double(t.real()), double(t.imag()));
+  }
+
+                                  // call std::copy, except for in
+                                  // the case where we want to copy
+                                  // from std::complex to a
+                                  // non-complex type
+  template <typename T, typename U>
+  void copy (const T *begin,
+            const T *end,
+            U       *dest)
+  {
+    std::copy (begin, end, dest);
+  }
+
+  template <typename T, typename U>
+  void copy (const std::complex<T> *begin,
+            const std::complex<T> *end,
+            std::complex<U>       *dest)
+  {
+    std::copy (begin, end, dest);
+  }
+
+  template <typename T, typename U>
+  void copy (const std::complex<T> *,
+            const std::complex<T> *,
+            U                     *)
+  {
+    Assert (false, ExcMessage ("Can't convert a vector of complex numbers "
+                              "into a vector of reals/doubles"));
   }
 }
 
 
+
+
 template <typename Number>
 Vector<Number>::Vector (const Vector<Number>& v)
                 :
@@ -119,6 +150,7 @@ Vector<Number>::Vector (const Vector<OtherNumber>& v)
 
 #ifdef DEAL_II_USE_PETSC
 
+
 template <typename Number>
 Vector<Number>::Vector (const PETScWrappers::Vector &v)
                 :
@@ -137,8 +169,8 @@ Vector<Number>::Vector (const PETScWrappers::Vector &v)
       PetscScalar *start_ptr;
       int ierr = VecGetArray (static_cast<const Vec&>(v), &start_ptr);
       AssertThrow (ierr == 0, ExcPETScError(ierr));
-      
-      std::copy (start_ptr, start_ptr+vec_size, begin());
+
+      internal::copy (start_ptr, start_ptr+vec_size, begin());
 
                                        // restore the representation of the
                                        // vector
@@ -198,7 +230,7 @@ Vector<Number>::Vector (const TrilinosWrappers::MPI::Vector &v)
 
       int ierr = localized_vector.trilinos_vector().ExtractView (&start_ptr);
       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-      
+
       std::copy (start_ptr[0], start_ptr[0]+vec_size, begin());
     }
 }
@@ -224,7 +256,7 @@ Vector<Number>::Vector (const TrilinosWrappers::Vector &v)
 
       int ierr = v.trilinos_vector().ExtractView (&start_ptr);
       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-      
+
       std::copy (start_ptr[0], start_ptr[0]+vec_size, begin());
     }
 }
@@ -272,7 +304,7 @@ bool
 Vector<Number>::is_non_negative () const
 {
   Assert (vec_size!=0, ExcEmptyObject());
-  
+
   for (unsigned int i=0; i<vec_size; ++i)
     if ( ! internal::is_non_negative (val[i]))
       return false;
@@ -287,13 +319,13 @@ template <typename Number2>
 Number Vector<Number>::operator * (const Vector<Number2>& v) const
 {
   Assert (vec_size!=0, ExcEmptyObject());
-  
+
   if (PointerComparison::equal (this, &v))
     return norm_sqr();
-  
+
   Assert (vec_size == v.size(),
          ExcDimensionMismatch(vec_size, v.size()));
-  
+
   Number sum = 0;
 
                                   // multiply the two vectors. we have to
@@ -304,7 +336,7 @@ Number Vector<Number>::operator * (const Vector<Number2>& v) const
                                   // is not defined by default
   for (unsigned int i=0; i<vec_size; ++i)
     sum += val[i] * Number(numbers::NumberTraits<Number2>::conjugate(v.val[i]));
-    
+
   return sum;
 }
 
@@ -319,7 +351,7 @@ Vector<Number>::norm_sqr () const
 
   for (unsigned int i=0; i<vec_size; ++i)
     sum += numbers::NumberTraits<Number>::abs_square(val[i]);
-  
+
   return sum;
 }
 
@@ -333,7 +365,7 @@ Number Vector<Number>::mean_value () const
 
   for (unsigned int i=0; i<vec_size; ++i)
     sum += val[i];
-  
+
   return sum / real_type(size());
 }
 
@@ -374,7 +406,7 @@ Vector<Number>::lp_norm (const real_type p) const
 
   for (unsigned int i=0; i<vec_size; ++i)
     sum += std::pow (numbers::NumberTraits<Number>::abs(val[i]), p);
-  
+
   return std::pow(sum, static_cast<real_type>(1./p));
 }
 
@@ -454,9 +486,9 @@ template <typename Number>
 void Vector<Number>::add (const Number a, const Vector<Number>& v,
                          const Number b, const Vector<Number>& w)
 {
-  Assert (numbers::is_finite(a), 
+  Assert (numbers::is_finite(a),
           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
-  Assert (numbers::is_finite(b), 
+  Assert (numbers::is_finite(b),
           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
   Assert (vec_size!=0, ExcEmptyObject());
@@ -477,7 +509,7 @@ template <typename Number>
 void Vector<Number>::sadd (const Number x,
                           const Vector<Number>& v)
 {
-  Assert (numbers::is_finite(x), 
+  Assert (numbers::is_finite(x),
           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
   Assert (vec_size!=0, ExcEmptyObject());
@@ -498,11 +530,11 @@ void Vector<Number>::sadd (const Number x, const Number a,
                           const Vector<Number>& v, const Number b,
                            const Vector<Number>& w)
 {
-  Assert (numbers::is_finite(x), 
+  Assert (numbers::is_finite(x),
           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
-  Assert (numbers::is_finite(a), 
+  Assert (numbers::is_finite(a),
           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
-  Assert (numbers::is_finite(b), 
+  Assert (numbers::is_finite(b),
           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
   Assert (vec_size!=0, ExcEmptyObject());
@@ -564,7 +596,7 @@ template <typename Number>
 void Vector<Number>::equ (const Number a,
                          const Vector<Number>& u)
 {
-  Assert (numbers::is_finite(a), 
+  Assert (numbers::is_finite(a),
          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
   Assert (vec_size!=0, ExcEmptyObject());
@@ -584,7 +616,7 @@ template <typename Number2>
 void Vector<Number>::equ (const Number a,
                          const Vector<Number2>& u)
 {
-  Assert (numbers::is_finite(a), 
+  Assert (numbers::is_finite(a),
          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
   Assert (vec_size!=0, ExcEmptyObject());
@@ -606,9 +638,9 @@ template <typename Number>
 void Vector<Number>::equ (const Number a, const Vector<Number>& u,
                          const Number b, const Vector<Number>& v)
 {
-  Assert (numbers::is_finite(a), 
+  Assert (numbers::is_finite(a),
           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
-  Assert (numbers::is_finite(b), 
+  Assert (numbers::is_finite(b),
           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
   Assert (vec_size!=0, ExcEmptyObject());
@@ -677,7 +709,7 @@ Vector<Number>::operator = (const BlockVector<Number>& v)
   for (unsigned int b=0; b<v.n_blocks(); ++b)
     for (unsigned int i=0; i<v.block(b).size(); ++i, ++this_index)
       val[this_index] = v.block(b)(i);
-  
+
   return *this;
 }
 
@@ -698,8 +730,8 @@ Vector<Number>::operator = (const PETScWrappers::Vector &v)
       PetscScalar *start_ptr;
       int ierr = VecGetArray (static_cast<const Vec&>(v), &start_ptr);
       AssertThrow (ierr == 0, ExcPETScError(ierr));
-      
-      std::copy (start_ptr, start_ptr+vec_size, begin());
+
+      internal::copy (start_ptr, start_ptr+vec_size, begin());
 
                                        // restore the representation of the
                                        // vector
@@ -758,7 +790,7 @@ Vector<Number>::operator = (const TrilinosWrappers::Vector &v)
       TrilinosScalar **start_ptr;
       int ierr = v.trilinos_vector().ExtractView (&start_ptr);
       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-      
+
       std::copy (start_ptr[0], start_ptr[0]+vec_size, begin());
     }
 
@@ -813,7 +845,7 @@ void Vector<Number>::print (std::ostream      &out,
 
   std::ios::fmtflags old_flags = out.flags();
   unsigned int old_precision = out.precision (precision);
-  
+
   out.precision (precision);
   if (scientific)
     out.setf (std::ios::scientific, std::ios::floatfield);
@@ -827,7 +859,7 @@ void Vector<Number>::print (std::ostream      &out,
     for (unsigned int i=0; i<size(); ++i)
       out << val[i] << std::endl;
   out << std::endl;
-  
+
   AssertThrow (out, ExcIO());
                                    // reset output format
   out.flags (old_flags);
@@ -849,19 +881,19 @@ void Vector<Number>::block_write (std::ostream &out) const
                                   // environment
   const unsigned int sz = size();
   char buf[16];
-  
+
   std::sprintf(buf, "%d", sz);
   std::strcat(buf, "\n[");
-  
+
   out.write(buf, std::strlen(buf));
   out.write (reinterpret_cast<const char*>(begin()),
             reinterpret_cast<const char*>(end())
             - reinterpret_cast<const char*>(begin()));
-  
+
                                   // out << ']';
   const char outro = ']';
   out.write (&outro, 1);
-  
+
   AssertThrow (out, ExcIO());
 }
 
@@ -875,24 +907,24 @@ void Vector<Number>::block_read (std::istream &in)
   unsigned int sz;
 
   char buf[16];
-  
+
 
   in.getline(buf,16,'\n');
   sz=std::atoi(buf);
-  
+
                                   // fast initialization, since the
                                   // data elements are overwritten anyway
-  reinit (sz, true);     
+  reinit (sz, true);
 
   char c;
                                   //  in >> c;
   in.read (&c, 1);
   AssertThrow (c=='[', ExcIO());
-  
+
   in.read (reinterpret_cast<char*>(begin()),
           reinterpret_cast<const char*>(end())
           - reinterpret_cast<const char*>(begin()));
-  
+
                                   //  in >> c;
   in.read (&c, 1);
   AssertThrow (c==']', ExcIO());
index b6d27f86481430955f9441724210437dc9540713..2508eaf432f849e3161e646437dc246f8af1d3c3 100644 (file)
@@ -41,7 +41,7 @@ namespace PETScWrappers
 
           return;
         }
-      
+
                                        // otherwise first flush PETSc caches
       matrix->compress ();
 
@@ -58,7 +58,7 @@ namespace PETScWrappers
       const int         *colnums;
       const PetscScalar *values;
 #endif
-      
+
       int ierr;
       ierr = MatGetRow(*matrix, this->a_row, &ncols, &colnums, &values);
       AssertThrow (ierr == 0, MatrixBase::ExcPETScError(ierr));
@@ -73,7 +73,7 @@ namespace PETScWrappers
       colnum_cache.reset (new std::vector<unsigned int> (colnums,
                                                          colnums+ncols));
       value_cache.reset (new std::vector<PetscScalar> (values, values+ncols));
-      
+
                                        // and finally restore the matrix
       ierr = MatRestoreRow(*matrix, this->a_row, &ncols, &colnums, &values);
       AssertThrow (ierr == 0, MatrixBase::ExcPETScError(ierr));
@@ -81,18 +81,18 @@ namespace PETScWrappers
   }
 
 
-  
+
   MatrixBase::MatrixBase ()
                   :
                   last_action (LastAction::none)
   {}
 
-  
+
 
   MatrixBase::~MatrixBase ()
   {
     const int ierr = MatDestroy (matrix);
-    AssertThrow (ierr == 0, ExcPETScError(ierr));    
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
 
 
@@ -102,7 +102,7 @@ namespace PETScWrappers
   {
                                      // destroy the matrix...
     int ierr = MatDestroy (matrix);
-    AssertThrow (ierr == 0, ExcPETScError(ierr));    
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
                                      // ...and replace it by an empty
                                      // sequential matrix
     const int m=0, n=0, n_nonzero_per_row=0;
@@ -110,20 +110,20 @@ namespace PETScWrappers
                            0, &matrix);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
-  
+
 
 
   MatrixBase &
-  MatrixBase::operator = (const double d)
+  MatrixBase::operator = (const value_type d)
   {
-    Assert (d==0, ExcScalarAssignmentOnlyForZeroValue());    
+    Assert (d==value_type(), ExcScalarAssignmentOnlyForZeroValue());
 
                                      // flush previously cached elements. this
                                      // seems to be necessary since petsc
                                      // 2.2.1, at least for parallel vectors
                                      // (see test bits/petsc_64)
     compress ();
-    
+
     const int ierr = MatZeroEntries (matrix);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
@@ -147,10 +147,10 @@ namespace PETScWrappers
 #else
     const PetscInt petsc_row = row;
 #endif
-    
+
     IS index_set;
     ISCreateGeneral (get_mpi_communicator(), 1, &petsc_row, &index_set);
-    
+
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR <= 2)
     const int ierr
       = MatZeroRows(matrix, index_set, &new_diag_value);
@@ -158,11 +158,11 @@ namespace PETScWrappers
     const int ierr
       = MatZeroRowsIS(matrix, index_set, new_diag_value);
 #endif
-    
+
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
     ISDestroy (index_set);
-    
+
     compress ();
   }
 
@@ -190,7 +190,7 @@ namespace PETScWrappers
     IS index_set;
     ISCreateGeneral (get_mpi_communicator(), rows.size(),
                      &petsc_rows[0], &index_set);
-    
+
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR <= 2)
     const int ierr
       = MatZeroRows(matrix, index_set, &new_diag_value);
@@ -198,15 +198,15 @@ namespace PETScWrappers
     const int ierr
       = MatZeroRowsIS(matrix, index_set, new_diag_value);
 #endif
-    
+
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
     ISDestroy (index_set);
-    
+
     compress ();
   }
-  
-  
+
+
 
   PetscScalar
   MatrixBase::el (const unsigned int i,
@@ -215,7 +215,7 @@ namespace PETScWrappers
     const signed int petsc_i = i;
     const signed int petsc_j = j;
     PetscScalar value;
-    
+
     const int ierr
       = MatGetValues (matrix, 1, &petsc_i, 1, &petsc_j,
                       &value);
@@ -230,14 +230,14 @@ namespace PETScWrappers
   MatrixBase::diag_element (const unsigned int i) const
   {
     Assert (m() == n(), ExcNotQuadratic());
-    
+
                                      // this doesn't seem to work any
                                      // different than any other element
     return el(i,i);
   }
-  
 
-  
+
+
   void
   MatrixBase::compress ()
   {
@@ -265,7 +265,7 @@ namespace PETScWrappers
 
     return n_rows;
   }
-  
+
 
 
   unsigned int
@@ -291,7 +291,7 @@ namespace PETScWrappers
   }
 
 
-  
+
   std::pair<unsigned int, unsigned int>
   MatrixBase::local_range () const
   {
@@ -317,7 +317,7 @@ namespace PETScWrappers
   }
 
 
-  
+
   unsigned int
   MatrixBase::
   row_length (const unsigned int row) const
@@ -364,21 +364,21 @@ namespace PETScWrappers
   MatrixBase::l1_norm () const
   {
     PetscReal result;
-    
+
     const int ierr
       = MatNorm (matrix, NORM_1, &result);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
     return result;
   }
-  
-  
+
+
 
   PetscReal
   MatrixBase::linfty_norm () const
   {
     PetscReal result;
-    
+
     const int ierr
       = MatNorm (matrix, NORM_INFINITY, &result);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
@@ -392,7 +392,7 @@ namespace PETScWrappers
   MatrixBase::frobenius_norm () const
   {
     PetscReal result;
-    
+
     const int ierr
       = MatNorm (matrix, NORM_FROBENIUS, &result);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
@@ -421,7 +421,7 @@ namespace PETScWrappers
   MatrixBase::operator /= (const PetscScalar a)
   {
     const PetscScalar factor = 1./a;
-    
+
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
     const int ierr = MatScale (&factor, matrix);
 #else
@@ -439,13 +439,13 @@ namespace PETScWrappers
                      const VectorBase &src) const
   {
     Assert (&src != &dst, ExcSourceEqualsDestination());
-    
+
     const int ierr = MatMult (matrix, src, dst);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
-  
 
-  
+
+
   void
   MatrixBase::Tvmult (VectorBase       &dst,
                       const VectorBase &src) const
@@ -455,15 +455,15 @@ namespace PETScWrappers
     const int ierr = MatMultTranspose (matrix, src, dst);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
-  
 
-  
+
+
   void
   MatrixBase::vmult_add (VectorBase       &dst,
                          const VectorBase &src) const
   {
     Assert (&src != &dst, ExcSourceEqualsDestination());
-    
+
     const int ierr = MatMultAdd (matrix, src, dst, dst);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
@@ -475,7 +475,7 @@ namespace PETScWrappers
                           const VectorBase &src) const
   {
     Assert (&src != &dst, ExcSourceEqualsDestination());
-    
+
     const int ierr = MatMultTransposeAdd (matrix, src, dst, dst);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
@@ -489,9 +489,9 @@ namespace PETScWrappers
     vmult (tmp, v);
     return tmp*v;
   }
-  
 
-  
+
+
   PetscScalar
   MatrixBase::matrix_scalar_product (const VectorBase &u,
                                      const VectorBase &v) const
@@ -500,9 +500,9 @@ namespace PETScWrappers
     vmult (tmp, v);
     return u*tmp;
   }
-  
 
-  
+
+
   PetscScalar
   MatrixBase::residual (VectorBase       &dst,
                         const VectorBase &x,
@@ -514,23 +514,23 @@ namespace PETScWrappers
     vmult (dst, x);
     dst -= b;
     dst *= -1;
-    
+
     return dst.l2_norm();
   }
-  
 
-  
+
+
   MatrixBase::operator const Mat () const
   {
     return matrix;
-  }  
+  }
 
   void
   MatrixBase::transpose ()
   {
     int ierr;
 
-#if (PETSC_VERSION_MAJOR <= 2) 
+#if (PETSC_VERSION_MAJOR <= 2)
     ierr = MatTranspose(matrix, PETSC_NULL);
 #else
     ierr = MatTranspose(matrix, MAT_REUSE_MATRIX, &matrix);
@@ -540,17 +540,17 @@ namespace PETScWrappers
   }
 
   PetscTruth
-  MatrixBase::is_symmetric (const double tolerance) 
+  MatrixBase::is_symmetric (const double tolerance)
   {
     PetscTruth truth;
                                        // First flush PETSc caches
     compress ();
     MatIsSymmetric (matrix, tolerance, &truth);
     return truth;
-  }  
+  }
 
   PetscTruth
-  MatrixBase::is_hermitian (const double tolerance) 
+  MatrixBase::is_hermitian (const double tolerance)
   {
     PetscTruth truth;
                                     // First flush PETSc caches
@@ -559,13 +559,13 @@ namespace PETScWrappers
 #if (PETSC_VERSION_MAJOR <= 2)
                                     // avoid warning about unused variables
     (void) tolerance;
-    
+
     MatIsHermitian (matrix, &truth);
 #else
     MatIsHermitian (matrix, tolerance, &truth);
 #endif
     return truth;
-  }  
+  }
 
   void
   MatrixBase::write_ascii ()
@@ -574,7 +574,7 @@ namespace PETScWrappers
     compress ();
 
                                        // Set options
-#if (PETSC_VERSION_MAJOR <= 2) 
+#if (PETSC_VERSION_MAJOR <= 2)
     PetscViewerSetFormat (PETSC_VIEWER_STDOUT_WORLD,
                          PETSC_VIEWER_ASCII_DEFAULT);
 #else
index 0b9c1c2a76360f426c9e8681003fa327cd8d7a00..0d804eb214602d049a459578fafcbd964d384e19 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004, 2005, 2006, 2008 by the deal.II authors
+//    Copyright (C) 2004, 2005, 2006, 2008, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -27,7 +27,7 @@ namespace PETScWrappers
 {
   namespace MPI
   {
-    
+
     SparseMatrix::SparseMatrix ()
     {
                                        // just like for vectors: since we
@@ -89,17 +89,17 @@ namespace PETScWrappers
                  local_columns_per_process, this_process,
                  preset_nonzero_locations);
     }
-    
+
 
 
     SparseMatrix &
-    SparseMatrix::operator = (const double d)
+    SparseMatrix::operator = (const value_type d)
     {
       MatrixBase::operator = (d);
       return *this;
     }
 
-    
+
 
     void
     SparseMatrix::reinit (const MPI_Comm    &communicator,
@@ -115,8 +115,8 @@ namespace PETScWrappers
                                        // get rid of old matrix and generate a
                                        // new one
       const int ierr = MatDestroy (matrix);
-      AssertThrow (ierr == 0, ExcPETScError(ierr));    
-    
+      AssertThrow (ierr == 0, ExcPETScError(ierr));
+
       do_reinit (m, n, local_rows, local_columns,
                  n_nonzero_per_row, is_symmetric);
     }
@@ -133,14 +133,14 @@ namespace PETScWrappers
                           const bool         is_symmetric)
     {
       this->communicator = communicator;
-      
+
                                        // get rid of old matrix and generate a
                                        // new one
       const int ierr = MatDestroy (matrix);
-      AssertThrow (ierr == 0, ExcPETScError(ierr));    
+      AssertThrow (ierr == 0, ExcPETScError(ierr));
 
       do_reinit (m, n, local_rows, local_columns, row_lengths, is_symmetric);
-    }  
+    }
 
 
 
@@ -155,18 +155,18 @@ namespace PETScWrappers
             const bool                       preset_nonzero_locations)
     {
       this->communicator = communicator;
-      
+
                                        // get rid of old matrix and generate a
                                        // new one
       const int ierr = MatDestroy (matrix);
-      AssertThrow (ierr == 0, ExcPETScError(ierr));    
+      AssertThrow (ierr == 0, ExcPETScError(ierr));
 
       do_reinit (sparsity_pattern, local_rows_per_process,
                  local_columns_per_process, this_process,
                  preset_nonzero_locations);
     }
 
-    
+
 
     void
     SparseMatrix::do_reinit (const unsigned int m,
@@ -177,7 +177,7 @@ namespace PETScWrappers
                              const bool         is_symmetric)
     {
       Assert (local_rows <= m, ExcLocalRowsTooLarge (local_rows, m));
-      
+
                                        // use the call sequence indicating only
                                        // a maximal number of elements per row
                                        // for all rows globally
@@ -192,7 +192,7 @@ namespace PETScWrappers
                                        // set symmetric flag, if so requested
       if (is_symmetric == true)
         {
-#if (PETSC_VERSION_MAJOR <= 2) 
+#if (PETSC_VERSION_MAJOR <= 2)
           const int ierr
             = MatSetOption (matrix, MAT_SYMMETRIC);
 #else
@@ -228,7 +228,7 @@ namespace PETScWrappers
       for (unsigned int i=0; i<row_lengths.size(); ++i)
        Assert(row_lengths[i]<=local_columns,
               ExcIndexRange(row_lengths[i], 1, local_columns+1));
-    
+
                                        // use the call sequence indicating a
                                        // maximal number of elements for each
                                        // row individually. annoyingly, we
@@ -252,7 +252,7 @@ namespace PETScWrappers
                                        // set symmetric flag, if so requested
       if (is_symmetric == true)
         {
-#if (PETSC_VERSION_MAJOR <= 2) 
+#if (PETSC_VERSION_MAJOR <= 2)
           const int ierr
             = MatSetOption (matrix, MAT_SYMMETRIC);
 #else
@@ -261,7 +261,7 @@ namespace PETScWrappers
 #endif
 
           AssertThrow (ierr == 0, ExcPETScError(ierr));
-        }    
+        }
     }
 
 
@@ -307,7 +307,7 @@ namespace PETScWrappers
         for (unsigned int c=0; c<sparsity_pattern.row_length(row); ++c)
           {
             const unsigned int column = sparsity_pattern.column_number(row,c);
-            
+
             if ((column >= local_col_start) &&
                 (column < local_col_end))
               ++row_lengths_in_window[row-local_row_start];
@@ -359,7 +359,7 @@ namespace PETScWrappers
                                        // the whole thing:
 #if (PETSC_VERSION_MAJOR <= 2) && \
     ((PETSC_VERSION_MINOR < 2) ||  \
-     ((PETSC_VERSION_MINOR == 2) && (PETSC_VERSION_SUBMINOR == 0)))      
+     ((PETSC_VERSION_MINOR == 2) && (PETSC_VERSION_SUBMINOR == 0)))
 
           std::vector<int> row_entries;
           std::vector<PetscScalar> row_values;
@@ -369,7 +369,7 @@ namespace PETScWrappers
               row_values.resize (sparsity_pattern.row_length(i), 0.0);
               for (unsigned int j=0; j<sparsity_pattern.row_length(i); ++j)
                 row_entries[j] = sparsity_pattern.column_number (i,j);
-              
+
               const int int_row = i;
               MatSetValues (matrix, 1, &int_row,
                             sparsity_pattern.row_length(i), &row_entries[0],
@@ -377,7 +377,7 @@ namespace PETScWrappers
             }
 
           compress ();
-          
+
 #else
 
                                            // first set up the column number
@@ -401,7 +401,7 @@ namespace PETScWrappers
               }
             colnums_in_window.resize (n_cols+1, -1);
           }
-          
+
                                            // now copy over the information
                                            // from the sparsity pattern.
           {
@@ -412,7 +412,7 @@ namespace PETScWrappers
               colnums_in_window[index] = sparsity_pattern.column_number(i,j);
             Assert (index == colnums_in_window.size()-1, ExcInternalError());
           }
-          
+
                                            // then call the petsc function
                                            // that summarily allocates these
                                            // entries:
@@ -441,7 +441,7 @@ namespace PETScWrappers
             const std::vector<PetscScalar>
               values (sparsity_pattern.max_entries_per_row(),
                       1.);
-            
+
             for (unsigned int i=local_row_start; i<local_row_end; ++i)
               {
                 const int petsc_i = i;
@@ -451,28 +451,28 @@ namespace PETScWrappers
                               &values[0], INSERT_VALUES);
               }
           }
-          
+
           compress ();
 
                                            // set the dummy entries set above
                                            // back to zero
           *this = 0;
           compress ();
-          
+
 #endif
 
                                           // Now we won't insert any
                                           // further entries, so PETSc can
                                           // internally optimize some data
                                           // structures.
-#if (PETSC_VERSION_MAJOR <= 2) 
+#if (PETSC_VERSION_MAJOR <= 2)
           const int ierr =
            MatSetOption (matrix, MAT_NO_NEW_NONZERO_LOCATIONS);
 #else
           const int ierr =
            MatSetOption (matrix, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE);
 #endif
-         
+
          AssertThrow (ierr == 0, ExcPETScError(ierr));
         }
     }
index 907837d42a29bcdcfef3b457d612be3b02caf764..04cc47144ba337d30af57cdf734207be522d52dd 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004, 2006, 2008 by the deal.II authors
+//    Copyright (C) 2004, 2006, 2008, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -35,8 +35,8 @@ namespace PETScWrappers
   {
                                      // destroy the solver object
     int ierr = KSPDestroy (ksp);
-    AssertThrow (ierr == 0, ExcPETScError(ierr)); 
-  
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
                                      // and destroy the solver object if we
                                      // are in old PETSc mode
 #if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR < 2)
@@ -45,8 +45,8 @@ namespace PETScWrappers
 #endif
   }
 
-  
-  
+
+
   SolverBase::SolverBase (SolverControl  &cn,
                           const MPI_Comm &mpi_communicator)
                   :
@@ -54,12 +54,12 @@ namespace PETScWrappers
                   mpi_communicator (mpi_communicator)
   {}
 
-  
+
 
   SolverBase::~SolverBase ()
   {}
 
-  
+
 
   void
   SolverBase::solve (const MatrixBase         &A,
@@ -80,7 +80,7 @@ namespace PETScWrappers
     SLES sles;
     ierr = SLESCreate (mpi_communicator, &sles);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
-    
+
                                      // set the matrices involved. the
                                      // last argument is irrelevant here,
                                      // since we use the solver only once
@@ -88,29 +88,29 @@ namespace PETScWrappers
     ierr = SLESSetOperators (sles, A, preconditioner,
                              SAME_PRECONDITIONER);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
-    
+
                                      // let derived classes set the solver
                                      // type, and the preconditioning object
                                      // set the type of preconditioner
     KSP ksp;
     ierr = SLESGetKSP (sles, &ksp);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
-    
+
     set_solver_type (ksp);
-    
+
     PC pc;
     ierr = SLESGetPC (sles, &pc);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
-    
+
     preconditioner.set_preconditioner_type (pc);
-  
+
                                      // then a convergence monitor
                                      // function. that function simply checks
                                      // with the solver_control object we have
                                      // in this object for convergence
     KSPSetConvergenceTest (ksp, &convergence_test,
                            reinterpret_cast<void *>(&solver_control));
-      
+
                                      // then do the real work: set up solver
                                      // internal data and solve the
                                      // system. this could be joined in one
@@ -120,7 +120,7 @@ namespace PETScWrappers
     int iterations = 0;
     ierr = SLESSetUp (sles, b, x);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
-      
+
     ierr = SLESSolve (sles, b, x, &iterations);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
@@ -135,7 +135,7 @@ namespace PETScWrappers
     if (solver_data == 0)
       {
         solver_data.reset (new SolverData());
-        
+
         ierr = KSPCreate (mpi_communicator, &solver_data->ksp);
         AssertThrow (ierr == 0, ExcPETScError(ierr));
 
@@ -152,7 +152,7 @@ namespace PETScWrappers
                                          // object set the type of
                                          // preconditioner
         set_solver_type (solver_data->ksp);
-    
+
         ierr = KSPGetPC (solver_data->ksp, &solver_data->pc);
         AssertThrow (ierr == 0, ExcPETScError(ierr));
 
@@ -163,7 +163,7 @@ namespace PETScWrappers
                                          // checks with the solver_control
                                          // object we have in this object for
                                          // convergence
-#if (PETSC_VERSION_MAJOR <= 2) 
+#if (PETSC_VERSION_MAJOR <= 2)
         KSPSetConvergenceTest (solver_data->ksp, &convergence_test,
                                reinterpret_cast<void *>(&solver_control));
 #else
@@ -173,7 +173,7 @@ namespace PETScWrappers
 #endif
 
       }
-    
+
                                      // then do the real work: set up solver
                                      // internal data and solve the
                                      // system. unfortunately, the call
@@ -186,7 +186,7 @@ namespace PETScWrappers
     AssertThrow (ierr == 0, ExcPETScError(ierr));
     ierr = KSPSetUp (solver_data->ksp);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
-    
+
     ierr = KSPSolve (solver_data->ksp);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 #  else
@@ -195,7 +195,7 @@ namespace PETScWrappers
 
     ierr = KSPSolve (solver_data->ksp, b, x);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
-#  endif 
+#  endif
 
 #endif
                                      // destroy solver object
@@ -216,38 +216,38 @@ namespace PETScWrappers
   {
     return solver_control;
   }
-  
+
 
 
   int
   SolverBase::convergence_test (KSP                 /*ksp*/,
                                 const int           iteration,
-                                const PetscScalar   residual_norm,
+                                const PetscReal     residual_norm,
                                 KSPConvergedReason *reason,
                                 void               *solver_control_x)
   {
     SolverControl &solver_control = *reinterpret_cast<SolverControl*>(solver_control_x);
-    
+
     const SolverControl::State state
         = solver_control.check (iteration, residual_norm);
-    
+
     switch (state)
       {
         case ::dealii::SolverControl::iterate:
               *reason = KSP_CONVERGED_ITERATING;
               break;
-              
+
         case ::dealii::SolverControl::success:
               *reason = static_cast<KSPConvergedReason>(1);
               break;
-              
+
         case ::dealii::SolverControl::failure:
               if (solver_control.last_step() > solver_control.max_steps())
                 *reason = KSP_DIVERGED_ITS;
               else
                 *reason = KSP_DIVERGED_DTOL;
               break;
-              
+
         default:
               Assert (false, ExcNotImplemented());
       }
@@ -255,8 +255,8 @@ namespace PETScWrappers
                                      // return without failure
     return 0;
   }
-  
-    
+
+
 
 /* ---------------------- SolverRichardson ------------------------ */
 
@@ -266,8 +266,8 @@ namespace PETScWrappers
                   omega (omega)
   {}
 
-  
-  
+
+
   SolverRichardson::SolverRichardson (SolverControl        &cn,
                                       const MPI_Comm       &mpi_communicator,
                                       const AdditionalData &data)
@@ -313,7 +313,7 @@ namespace PETScWrappers
     KSPSetTolerances(ksp, PETSC_DEFAULT, this->solver_control.tolerance(),
                     PETSC_DEFAULT, this->solver_control.max_steps()+1);
   }
-  
+
 
 /* ---------------------- SolverChebychev ------------------------ */
 
@@ -343,7 +343,7 @@ namespace PETScWrappers
     KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
   }
 
-  
+
 /* ---------------------- SolverCG ------------------------ */
 
   SolverCG::SolverCG (SolverControl        &cn,
@@ -371,7 +371,7 @@ namespace PETScWrappers
                                      // solution vector. do so here as well:
     KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
   }
-  
+
 
 /* ---------------------- SolverBiCG ------------------------ */
 
@@ -401,7 +401,7 @@ namespace PETScWrappers
     KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
   }
 
-  
+
 /* ---------------------- SolverGMRES ------------------------ */
 
   SolverGMRES::AdditionalData::
@@ -412,8 +412,8 @@ namespace PETScWrappers
                  right_preconditioning (right_preconditioning)
   {}
 
-  
-  
+
+
   SolverGMRES::SolverGMRES (SolverControl        &cn,
                             const MPI_Comm       &mpi_communicator,
                             const AdditionalData &data)
@@ -454,7 +454,7 @@ namespace PETScWrappers
                                      // so rather expand their macros by hand,
                                      // and do some equally nasty stuff that at
                                      // least doesn't yield warnings...
-    int (*fun_ptr)(KSP,int); 
+    int (*fun_ptr)(KSP,int);
     ierr = PetscObjectQueryFunction((PetscObject)(ksp),
                                     "KSPGMRESSetRestart_C",
                                     (void (**)())&fun_ptr);
@@ -476,7 +476,7 @@ namespace PETScWrappers
                                      // solution vector. do so here as well:
     KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
   }
-  
+
 
 /* ---------------------- SolverBicgstab ------------------------ */
 
@@ -506,7 +506,7 @@ namespace PETScWrappers
     KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
   }
 
-  
+
 /* ---------------------- SolverCGS ------------------------ */
 
   SolverCGS::SolverCGS (SolverControl        &cn,
@@ -534,7 +534,7 @@ namespace PETScWrappers
                                      // solution vector. do so here as well:
     KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
   }
-  
+
 
 /* ---------------------- SolverTFQMR ------------------------ */
 
@@ -564,7 +564,7 @@ namespace PETScWrappers
     KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
   }
 
-  
+
 /* ---------------------- SolverTCQMR ------------------------ */
 
   SolverTCQMR::SolverTCQMR (SolverControl        &cn,
@@ -593,7 +593,7 @@ namespace PETScWrappers
     KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
   }
 
-  
+
 /* ---------------------- SolverCR ------------------------ */
 
   SolverCR::SolverCR (SolverControl        &cn,
@@ -622,7 +622,7 @@ namespace PETScWrappers
     KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
   }
 
-  
+
 /* ---------------------- SolverLSQR ------------------------ */
 
   SolverLSQR::SolverLSQR (SolverControl        &cn,
@@ -650,7 +650,7 @@ namespace PETScWrappers
                                      // solution vector. do so here as well:
     KSPSetInitialGuessNonzero (ksp, PETSC_TRUE);
   }
+
 
 /* ---------------------- SolverPreOnly ------------------------ */
 
@@ -676,11 +676,11 @@ namespace PETScWrappers
 
                                      // The KSPPREONLY solver of
                                      // PETSc never calls the convergence
-                                     // monitor, which leads to failure 
+                                     // monitor, which leads to failure
                                      // even when everything was ok.
-                                     // Therefore the SolverControl status 
+                                     // Therefore the SolverControl status
                                      // is set to some nice values, which
-                                     // guarantee a nice result at the end 
+                                     // guarantee a nice result at the end
                                      // of the solution process.
     solver_control.check (1, 0.0);
 
index 95436b78e8aaa4105ad56c9a43a65127d62986e4..a359b72f8715a35b988d523f313924b7ee3c10ee 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004, 2005, 2006, 2007, 2008 by the deal.II authors
+//    Copyright (C) 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -30,7 +30,7 @@ namespace PETScWrappers
     {
       Assert (index < vector.size(),
               ExcIndexRange (index, 0, vector.size()));
-              
+
                                        // this is clumsy: there is no simple
                                        // way in PETSc to read an element from
                                        // a vector, i.e. there is no function
@@ -49,12 +49,12 @@ namespace PETScWrappers
           int ierr
             = VecGetArray (static_cast<const Vec &>(vector), &ptr);
           AssertThrow (ierr == 0, ExcPETScError(ierr));
-          
+
           const PetscScalar value = *(ptr+index);
-          
+
           ierr = VecRestoreArray (static_cast<const Vec &>(vector), &ptr);
           AssertThrow (ierr == 0, ExcPETScError(ierr));
-          
+
           return value;
         }
       else if (dynamic_cast<const PETScWrappers::MPI::Vector *>(&vector) != 0)
@@ -76,12 +76,12 @@ namespace PETScWrappers
           PetscScalar *ptr;
           ierr = VecGetArray (static_cast<const Vec &>(vector), &ptr);
           AssertThrow (ierr == 0, ExcPETScError(ierr));
-          
+
           const PetscScalar value = *(ptr+index-begin);
-          
+
           ierr = VecRestoreArray (static_cast<const Vec &>(vector), &ptr);
           AssertThrow (ierr == 0, ExcPETScError(ierr));
-          
+
           return value;
         }
       else
@@ -89,16 +89,16 @@ namespace PETScWrappers
                                          // exists there?
         Assert (false, ExcInternalError());
       return -1e20;
-    }  
+    }
   }
-  
+
   VectorBase::VectorBase ()
                   :
                   last_action (LastAction::none)
   {}
-  
 
-  
+
+
   VectorBase::VectorBase (const VectorBase &v)
                   :
                   last_action (LastAction::none)
@@ -110,7 +110,7 @@ namespace PETScWrappers
     AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
 
-  
+
   VectorBase::~VectorBase ()
   {
     const int ierr = VecDestroy (vector);
@@ -123,7 +123,7 @@ namespace PETScWrappers
   VectorBase::operator = (const PetscScalar s)
   {
 
-    Assert (numbers::is_finite(s), 
+    Assert (numbers::is_finite(s),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
                                      // flush previously cached elements. this
@@ -131,13 +131,13 @@ namespace PETScWrappers
                                      // 2.2.1, at least for parallel vectors
                                      // (see test bits/petsc_65)
     compress ();
-    
+
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
     const int ierr = VecSet (&s, vector);
 #else
     const int ierr = VecSet (vector, s);
 #endif
-    
+
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
     return *this;
@@ -150,9 +150,9 @@ namespace PETScWrappers
   {
     Assert (size() == v.size(),
             ExcDimensionMismatch(size(), v.size()));
-    
+
     PetscTruth flag;
-    
+
     const int ierr = VecEqual (vector, v.vector, &flag);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
@@ -160,7 +160,7 @@ namespace PETScWrappers
   }
 
 
-  
+
   bool
   VectorBase::operator != (const VectorBase &v) const
   {
@@ -168,7 +168,7 @@ namespace PETScWrappers
             ExcDimensionMismatch(size(), v.size()));
 
     PetscTruth flag;
-    
+
     const int ierr = VecEqual (vector, v.vector, &flag);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
@@ -176,7 +176,7 @@ namespace PETScWrappers
   }
 
 
-  
+
   unsigned int
   VectorBase::size () const
   {
@@ -229,7 +229,7 @@ namespace PETScWrappers
   {
     do_set_add_operation(indices, values, true);
   }
-  
+
 
 
   PetscScalar
@@ -245,7 +245,7 @@ namespace PETScWrappers
 
     return result;
   }
-  
+
 
 
   void
@@ -263,7 +263,7 @@ namespace PETScWrappers
   VectorBase::real_type
   VectorBase::norm_sqr () const
   {
-    const PetscScalar d = l2_norm();
+    const real_type d = l2_norm();
     return d*d;
   }
 
@@ -300,48 +300,48 @@ namespace PETScWrappers
                                        // add up remaining elements
       while (ptr != start_ptr+size())
         sum0 += *ptr++;
-  
-      mean = (sum0+sum1+sum2+sum3)/size();
+
+      mean = (sum0+sum1+sum2+sum3)/static_cast<PetscScalar>(size());
     }
-    
+
                                      // restore the representation of the
                                      // vector
     ierr = VecRestoreArray (vector, &start_ptr);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
-      
+
     return mean;
   }
-  
+
 
 
   VectorBase::real_type
   VectorBase::l1_norm () const
   {
-    PetscScalar d;
+    real_type d;
 
     const int ierr = VecNorm (vector, NORM_1, &d);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
     return d;
   }
-  
+
 
 
   VectorBase::real_type
   VectorBase::l2_norm () const
   {
-    PetscScalar d;
+    real_type d;
 
     const int ierr = VecNorm (vector, NORM_2, &d);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
     return d;
   }
-  
+
 
 
   VectorBase::real_type
-  VectorBase::lp_norm (const PetscScalar p) const
+  VectorBase::lp_norm (const real_type p) const
   {
                                      // get a representation of the vector and
                                      // loop over all the elements
@@ -349,12 +349,12 @@ namespace PETScWrappers
     int ierr = VecGetArray (vector, &start_ptr);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
-    PetscScalar norm = 0;
+    real_type norm = 0;
     {
-      PetscScalar sum0 = 0,
-                  sum1 = 0,
-                  sum2 = 0,
-                  sum3 = 0;
+      real_type sum0 = 0,
+               sum1 = 0,
+               sum2 = 0,
+               sum3 = 0;
 
                                        // use modern processors better by
                                        // allowing pipelined commands to be
@@ -363,24 +363,24 @@ namespace PETScWrappers
       const PetscScalar * eptr = ptr + (size()/4)*4;
       while (ptr!=eptr)
         {
-          sum0 += std::pow(std::fabs(*ptr++), p);
-          sum1 += std::pow(std::fabs(*ptr++), p);
-          sum2 += std::pow(std::fabs(*ptr++), p);
-          sum3 += std::pow(std::fabs(*ptr++), p);
-        };
+          sum0 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
+          sum1 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
+          sum2 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
+          sum3 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
+        }
                                        // add up remaining elements
       while (ptr != start_ptr+size())
-        sum0 += std::pow(std::fabs(*ptr++), p);
-  
+        sum0 += std::pow(numbers::NumberTraits<value_type>::abs(*ptr++), p);
+
       norm = std::pow(sum0+sum1+sum2+sum3,
-                      static_cast<PetscScalar>(1./p));
+                      1./p);
     }
-    
+
                                      // restore the representation of the
                                      // vector
     ierr = VecRestoreArray (vector, &start_ptr);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
-      
+
     return norm;
   }
 
@@ -389,7 +389,7 @@ namespace PETScWrappers
   VectorBase::real_type
   VectorBase::linfty_norm () const
   {
-    PetscScalar d;
+    real_type d;
 
     const int ierr = VecNorm (vector, NORM_INFINITY, &d);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
@@ -402,7 +402,7 @@ namespace PETScWrappers
   VectorBase::real_type
   VectorBase::normalize () const
   {
-    PetscScalar d;
+    real_type d;
 
     const int ierr = VecNormalize (vector, &d);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
@@ -415,8 +415,8 @@ namespace PETScWrappers
   VectorBase::real_type
   VectorBase::max ()  const
   {
-    PetscInt    p;
-    PetscScalar d;
+    PetscInt  p;
+    real_type d;
 
     const int ierr = VecMax (vector, &p, &d);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
@@ -427,7 +427,7 @@ namespace PETScWrappers
 
 
   VectorBase &
-  VectorBase::abs () 
+  VectorBase::abs ()
   {
     const int ierr = VecAbs (vector);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
@@ -438,7 +438,7 @@ namespace PETScWrappers
 
 
   VectorBase &
-  VectorBase::conjugate () 
+  VectorBase::conjugate ()
   {
     const int ierr = VecConjugate (vector);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
@@ -449,7 +449,7 @@ namespace PETScWrappers
 
 
   VectorBase &
-  VectorBase::mult () 
+  VectorBase::mult ()
   {
     const int ierr = VecPointwiseMult (vector,vector,vector);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
@@ -473,23 +473,44 @@ namespace PETScWrappers
     bool flag = true;
     while (ptr != eptr)
       {
-        if (*ptr != 0)
+        if (*ptr != value_type())
           {
             flag = false;
             break;
           }
         ++ptr;
       }
-    
+
                                      // restore the representation of the
                                      // vector
     ierr = VecRestoreArray (vector, &start_ptr);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
-      
+
     return flag;
   }
 
 
+  namespace internal
+  {
+    template <typename T>
+    bool is_non_negative (const T &t)
+    {
+      return t >= 0;
+    }
+
+
+
+    template <typename T>
+    bool is_non_negative (const std::complex<T> &)
+    {
+      Assert (false,
+             ExcMessage ("You can't ask a complex value "
+                         "whether it is non-negative."))
+      return true;
+    }
+  }
+
+
 
   bool
   VectorBase::is_non_negative () const
@@ -505,29 +526,29 @@ namespace PETScWrappers
     bool flag = true;
     while (ptr != eptr)
       {
-        if (*ptr < 0.0)
+        if (! internal::is_non_negative(*ptr))
           {
             flag = false;
             break;
           }
         ++ptr;
       }
-    
+
                                      // restore the representation of the
                                      // vector
     ierr = VecRestoreArray (vector, &start_ptr);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
-      
+
     return flag;
   }
-  
+
 
 
   VectorBase &
   VectorBase::operator *= (const PetscScalar a)
   {
 
-    Assert (numbers::is_finite(a), 
+    Assert (numbers::is_finite(a),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
@@ -535,7 +556,7 @@ namespace PETScWrappers
 #else
     const int ierr = VecScale (vector, a);
 #endif
-    
+
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
     return *this;
@@ -547,12 +568,12 @@ namespace PETScWrappers
   VectorBase::operator /= (const PetscScalar a)
   {
 
-    Assert (numbers::is_finite(a), 
+    Assert (numbers::is_finite(a),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
     const PetscScalar factor = 1./a;
-    
-    Assert (numbers::is_finite(factor), 
+
+    Assert (numbers::is_finite(factor),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
@@ -572,12 +593,12 @@ namespace PETScWrappers
   VectorBase::operator += (const VectorBase &v)
   {
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
-    const PetscScalar one = 1.0;    
+    const PetscScalar one = 1.0;
     const int ierr = VecAXPY (&one, v, vector);
 #else
     const int ierr = VecAXPY (vector, 1, v);
 #endif
-    
+
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
     return *this;
@@ -589,12 +610,12 @@ namespace PETScWrappers
   VectorBase::operator -= (const VectorBase &v)
   {
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
-    const PetscScalar minus_one = -1.0;    
+    const PetscScalar minus_one = -1.0;
     const int ierr = VecAXPY (&minus_one, v, vector);
 #else
     const int ierr = VecAXPY (vector, -1, v);
 #endif
-    
+
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
     return *this;
@@ -606,7 +627,7 @@ namespace PETScWrappers
   VectorBase::add (const PetscScalar s)
   {
 
-    Assert (numbers::is_finite(s), 
+    Assert (numbers::is_finite(s),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
@@ -614,10 +635,10 @@ namespace PETScWrappers
 #else
     const int ierr = VecShift (vector, s);
 #endif
-    
+
     AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
-  
+
 
 
   void
@@ -625,7 +646,7 @@ namespace PETScWrappers
   {
     *this += v;
   }
-  
+
 
 
   void
@@ -633,7 +654,7 @@ namespace PETScWrappers
                    const VectorBase     &v)
   {
 
-    Assert (numbers::is_finite(a), 
+    Assert (numbers::is_finite(a),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
@@ -641,10 +662,10 @@ namespace PETScWrappers
 #else
     const int ierr = VecAXPY (vector, a, v);
 #endif
-    
+
     AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
-  
+
 
 
   void
@@ -654,14 +675,14 @@ namespace PETScWrappers
                    const VectorBase &w)
   {
 
-    Assert (numbers::is_finite(a), 
+    Assert (numbers::is_finite(a),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
-    Assert (numbers::is_finite(b), 
+    Assert (numbers::is_finite(b),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
     const PetscScalar weights[2] = {a,b};
     Vec               addends[2] = {v.vector, w.vector};
-    
+
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
     const int ierr = VecMAXPY (2, weights, vector, addends);
 #else
@@ -670,7 +691,7 @@ namespace PETScWrappers
 
     AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
-  
+
 
 
   void
@@ -678,7 +699,7 @@ namespace PETScWrappers
                     const VectorBase &v)
   {
 
-    Assert (numbers::is_finite(s), 
+    Assert (numbers::is_finite(s),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
@@ -686,10 +707,10 @@ namespace PETScWrappers
 #else
     const int ierr = VecAYPX (vector, s, v);
 #endif
-    
+
     AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
-  
+
 
 
   void
@@ -698,9 +719,9 @@ namespace PETScWrappers
                     const VectorBase     &v)
   {
 
-    Assert (numbers::is_finite(s), 
+    Assert (numbers::is_finite(s),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
-    Assert (numbers::is_finite(a), 
+    Assert (numbers::is_finite(a),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
                                      // there is nothing like a AXPAY
@@ -709,7 +730,7 @@ namespace PETScWrappers
     *this *= s;
     add (a,v);
   }
-  
+
 
 
   void
@@ -719,27 +740,27 @@ namespace PETScWrappers
                     const PetscScalar b,
                     const VectorBase     &w)
   {
-    Assert (numbers::is_finite(s), 
+
+    Assert (numbers::is_finite(s),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
-    Assert (numbers::is_finite(a), 
+    Assert (numbers::is_finite(a),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
-    Assert (numbers::is_finite(b), 
+    Assert (numbers::is_finite(b),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
                                      // there is no operation like MAXPAY, so
                                      // do it in two steps
     *this *= s;
-    
+
     const PetscScalar weights[2] = {a,b};
     Vec               addends[2] = {v.vector,w.vector};
-    
+
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
     const int ierr = VecMAXPY (2, weights, vector, addends);
 #else
     const int ierr = VecMAXPY (vector, 2, weights, addends);
 #endif
-    
+
     AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
 
@@ -754,14 +775,14 @@ namespace PETScWrappers
                     const PetscScalar c,
                     const VectorBase     &x)
   {
-    Assert (numbers::is_finite(s), 
+
+    Assert (numbers::is_finite(s),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
-    Assert (numbers::is_finite(a), 
+    Assert (numbers::is_finite(a),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
-    Assert (numbers::is_finite(b), 
+    Assert (numbers::is_finite(b),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
-    Assert (numbers::is_finite(c), 
+    Assert (numbers::is_finite(c),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
                                      // there is no operation like MAXPAY, so
@@ -770,13 +791,13 @@ namespace PETScWrappers
 
     const PetscScalar weights[3] = {a,b,c};
     Vec               addends[3] = {v.vector, w.vector, x.vector};
-    
+
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
     const int ierr = VecMAXPY (3, weights, vector, addends);
 #else
     const int ierr = VecMAXPY (vector, 3, weights, addends);
 #endif
-    
+
     AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
 
@@ -797,7 +818,7 @@ namespace PETScWrappers
                    const VectorBase &v)
   {
 
-    Assert (numbers::is_finite(a), 
+    Assert (numbers::is_finite(a),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
     Assert (size() == v.size(),
@@ -811,7 +832,7 @@ namespace PETScWrappers
 
     *this *= a;
   }
-  
+
 
 
   void
@@ -821,9 +842,9 @@ namespace PETScWrappers
                    const VectorBase &w)
   {
 
-    Assert (numbers::is_finite(a), 
+    Assert (numbers::is_finite(a),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
-    Assert (numbers::is_finite(b), 
+    Assert (numbers::is_finite(b),
            ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
 
     Assert (size() == v.size(),
@@ -877,29 +898,29 @@ namespace PETScWrappers
 
     if (across)
       for (unsigned int i=0; i<size(); ++i)
-        out << static_cast<double>(val[i]) << ' ';
+        out << val[i] << ' ';
     else
       for (unsigned int i=0; i<size(); ++i)
-        out << static_cast<double>(val[i]) << std::endl;
+        out << val[i] << std::endl;
     out << std::endl;
 
                                      // restore the representation of the
                                      // vector
     ierr = VecRestoreArray (vector, &val);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
-    
+
     AssertThrow (out, ExcIO());
   }
 
 
-  
+
   void
   VectorBase::swap (VectorBase &v)
   {
     const int ierr = VecSwap (vector, v.vector);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
-  
+
 
 
   VectorBase::operator const Vec & () const
@@ -915,7 +936,7 @@ namespace PETScWrappers
     AssertThrow(false, ExcNotImplemented() );
     return 0;
   }
-  
+
 
 
   void
@@ -966,7 +987,7 @@ namespace PETScWrappers
 
     last_action = LastAction::insert;
   }
-  
+
 }
 
 DEAL_II_NAMESPACE_CLOSE

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.