]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Make the colnums array an array of unsigned ints instead of signed ones, and let...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 10 Feb 2000 14:08:20 +0000 (14:08 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 10 Feb 2000 14:08:20 +0000 (14:08 +0000)
git-svn-id: https://svn.dealii.org/trunk@2368 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_constraints.templates.h
deal.II/deal.II/source/dofs/dof_constraints.cc
deal.II/deal.II/source/dofs/dof_renumbering.cc
deal.II/deal.II/source/numerics/matrices.cc
deal.II/lac/Todo
deal.II/lac/include/lac/precondition_block.templates.h
deal.II/lac/include/lac/sparse_ilu.templates.h
deal.II/lac/include/lac/sparse_matrix.h
deal.II/lac/include/lac/sparse_matrix.templates.h
deal.II/lac/include/lac/sparsity_pattern.h
deal.II/lac/source/sparsity_pattern.cc

index 854aa7e35b63a661a377e7ab307c756895a30542..cce8e27ab8c4eb90f23beb0cd7bc210b92308ab2 100644 (file)
@@ -87,7 +87,7 @@ ConstraintMatrix::condense (const SparseMatrix<number> &uncondensed,
                                             // let c point to the constraint
                                             // of this column
            vector<ConstraintLine>::const_iterator c = lines.begin();
-           while (c->line != (unsigned int)uncondensed_struct.get_column_numbers()[j])
+           while (c->line != uncondensed_struct.get_column_numbers()[j])
              ++c;
 
            for (unsigned int q=0; q!=c->entries.size(); ++q)
@@ -117,7 +117,7 @@ ConstraintMatrix::condense (const SparseMatrix<number> &uncondensed,
                                               // let c point to the constraint
                                               // of this column
              vector<ConstraintLine>::const_iterator c = lines.begin();
-             while (c->line != (unsigned int)uncondensed_struct.get_column_numbers()[j])
+             while (c->line != uncondensed_struct.get_column_numbers()[j])
                ++c;
              
              for (unsigned int p=0; p!=c->entries.size(); ++p)
@@ -157,15 +157,15 @@ ConstraintMatrix::condense (SparseMatrix<number> &uncondensed) const
   for (unsigned int c=0; c<lines.size(); ++c)
     distribute[lines[c].line] = static_cast<signed int>(c);
 
-  int n_rows = sparsity.n_rows();
-  for (int row=0; row<n_rows; ++row)
+  const unsigned int n_rows = sparsity.n_rows();
+  for (unsigned int row=0; row<n_rows; ++row)
     {
       if (distribute[row] == -1)
                                         // regular line. loop over cols
        for (unsigned int j=sparsity.get_rowstart_indices()[row];
             j<sparsity.get_rowstart_indices()[row+1]; ++j)
                                           // end of row reached?
-         if (sparsity.get_column_numbers()[j] == -1)
+         if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
            {
                                               // this should not happen, since
                                               // we only operate on compressed
@@ -199,7 +199,7 @@ ConstraintMatrix::condense (SparseMatrix<number> &uncondensed) const
        for (unsigned int j=sparsity.get_rowstart_indices()[row];
             j<sparsity.get_rowstart_indices()[row+1]; ++j)
                                           // end of row reached?
-         if (sparsity.get_column_numbers()[j] == -1)
+         if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
            {
                                               // this should not happen, since
                                               // we only operate on compressed
index 9e1cf51df23d484747f75fbc966f1b0a8e8cc24a..6c6315daf7002b8b1676dc3918e3d1dfe7f7c69e 100644 (file)
@@ -190,7 +190,7 @@ void ConstraintMatrix::condense (const SparsityPattern &uncondensed,
                                             // let c point to the constraint
                                             // of this column
            vector<ConstraintLine>::const_iterator c = lines.begin();
-           while (c->line != (unsigned int)uncondensed.get_column_numbers()[j])
+           while (c->line != uncondensed.get_column_numbers()[j])
              ++c;
 
            for (unsigned int q=0; q!=c->entries.size(); ++q) 
@@ -215,7 +215,7 @@ void ConstraintMatrix::condense (const SparsityPattern &uncondensed,
                                               // let c point to the constraint
                                               // of this column
              vector<ConstraintLine>::const_iterator c = lines.begin();
-             while (c->line != (unsigned int)uncondensed.get_column_numbers()[j]) ++c;
+             while (c->line != uncondensed.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)
@@ -256,10 +256,10 @@ void ConstraintMatrix::condense (SparsityPattern &sparsity) const {
        for (unsigned int j=sparsity.get_rowstart_indices()[row];
             j<sparsity.get_rowstart_indices()[row+1]; ++j)
          {
-           const int column = sparsity.get_column_numbers()[j];
+           const unsigned int column = sparsity.get_column_numbers()[j];
            
                                             // end of row reached?
-           if (column == -1)
+           if (column == SparsityPattern::invalid_entry)
              break;
            else
              if (distribute[column] != -1)
@@ -279,7 +279,7 @@ void ConstraintMatrix::condense (SparsityPattern &sparsity) const {
        for (unsigned int j=sparsity.get_rowstart_indices()[row];
             j<sparsity.get_rowstart_indices()[row+1]; ++j)
                                           // end of row reached?
-         if (sparsity.get_column_numbers()[j] == -1)
+         if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
            break;
          else
            {
index 2e95645a98a43702b36dfc20edfc996f8a26124a..afbd6653fc0ac6cee80775ae9ef0bd679278318c 100644 (file)
@@ -69,7 +69,7 @@ void DoFRenumbering::Cuthill_McKee (DoFHandler<dim>   &dof_handler,
                                           // of this row's entries
          for (j=sparsity.get_rowstart_indices()[row];
               j<sparsity.get_rowstart_indices()[row+1]; ++j)
-           if (sparsity.get_column_numbers()[j] == -1)
+           if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
              break;
                                           // post-condition after loop:
                                           // coordination, i.e. the number
@@ -123,7 +123,7 @@ void DoFRenumbering::Cuthill_McKee (DoFHandler<dim>   &dof_handler,
       for (unsigned int i=0; i<last_round_dofs.size(); ++i)
        for (unsigned int j=sparsity.get_rowstart_indices()[last_round_dofs[i]];
             j<sparsity.get_rowstart_indices()[last_round_dofs[i]+1]; ++j)
-         if (sparsity.get_column_numbers()[j] == -1)
+         if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
            break;
          else
            next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
@@ -163,7 +163,7 @@ void DoFRenumbering::Cuthill_McKee (DoFHandler<dim>   &dof_handler,
          unsigned int coordination = 0;
          for (unsigned int j=sparsity.get_rowstart_indices()[*s];
               j<sparsity.get_rowstart_indices()[*s+1]; ++j)
-           if (sparsity.get_column_numbers()[j] == -1)
+           if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
              break;
            else
              ++coordination;
@@ -248,7 +248,7 @@ void DoFRenumbering::Cuthill_McKee (MGDoFHandler<dim>      &dof_handler,
          unsigned int j;
          for (j=sparsity.get_rowstart_indices()[row];
               j<sparsity.get_rowstart_indices()[row+1]; ++j)
-           if (sparsity.get_column_numbers()[j] == -1)
+           if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
              break;
                                           // post-condition after loop:
                                           // coordination is now
@@ -287,7 +287,7 @@ void DoFRenumbering::Cuthill_McKee (MGDoFHandler<dim>      &dof_handler,
       for (unsigned int i=0; i<last_round_dofs.size(); ++i)
        for (unsigned int j=sparsity.get_rowstart_indices()[last_round_dofs[i]];
             j<sparsity.get_rowstart_indices()[last_round_dofs[i]+1]; ++j)
-         if (sparsity.get_column_numbers()[j] == -1)
+         if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
            break;
          else
            next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
@@ -327,7 +327,7 @@ void DoFRenumbering::Cuthill_McKee (MGDoFHandler<dim>      &dof_handler,
          unsigned int coordination = 0;
          for (unsigned int j=sparsity.get_rowstart_indices()[*s];
               j<sparsity.get_rowstart_indices()[*s+1]; ++j)
-           if (sparsity.get_column_numbers()[j] == -1)
+           if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
              break;
            else
              ++coordination;
index 90c482d3730e0cf58e07f0b4d9edf96c39a5c82f..e03d7ee116d1f662c2db9daed98b12bd3b9f0b12 100644 (file)
@@ -512,11 +512,13 @@ void MatrixTools<dim>::apply_boundary_values (const map<int,double> &boundary_va
   const unsigned int n_dofs             = matrix.m();
   const SparsityPattern    &sparsity    = matrix.get_sparsity_pattern();
   const unsigned int *sparsity_rowstart = sparsity.get_rowstart_indices();
-  const int          *sparsity_colnums  = sparsity.get_column_numbers();
+  const unsigned int *sparsity_colnums  = sparsity.get_column_numbers();
 
   for (; dof != endd; ++dof)
     {
-      const int dof_number = dof->first;
+      Assert (dof->first >= 0, ExcInternalError());
+      
+      const unsigned int dof_number = static_cast<unsigned int>(dof->first);
                                       // for each boundary dof:
       
                                       // set entries of this line
@@ -577,7 +579,7 @@ void MatrixTools<dim>::apply_boundary_values (const map<int,double> &boundary_va
        {
                                           // we need not handle the
                                           // row we have already cleared
-         if ((signed int)row == dof_number)
+         if (row == dof_number)
            continue;
 
                                           // check whether the line has
@@ -592,9 +594,9 @@ void MatrixTools<dim>::apply_boundary_values (const map<int,double> &boundary_va
                                           //
                                           // if this row contains an element
                                           // for this dof, *p==dof_number
-         const int * p = lower_bound(&sparsity_colnums[sparsity_rowstart[row]+1],
-                                     &sparsity_colnums[sparsity_rowstart[row+1]],
-                                     dof_number);
+         const unsigned int * p = lower_bound(&sparsity_colnums[sparsity_rowstart[row]+1],
+                                              &sparsity_colnums[sparsity_rowstart[row+1]],
+                                              dof_number);
                                           // check whether this line has
                                           // an entry in the regarding column
                                           // (check for ==dof_number and
index bf2b8a23f49bf553b51556276058756f333144a4..5a69843553ea55cb57f0692460fbfe93ca239653 100644 (file)
@@ -2,13 +2,6 @@
 ====== Things to be done for the lac library     ================
 =================================================================
 
-Use unsigned integers for the colnums array in dSMatrixStruct. This
-  would enhance safety since colnum=-1 would no longer point to a
-  valid address. Mark unused values by numeric_limits<uint>::max
-  (or some such constant) instead, which probably has the same
-  bit representation as (uint)-1, but is an unsigned, therefore
-  eliminating the need for conversions.
-
 
 Make ILU decomposition faster by following the comment in the
   innermost loop (about the better use of the index j).
index 386549ac3f039e90ad5044a0b251cc7d642339bf..6f373177be998429a81c1bdb9c61b8fdc8476a06 100644 (file)
@@ -135,7 +135,7 @@ PreconditionBlockSOR<number,inverse_type>::~PreconditionBlockSOR(){}
 
 template <typename number, typename inverse_type>
 template <typename number2>
-void PreconditionBlockSOR<number,inverse_type>::operator() (Vector<number2> &dst,
+void PreconditionBlockSOR<number,inverse_type>::operator() (Vector<number2>       &dst,
                                                            const Vector<number2> &src) const
 {
   Assert(A!=0, ExcNoMatrixGivenToUse());
@@ -149,7 +149,7 @@ void PreconditionBlockSOR<number,inverse_type>::operator() (Vector<number2> &dst
 
   const SparsityPattern &spars    = M.get_sparsity_pattern();
   const unsigned int    *rowstart = spars.get_rowstart_indices();
-  const int             *columns  = spars.get_column_numbers();
+  const unsigned int    *columns  = spars.get_column_numbers();
 
   Vector<number2> b_cell(blocksize), x_cell(blocksize);
 
@@ -172,8 +172,7 @@ void PreconditionBlockSOR<number,inverse_type>::operator() (Vector<number2> &dst
            {
              b_cell_row=src(row);
              for (unsigned int j=rowstart[row]; j<rowstart[row+1]; ++j)
-               if ((column=static_cast<unsigned int>(columns[j]))
-                   < begin_diag_block)
+               if ((column=columns[j]) < begin_diag_block)
                    b_cell_row -= M.global_entry(j) * dst(column);
              b_cell(row_cell)=b_cell_row;
              for (unsigned int column_cell=0, column=cell*blocksize;
@@ -196,7 +195,7 @@ void PreconditionBlockSOR<number,inverse_type>::operator() (Vector<number2> &dst
          {
            b_cell_row=src(row);
            for (unsigned int j=rowstart[row]; j<rowstart[row+1]; ++j)
-             if ((column=static_cast<unsigned int>(columns[j])) < begin_diag_block)
+             if ((column=columns[j]) < begin_diag_block)
                {
                  b_cell_row -= M.global_entry(j) * dst(column);
                }
index c7b91ecb97612536e207fe2ad46be93d74cc6f2b..e3569dbdf2c6c7bd0c7dda6c060cd973083028e4 100644 (file)
@@ -73,11 +73,11 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
                                       // pattern of the old matrix!
       const unsigned int * const rowstart_indices
        = matrix.get_sparsity_pattern().get_rowstart_indices();
-      const int * const          column_numbers
+      const unsigned int * const column_numbers
        = matrix.get_sparsity_pattern().get_column_numbers();
       
       for (unsigned int row=0; row<m(); ++row)
-       for (const int * col = &column_numbers[rowstart_indices[row]];
+       for (const unsigned int * col = &column_numbers[rowstart_indices[row]];
             col != &column_numbers[rowstart_indices[row+1]]; ++col)
          set (row, *col, matrix.global_entry(col-column_numbers));
     };
@@ -110,7 +110,7 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
                                   // matrix
   const SparsityPattern             &sparsity = get_sparsity_pattern();
   const unsigned int * const rowstart_indices = sparsity.get_rowstart_indices();
-  const int * const          column_numbers   = sparsity.get_column_numbers();
+  const unsigned int * const column_numbers   = sparsity.get_column_numbers();
   
 /*
   PSEUDO-ALGORITHM
@@ -148,15 +148,15 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
                                       // let k run over all lower-left
                                       // elements of row i; skip
                                       // diagonal element at start
-      const int * first_of_row
+      const unsigned int * first_of_row
        = &column_numbers[rowstart_indices[row]+1];
-      const int * first_after_diagonal
+      const unsigned int * first_after_diagonal
        = lower_bound (&column_numbers[rowstart_indices[row]+1],
                       &column_numbers[rowstart_indices[row+1]],
-                      static_cast<int>(row));
+                      row);
 
                                       // k := *col_ptr
-      for (const int * col_ptr = first_of_row; col_ptr!=first_after_diagonal; ++col_ptr)
+      for (const unsigned int * col_ptr = first_of_row; col_ptr!=first_after_diagonal; ++col_ptr)
        {
          const unsigned int global_index_ik = col_ptr-column_numbers;
          global_entry(global_index_ik) *= diag_element(*col_ptr);
@@ -169,7 +169,7 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
                                           // row correctly, but excluding
                                           // the main diagonal entry
          bool left_of_diagonal = true;
-         for (const int * j = col_ptr+1;
+         for (const unsigned int * j = col_ptr+1;
               j<&column_numbers[rowstart_indices[row+1]];
               ++j)
            {
@@ -189,9 +189,9 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
                                               //
                                               // note that j should never point
                                               // to the diagonal itself!
-             if (left_of_diagonal && (*j > static_cast<int>(row)))
+             if (left_of_diagonal && (*j > row))
                {
-                 Assert (*j != static_cast<int>(row), ExcInternalError());
+                 Assert (*j != row, ExcInternalError());
                  
                  left_of_diagonal = false;
                                                   // a[i,i] -= a[i,k]*a[k,i]
@@ -271,10 +271,8 @@ void SparseILU<number>::apply_decomposition (Vector<somenumber>       &dst,
   Assert (dst.size() == m(), ExcSizeMismatch(dst.size(), m()));
   
   const unsigned int N=dst.size();
-  const unsigned int * const rowstart_indices
-    = get_sparsity_pattern().get_rowstart_indices();
-  const int * const          column_numbers
-    = get_sparsity_pattern().get_column_numbers();
+  const unsigned int * const rowstart_indices = get_sparsity_pattern().get_rowstart_indices();
+  const unsigned int * const column_numbers   = get_sparsity_pattern().get_column_numbers();
                                   // solve LUx=b in two steps:
                                   // first Ly = b, then
                                   //       Ux = y
@@ -292,15 +290,15 @@ void SparseILU<number>::apply_decomposition (Vector<somenumber>       &dst,
     {
                                       // get start of this row. skip the
                                       // diagonal element
-      const int * const rowstart = &column_numbers[rowstart_indices[row]+1];
+      const unsigned int * const rowstart = &column_numbers[rowstart_indices[row]+1];
                                       // find the position where the part
                                       // right of the diagonal starts
-      const int * const first_after_diagonal
+      const unsigned int * const first_after_diagonal
        = lower_bound (rowstart,
                       &column_numbers[rowstart_indices[row+1]],
-                      static_cast<int>(row));
+                      row);
       
-      for (const int * col=rowstart; col!=first_after_diagonal; ++col)
+      for (const unsigned int * col=rowstart; col!=first_after_diagonal; ++col)
        dst(row) -= global_entry (col-column_numbers) * dst(*col);
     };
 
@@ -315,15 +313,15 @@ void SparseILU<number>::apply_decomposition (Vector<somenumber>       &dst,
   for (int row=N-1; row>=0; --row)
     {
                                       // get end of this row
-      const int * const rowend = &column_numbers[rowstart_indices[row+1]];
+      const unsigned int * const rowend = &column_numbers[rowstart_indices[row+1]];
                                       // find the position where the part
                                       // right of the diagonal starts
-      const int * const first_after_diagonal
+      const unsigned int * const first_after_diagonal
        = lower_bound (&column_numbers[rowstart_indices[row]+1],
                       &column_numbers[rowstart_indices[row+1]],
-                      row);
+                      static_cast<unsigned int>(row));
       
-      for (const int * col=first_after_diagonal; col!=rowend; ++col)
+      for (const unsigned int * col=first_after_diagonal; col!=rowend; ++col)
        dst(row) -= global_entry (col-column_numbers) * dst(*col);
 
                                       // scale by the diagonal element.
index d93459ae6705c2eafe3987913a234864e0fc4e34..350bae412fa052c84af6edc8c85818c91f4ba4f6 100644 (file)
@@ -680,16 +680,17 @@ unsigned int SparseMatrix<number>::n () const
 
 template <typename number>
 inline
-void SparseMatrix<number>::set (const unsigned int i, const unsigned int j,
-                               const number value) {
+void SparseMatrix<number>::set (const unsigned int i,
+                               const unsigned int j,
+                               const number value)
+{
   Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert ((cols->operator()(i,j) != -1) || (value == 0.),
+  Assert ((cols->operator()(i,j) != SparsityPattern::invalid_entry) ||
+         (value == 0.),
          ExcInvalidIndex(i,j));
 
-  const int index = cols->operator()(i,j);
-
-  if (index >= 0)
-    val[index] = value;
+  if (value != 0.)
+    val[cols->operator()(i,j)] = value;
 };
 
 
@@ -699,13 +700,12 @@ inline
 void SparseMatrix<number>::add (const unsigned int i, const unsigned int j,
                                const number value) {
   Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert ((cols->operator()(i,j) != -1) || (value == 0.),
+  Assert ((cols->operator()(i,j) != SparsityPattern::invalid_entry) ||
+         (value == 0.),
          ExcInvalidIndex(i,j));
 
-  const int index = cols->operator()(i,j);
-  
-  if (index >= 0)
-    val[index] += value;
+  if (value != 0.)
+    val[cols->operator()(i,j)] += value;
 };
 
 
@@ -714,9 +714,11 @@ void SparseMatrix<number>::add (const unsigned int i, const unsigned int j,
 
 template <typename number>
 inline
-number SparseMatrix<number>::operator () (const unsigned int i, const unsigned int j) const {
+number SparseMatrix<number>::operator () (const unsigned int i,
+                                         const unsigned int j) const
+{
   Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (cols->operator()(i,j) != -1,
+  Assert (cols->operator()(i,j) != SparsityPattern::invalid_entry,
          ExcInvalidIndex(i,j));
   return val[cols->operator()(i,j)];
 };
@@ -725,7 +727,8 @@ number SparseMatrix<number>::operator () (const unsigned int i, const unsigned i
 
 template <typename number>
 inline
-number SparseMatrix<number>::diag_element (const unsigned int i) const {
+number SparseMatrix<number>::diag_element (const unsigned int i) const
+{
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert (m() == n(), ExcMatrixNotSquare());
   Assert (i<max_len, ExcInvalidIndex1(i));
@@ -757,7 +760,8 @@ number & SparseMatrix<number>::diag_element (const unsigned int i)
 template <typename number>
 inline
 number
-SparseMatrix<number>::raw_entry(unsigned int row, unsigned int index) const
+SparseMatrix<number>::raw_entry (const unsigned int row,
+                                const unsigned int index) const
 {
   Assert(row<cols->rows, ExcIndexRange(row,0,cols->rows));
   Assert(index<cols->row_length(row), ExcIndexRange(index,0,cols->row_length(row)));
@@ -769,7 +773,8 @@ SparseMatrix<number>::raw_entry(unsigned int row, unsigned int index) const
 
 template <typename number>
 inline
-number SparseMatrix<number>::global_entry (const unsigned int j) const {
+number SparseMatrix<number>::global_entry (const unsigned int j) const
+{
   Assert (cols != 0, ExcMatrixNotInitialized());
   return val[j];
 };
@@ -778,7 +783,8 @@ number SparseMatrix<number>::global_entry (const unsigned int j) const {
 
 template <typename number>
 inline
-number & SparseMatrix<number>::global_entry (const unsigned int j) {
+number & SparseMatrix<number>::global_entry (const unsigned int j)
+{
   Assert (cols != 0, ExcMatrixNotInitialized());
   return val[j];
 };
index b7364fe005e6d7ee9e388bd45d224cae42d31e2d..d80666f0204706f1bd41c59dd877e9e07bd052a5 100644 (file)
@@ -260,9 +260,9 @@ SparseMatrix<number>::vmult (Vector<somenumber>& dst, const Vector<somenumber>&
 
                                   // if not in MT mode or size<2000
                                   // do it in an oldfashioned way
-  const number *val_ptr = &val[cols->rowstart[0]];
-  const int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
-  somenumber   *dst_ptr = &dst(0);
+  const number       *val_ptr    = &val[cols->rowstart[0]];
+  const unsigned int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
+  somenumber         *dst_ptr    = &dst(0);
   for (unsigned int row=0; row<n_rows; ++row)
     {
       somenumber s = 0.;
@@ -284,9 +284,9 @@ SparseMatrix<number>::threaded_vmult (Vector<somenumber>       &dst,
                                      const unsigned int        end_row) const
 {
 #ifdef DEAL_II_USE_MT
-  const number *val_ptr = &val[cols->rowstart[begin_row]];
-  const int *colnum_ptr = &cols->colnums[cols->rowstart[begin_row]];
-  somenumber   *dst_ptr = &dst(begin_row);
+  const number       *val_ptr    = &val[cols->rowstart[begin_row]];
+  const unsigned int *colnum_ptr = &cols->colnums[cols->rowstart[begin_row]];
+  somenumber         *dst_ptr    = &dst(begin_row);
   for (unsigned int row=begin_row; row<end_row; ++row)
     {
       somenumber s = 0.;
@@ -322,7 +322,7 @@ SparseMatrix<number>::Tvmult (Vector<somenumber>& dst, const Vector<somenumber>&
     {
       for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
        {
-         int p = cols->colnums[j];
+         const unsigned int p = cols->colnums[j];
          dst(p) += val[j] * src(i);
        }
     }
@@ -340,10 +340,10 @@ SparseMatrix<number>::vmult_add (Vector<somenumber>& dst, const Vector<somenumbe
   Assert(m() == dst.size(), ExcDimensionsDontMatch(m(),dst.size()));
   Assert(n() == src.size(), ExcDimensionsDontMatch(n(),src.size()));
 
-  const unsigned int n_rows = m();
-  const number *val_ptr = &val[cols->rowstart[0]];
-  const int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
-  somenumber   *dst_ptr = &dst(0);
+  const unsigned int  n_rows     = m();
+  const number       *val_ptr    = &val[cols->rowstart[0]];
+  const unsigned int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
+  somenumber         *dst_ptr    = &dst(0);
   for (unsigned int row=0; row<n_rows; ++row)
     {
       somenumber s = 0.;
@@ -369,7 +369,7 @@ SparseMatrix<number>::Tvmult_add (Vector<somenumber>& dst, const Vector<somenumb
     {
       for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
        {
-         int p = cols->colnums[j];
+         const unsigned int p = cols->colnums[j];
          dst(p) += val[j] * src(i);
        }
     }
@@ -438,9 +438,9 @@ SparseMatrix<number>::matrix_norm (const Vector<somenumber>& v) const
 #endif
                                   // if not in MT mode or the matrix is
                                   // too small: do it one-by-one
-  somenumber sum = 0.;
-  const number *val_ptr = &val[cols->rowstart[0]];
-  const int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
+  somenumber          sum        = 0.;
+  const number       *val_ptr    = &val[cols->rowstart[0]];
+  const unsigned int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
   for (unsigned int row=0; row<n_rows; ++row)
     {
       somenumber s = 0.;
@@ -466,8 +466,8 @@ SparseMatrix<number>::threaded_matrix_norm (const Vector<somenumber> &v,
 {
 #ifdef DEAL_II_USE_MT
   somenumber sum = 0.;
-  const number *val_ptr = &val[cols->rowstart[begin_row]];
-  const int *colnum_ptr = &cols->colnums[cols->rowstart[begin_row]];
+  const number      *val_ptr    = &val[cols->rowstart[begin_row]];
+  const unsignedint *colnum_ptr = &cols->colnums[cols->rowstart[begin_row]];
   for (unsigned int row=begin_row; row<end_row; ++row)
     {
       somenumber s = 0.;
@@ -600,7 +600,7 @@ SparseMatrix<number>::residual (Vector<somenumber>       &dst,
       somenumber s = b(i);
       for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
        {
-         int p = cols->colnums[j];
+         const unsigned int p = cols->colnums[j];
          s -= val[j] * u(p);
        }
       dst(i) = s;
@@ -628,7 +628,7 @@ SparseMatrix<number>::threaded_residual (Vector<somenumber>       &dst,
       somenumber s = b(i);
       for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
        {
-         int p = cols->colnums[j];
+         const unsigned int p = cols->colnums[j];
          s -= val[j] * u(p);
        }
       dst(i) = s;
@@ -704,7 +704,7 @@ SparseMatrix<number>::precondition_SSOR (Vector<somenumber>& dst,
       const unsigned int first_right_of_diagonal_index
        = (lower_bound (&cols->colnums[*rowstart_ptr+1],
                        &cols->colnums[*(rowstart_ptr+1)],
-                       static_cast<signed int>(row)) -
+                       row) -
           &cols->colnums[0]);
                                       
       for (unsigned int j=(*rowstart_ptr)+1; j<first_right_of_diagonal_index; ++j)
@@ -727,10 +727,10 @@ SparseMatrix<number>::precondition_SSOR (Vector<somenumber>& dst,
       const unsigned int first_right_of_diagonal_index
        = (lower_bound (&cols->colnums[*rowstart_ptr+1],
                        &cols->colnums[*(rowstart_ptr+1)],
-                       static_cast<signed int>(row)) -
+                       static_cast<unsigned int>(row)) -
           &cols->colnums[0]);
       for (unsigned int j=first_right_of_diagonal_index; j<*(rowstart_ptr+1); ++j)
-       if (cols->colnums[j] > row)
+       if (cols->colnums[j] > static_cast<unsigned int>(row))
          *dst_ptr -= om* val[j] * dst(cols->colnums[j]);
       
       *dst_ptr /= val[*rowstart_ptr];
@@ -769,7 +769,7 @@ SparseMatrix<number>::SOR (Vector<somenumber>& dst, const number om) const
     {
       somenumber s = dst(row);
       for (unsigned int j=cols->rowstart[row]; j<cols->rowstart[row+1]; ++j)
-       if (static_cast<unsigned int>(cols->colnums[j]) < row)
+       if (cols->colnums[j] < row)
          s -= val[j] * dst(cols->colnums[j]);
 
       dst(row) = s * om / val[cols->rowstart[row]];
@@ -847,7 +847,6 @@ SparseMatrix<number>::SSOR (Vector<somenumber>& dst, const number om) const
   Assert (m() == n(), ExcMatrixNotSquare());
   Assert (m() == dst.size(), ExcDimensionsDontMatch(m(),dst.size()));
 
-  int p;
   const unsigned int  n = dst.size();
   unsigned int  j;
   somenumber s;
@@ -857,8 +856,8 @@ SparseMatrix<number>::SSOR (Vector<somenumber>& dst, const number om) const
       s = 0.;
       for (j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
        {
-         p = cols->colnums[j];
-         if (p>=0)
+         const unsigned int p = cols->colnums[j];
+         if (p != SparsityPattern::invalid_entry)
            {
              if (i>j) s += val[j] * dst(p);
            }
@@ -872,8 +871,8 @@ SparseMatrix<number>::SSOR (Vector<somenumber>& dst, const number om) const
       s = 0.;
       for (j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
        {
-         p = cols->colnums[j];
-         if (p>=0)
+         const unsigned int p = cols->colnums[j];
+         if (p != SparsityPattern::invalid_entry)
            {
              if (static_cast<unsigned int>(i)<j) s += val[j] * dst(p);
            }
@@ -933,7 +932,7 @@ void SparseMatrix<number>::print_formatted (ostream &out,
   for (unsigned int i=0; i<m(); ++i)
     {
       for (unsigned int j=0; j<n(); ++j)
-       if ((*cols)(i,j) != -1)
+       if ((*cols)(i,j) != SparsityPattern::invalid_entry)
          out << setw(width)
              << val[cols->operator()(i,j)] << ' ';
        else
index 41a65bb1c735b0df458c59b48f99d268ec8f0b5f..992a1e48c4f159e0a5ae7d534754e06370b2cd7c 100644 (file)
  * needed. Then #a_{24}=val[number of the found element] = 3#.
  *
  *
- * @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth; some documentation by Rasched Zamni
+ * @author Wolfgang Bangerth and others
  */
 class SparsityPattern : public Subscriptor
 {
   public:
+                                    /**
+                                     * Define a value which is used
+                                     * to indicate that a certain
+                                     * value in the #colnums# array
+                                     * is unused, i.e. does not
+                                     * represent a certain column
+                                     * number index.
+                                     *
+                                     * Indices with this invalid
+                                     * value are used to insert new
+                                     * entries to the sparsity
+                                     * pattern using the #add# member
+                                     * function, and are removed when
+                                     * calling #compress#.
+                                     *
+                                     * You should not assume that the
+                                     * variable declared here has a
+                                     * certain value. The
+                                     * initialization is given here
+                                     * only to enable the compiler to
+                                     * perform some optimizations,
+                                     * but the actual value of the
+                                     * variable may change over time.
+                                     */
+    static const unsigned int invalid_entry = static_cast<unsigned int>(-1);
+    
                                     /**
                                      * Initialize the matrix empty,
                                      * that is with no memory
@@ -305,7 +331,7 @@ class SparsityPattern : public Subscriptor
                                      * element with row number #i# and
                                      * column number #j#. If the matrix
                                      * element is not a nonzero one,
-                                     * return -1.
+                                     * return #SparsityPattern::invalid_entry#.
                                      *
                                      * This function is usually called
                                      * by the #operator()# of the
@@ -318,7 +344,7 @@ class SparsityPattern : public Subscriptor
                                      * because the column numbers are
                                      * sorted.
                                      */
-    int operator() (const unsigned int i, const unsigned int j) const;
+    unsigned int operator() (const unsigned int i, const unsigned int j) const;
 
                                     /**
                                      * Add a nonzero entry to the matrix.
@@ -330,25 +356,6 @@ class SparsityPattern : public Subscriptor
                                      */
     void add (const unsigned int i, const unsigned int j);
     
-                                    /**
-                                     * This matrix adds a whole connectivity
-                                     * list to the sparsity structure
-                                     * respresented by this object. It assumes
-                                     * the #rowcols# array to be a list of
-                                     * indices which are all linked together,
-                                     * i.e. all entries
-                                     * #(rowcols[i], rowcols[j])# for all
-                                     * #i,j=0...n# for this sparsity pattern
-                                     * are created. #n# is assumed to be the
-                                     * number of elements pointed to by
-                                     * #rowcols#.
-                                     */
-    void add_matrix (const unsigned int n, const int* rowcols);
-
-                                    //////////
-    void add_matrix (const unsigned int m, const unsigned int n,
-                    const int* rows, const int* cols);
-
                                     /**
                                      * Print the sparsity of the matrix
                                      * in a format that #gnuplot# understands
@@ -483,7 +490,7 @@ class SparsityPattern : public Subscriptor
                                      * avoid programs relying on outdated
                                      * information!
                                      */
-    const int * get_column_numbers () const;
+    const unsigned int * get_column_numbers () const;
     
     
                                     /**
@@ -655,7 +662,7 @@ class SparsityPattern : public Subscriptor
                                      * element, i.e.
                                      * #colnums[rowstart[r]]==r#.
                                      */
-    int *colnums;
+    unsigned int *colnums;
 
                                     /**
                                      * Store whether the #compress# function
@@ -711,13 +718,14 @@ SparsityPattern::get_rowstart_indices () const
 
 
 inline
-const int *
+const unsigned int *
 SparsityPattern::get_column_numbers () const
 {
   return colnums;
 };
 
 
+
 inline
 unsigned int
 SparsityPattern::row_length (const unsigned int row) const
@@ -726,10 +734,12 @@ SparsityPattern::row_length (const unsigned int row) const
   return rowstart[row+1]-rowstart[row];
 }
 
+
+
 inline
 unsigned int
 SparsityPattern::column_number (const unsigned int row,
-                                  const unsigned int index) const
+                               const unsigned int index) const
 {
   Assert(row<rows, ExcIndexRange(row,0,rows));
   Assert(index<row_length(row), ExcIndexRange(index,0,row_length(row)));
@@ -738,6 +748,7 @@ SparsityPattern::column_number (const unsigned int row,
 }
 
 
+
 inline
 unsigned int
 SparsityPattern::n_nonzero_elements () const
index 6785747d26217a8d90f04b99fcd63d5cf292c5d9..09170a04590f749389a08547e18c9c0be0aa9195 100644 (file)
@@ -92,6 +92,7 @@ SparsityPattern::SparsityPattern (const unsigned int          m,
 };
 
 
+
 SparsityPattern::SparsityPattern (const SparsityPattern &original,
                                  const unsigned int        max_per_row,
                                  const unsigned int        extra_off_diagonals)
@@ -126,28 +127,24 @@ SparsityPattern::SparsityPattern (const SparsityPattern &original,
                                       // the rest from that element onwards
                                       // which is not a side-diagonal any
                                       // more.
-      const int * const
+      const unsigned int * const
        original_row_start = &original.colnums[original.rowstart[row]] + 1;
-                                      // thw following requires that
+                                      // the following requires that
                                       // #original# be compressed since
-                                      // otherwise there might be -1's
-      const int * const
+                                      // otherwise there might be invalid_entry's
+      const unsigned int * const
        original_row_end   = &original.colnums[original.rowstart[row+1]];
 
-      const int * const
+      const unsigned int * const
        original_last_before_side_diagonals = lower_bound (original_row_start,
                                                           original_row_end,
-                                                          static_cast<int>
-                                                          (row
-                                                           -extra_off_diagonals));
-      const int * const
+                                                          row-extra_off_diagonals);
+      const unsigned int * const
        original_first_after_side_diagonals = upper_bound (original_row_start,
                                                           original_row_end,
-                                                          static_cast<int>
-                                                          (row
-                                                           +extra_off_diagonals));
+                                                          row+extra_off_diagonals);
 
-      int * next_free_slot = &colnums[rowstart[row]] + 1;
+      unsigned int * next_free_slot = &colnums[rowstart[row]] + 1;
 
                                       // copy elements before side-diagonals
       next_free_slot = copy (original_row_start,
@@ -274,7 +271,7 @@ SparsityPattern::reinit (const unsigned int m,
     {
       if (colnums) delete[] colnums;
       max_vec_len = vec_len;
-      colnums = new int[max_vec_len];
+      colnums = new unsigned int[max_vec_len];
     };
 
                                   // set the rowstart array 
@@ -286,7 +283,7 @@ SparsityPattern::reinit (const unsigned int m,
                                   // preset the column numbers by a
                                   // value indicating it is not in
                                   // use
-  fill_n (&colnums[0], vec_len, -1);
+  fill_n (&colnums[0], vec_len, invalid_entry);
 
                                   // if the matrix is square: let the
                                   // first entry in each row be the
@@ -303,11 +300,14 @@ void
 SparsityPattern::compress ()
 {
   Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
+
+                                  // do nothing if already compressed
+  if (compressed)
+    return;
   
-  if (compressed) return;
   unsigned int next_free_entry = 0,
-               next_row_start = 0,
-                   row_length = 0;
+              next_row_start  = 0,
+              row_length      = 0;
 
                                   // first find out how many non-zero
                                   // elements there are, in order to
@@ -316,14 +316,14 @@ SparsityPattern::compress ()
   const unsigned int
     nonzero_elements = count_if (&colnums[rowstart[0]],
                                 &colnums[rowstart[rows]],
-                                bind2nd(not_equal_to<int>(), -1));
+                                bind2nd(not_equal_to<unsigned int>(), invalid_entry));
                                   // now allocate the respective memory
-  int *new_colnums = new int[nonzero_elements];
+  unsigned int *new_colnums = new unsigned int[nonzero_elements];
 
   
                                   // reserve temporary storage to
                                   // store the entries of one row
-  vector<int> tmp_entries (max_row_length);
+  vector<unsigned int> tmp_entries (max_row_length);
   
                                   // Traverse all rows
   for (unsigned int line=0; line<rows; ++line)
@@ -332,7 +332,7 @@ SparsityPattern::compress ()
                                       // first unused entry is reached
       row_length = 0;
       for (unsigned int j=rowstart[line]; j<rowstart[line+1]; ++j,++row_length)
-       if (colnums[j] != -1)
+       if (colnums[j] != invalid_entry)
          tmp_entries[row_length] = colnums[j];
        else
          break;
@@ -361,7 +361,7 @@ SparsityPattern::compress ()
 
                                       // some internal checks
       Assert ((rows!=cols) ||
-             (new_colnums[rowstart[line]] == static_cast<signed int>(line)),
+             (new_colnums[rowstart[line]] == line),
              ExcInternalError());
                                       // assert that the first entry
                                       // does not show up in
@@ -452,8 +452,9 @@ SparsityPattern::max_entries_per_row () const
 
 
 
-int
-SparsityPattern::operator () (const unsigned int i, const unsigned int j) const
+unsigned int
+SparsityPattern::operator () (const unsigned int i,
+                             const unsigned int j) const
 {
   Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
   Assert (i<rows, ExcInvalidIndex(i,rows));
@@ -466,13 +467,13 @@ SparsityPattern::operator () (const unsigned int i, const unsigned int j) const
                                   // if a first entry exists)
   if (rowstart[i] != rowstart[i+1]) 
     {
-      if (static_cast<signed int>(j) == colnums[rowstart[i]])
+      if (j == colnums[rowstart[i]])
        return rowstart[i];
     }
   else
                                     // no first entry exists for this
                                     // line
-    return -1;
+    return invalid_entry;
 
                                   // all other entries are sorted, so
                                   // we can use a binary seach algorithm
@@ -484,19 +485,21 @@ SparsityPattern::operator () (const unsigned int i, const unsigned int j) const
                                   // at the top of this function, so it
                                   // may not be called for noncompressed
                                   // structures.
-  const int * const p = lower_bound (&colnums[rowstart[i]+1],
-                                    &colnums[rowstart[i+1]],
-                                    static_cast<signed int>(j));
-  if ((*p == static_cast<signed int>(j)) &&
+  const unsigned int * const p = lower_bound (&colnums[rowstart[i]+1],
+                                             &colnums[rowstart[i+1]],
+                                             j);
+  if ((*p == j) &&
       (p != &colnums[rowstart[i+1]]))
     return (p - &colnums[0]);
   else
-    return -1;
+    return invalid_entry;
 }
 
 
+
 void
-SparsityPattern::add (const unsigned int i, const unsigned int j)
+SparsityPattern::add (const unsigned int i,
+                     const unsigned int j)
 {
   Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
   Assert (i<rows, ExcInvalidIndex(i,rows));
@@ -506,10 +509,10 @@ SparsityPattern::add (const unsigned int i, const unsigned int j)
   for (unsigned int k=rowstart[i]; k<rowstart[i+1]; k++)
     {
                                       // entry already exists
-      if (colnums[k] == (signed int)j) return;
+      if (colnums[k] == j) return;
                                       // empty entry found, put new
                                       // entry here
-      if (colnums[k] == -1)
+      if (colnums[k] == invalid_entry)
        {
          colnums[k] = j;
          return;
@@ -524,36 +527,13 @@ SparsityPattern::add (const unsigned int i, const unsigned int j)
 
 
 
-void
-SparsityPattern::add_matrix (const unsigned int n, const int* rowcols)
-{
-  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
-  for (unsigned int i=0; i<n; ++i)
-    for (unsigned int j=0; j<n; ++j)
-      add(rowcols[i], rowcols[j]);
-}
-
-
-
-void
-SparsityPattern::add_matrix (const unsigned int m, const unsigned int n,
-                            const int* rows, const int* cols)
-{
-  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
-  for (unsigned i=0; i<m; ++i)
-    for (unsigned j=0; j<n; ++j)
-      add(rows[i], cols[j]);
-}
-
-
-
 void
 SparsityPattern::print_gnuplot (ostream &out) const
 {
   Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
   for (unsigned int i=0; i<rows; ++i)
     for (unsigned int j=rowstart[i]; j<rowstart[i+1]; ++j)
-      if (colnums[j]>=0)
+      if (colnums[j] != invalid_entry)
        out << i << " " << -colnums[j] << endl;
 
   AssertThrow (out, ExcIO());
@@ -568,7 +548,7 @@ SparsityPattern::bandwidth () const
   unsigned int b=0;
   for (unsigned int i=0; i<rows; ++i)
     for (unsigned int j=rowstart[i]; j<rowstart[i+1]; ++j)
-      if (colnums[j]>=0) 
+      if (colnums[j] != invalid_entry)
        {
          if (static_cast<unsigned int>(abs(static_cast<int>(i-colnums[j]))) > b)
            b = abs(static_cast<int>(i-colnums[j]));

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.