From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d> Date: Mon, 4 May 1998 08:36:03 +0000 (+0000) Subject: Optimisations in the compress function by using standard algorithms and thinking... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b1d10f5db1984474a3845b77692c71ab6a88f034;p=dealii-svn.git Optimisations in the compress function by using standard algorithms and thinking about what needs sorting and copying and what not. git-svn-id: https://svn.dealii.org/trunk@235 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/source/dsmatrix.cc b/deal.II/lac/source/dsmatrix.cc index 2c56b9bf8a..6ca5d49372 100644 --- a/deal.II/lac/source/dsmatrix.cc +++ b/deal.II/lac/source/dsmatrix.cc @@ -7,112 +7,7 @@ #include <lac/dsmatrix.h> #include <iostream> #include <iomanip> - -/*----------------- from sort.h -------------------------*/ - - -////////// -template<class T> -inline void swap (T* a, T* b) -{ - T x = *a; - *a = *b; - *b = x; -} - -////////// -template<class T> -inline void simple_sort (const long n, T* field) -{ - long i,j; - for (i=1;i<n;i++) - { - for (j=i+1;j<=n;j++) - { - if (field[j] < field[i]) - { - swap(field+i,field+j); - } - } - } -} - -template<class T> -inline void heapsort_sift (T* a, const long l, const long r) -{ - long i = l; - long j = 2*i; - T x = a[i]; - - while (j<=r) - { - if (j<r) - { - if (a[j] < a[j+1]) j++; - } - if (!(x < a[j])) break; - a[i] = a[j]; - i = j; - j = 2*i; - } - a[i] = x; -} - - -////////// -template<class T> -inline void heapsort (const int n, T* field) -{ - field--; - - long l =(n/2)+1; - long r = n; - - while (l>1) - { - l--; - heapsort_sift(field,l,r); - } - while (r>1) - { - swap(field+l,field+r); - r--; - heapsort_sift(field,l,r); - } -} - -////////// -template<class T> -inline void _quicksort (const long r, T* a, const long l) -{ - long i = l; - long j = r; - T* x = &a[(l+r)/2]; - do - { - while (a[i] < *x) i++; - while (*x < a[j]) j--; - if (i<=j) - { - swap(a+i,a+j); - i++; - j--; - } - } - while (i<=j); - if (l<j) _quicksort(j,a,l); - if (i<r) _quicksort(r,a,i); -} - -template<class T> -inline void quicksort (const long r, T* a) -{ - _quicksort(r,a,1); -} - - -/*----------------- end: from sort.h -------------------------*/ - +#include <algorithm> @@ -179,14 +74,13 @@ dSMatrixStruct::reinit (const unsigned int m, const unsigned int n, for (unsigned int i=0; i<=rows; i++) rowstart[i] = i * max_per_row; - for (int i = vec_len-1; i>=0; i--) - colnums[i] = -1; + fill_n (&colnums[0], vec_len, -1); if (rows == cols) for (unsigned int i=0;i<rows;i++) colnums[rowstart[i]] = i; - compressed = 0; + compressed = false; } @@ -194,41 +88,76 @@ void dSMatrixStruct::compress () { if (compressed) return; - unsigned int i,j,k,l,s; - int* entries; - - entries = new int[max_row_len]; + unsigned int next_free_entry = 0, + next_row_start = 0, + row_length = 0; - // Traverse all rows + // reserve temporary storage to + // store the entries of one wor + int *tmp_entries = new int[max_row_len]; - for (i=0,k=0,s=0 ; i<rows ; ++i) + // Traverse all rows + for (unsigned int line=0; line<rows; ++line) { - - // Sort entries in ascending order - - for (j=rowstart[i],l=0; j<rowstart[i+1]; ++j,++l) - entries[l] = colnums[j]; - heapsort(max_row_len, entries); - - // Re-insert column numbers into the field - - // Ensure diagonal entry first in quadratic matrix - - if (cols == rows) colnums[k++] = i; - - for (l=0; l < max_row_len ; ++l) - if (entries[l]>=0) - if (((signed int)i!=entries[l]) || (rows!=cols)) - colnums[k++] = entries[l]; - - rowstart[i] = s; - s = k; - } - vec_len = rowstart[i] = s; + // copy used entries, break if + // 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) + tmp_entries[row_length] = colnums[j]; + else + break; + // now #rowstart# is + // the number of entries in + // this line + + // for square matrices, the + // first entry in each row + // is the diagonal one. In + // this case only sort the + // remaining entries, otherwise + // sort all + sort ((rows==cols) ? &tmp_entries[1] : &tmp_entries[0], + &tmp_entries[row_length]); + + // Re-insert column numbers + // into the field + for (unsigned int j=0; j<row_length; ++j) + colnums[next_free_entry++] = tmp_entries[j]; + + // note new start of this and + // the next row + rowstart[line] = next_row_start; + next_row_start = next_free_entry; + + // some internal checks + Assert ((rows!=cols) || + (colnums[rowstart[line]] == static_cast<signed int>(line)), + ExcInternalError()); + // assert that the first entry + // does not show up in + // the remaining ones and that + // the remaining ones are unique + // among themselves (this handles + // both cases, quadratic and + // rectangular matrices) + Assert (find (&colnums[rowstart[line]+1], + &colnums[next_row_start], + colnums[rowstart[line]]) == + &colnums[next_row_start], + ExcInternalError()); + Assert (adjacent_find(&colnums[rowstart[line]+1], + &colnums[next_row_start]) == + &colnums[next_row_start], + ExcInternalError()); + }; + + vec_len = rowstart[rows] = next_row_start; compressed = true; - delete[] entries; -} + delete[] tmp_entries; +}; + int @@ -236,10 +165,31 @@ dSMatrixStruct::operator () (const unsigned int i, const unsigned int j) const { Assert (i<rows, ExcInvalidIndex(i,rows)); Assert (j<cols, ExcInvalidIndex(j,cols)); + Assert (compressed, ExcNotCompressed()); - for (unsigned int k=rowstart[i] ; k<rowstart[i+1] ; k++) - if (colnums[k] == (signed int)j) return k; - return -1; + // check first entry separately, since + // for square matrices this is + // the diagonal entry (check only + // if a first entry exists) + if (rowstart[i] != rowstart[i+1]) + { + if (static_cast<signed int>(j) == colnums[rowstart[i]]) + return rowstart[i]; + } + else + // no first entry exists for this + // line + return -1; + + // all other entries are sorted, so + // we can use a binary seach algorithm + const int* p = lower_bound (&colnums[rowstart[i]+1], + &colnums[rowstart[i+1]], + static_cast<signed int>(j)); + if (*p == static_cast<signed int>(j)) + return (p - &colnums[0]); + else + return -1; } @@ -248,6 +198,7 @@ dSMatrixStruct::add (const unsigned int i, const unsigned int j) { Assert (i<rows, ExcInvalidIndex(i,rows)); Assert (j<cols, ExcInvalidIndex(j,cols)); + Assert (compressed==false, ExcMatrixIsCompressed()); for (unsigned int k=rowstart[i]; k<rowstart[i+1]; k++) { @@ -273,9 +224,8 @@ dSMatrixStruct::add (const unsigned int i, const unsigned int j) void dSMatrixStruct::add_matrix (const unsigned int n, const int* rowcols) { - unsigned int i,j; - for (i=0;i<n;i++) - for (j=0;j<n;j++) + for (unsigned int i=0; i<n; ++i) + for (unsigned int j=0; j<n; ++j) add(rowcols[i], rowcols[j]); } @@ -285,9 +235,8 @@ void dSMatrixStruct::add_matrix (const unsigned int m, const unsigned int n, const int* rows, const int* cols) { - unsigned int i,j; - for (i=0;i<m;i++) - for (j=0;j<n;j++) + for (unsigned i=0; i<m; ++i) + for (unsigned j=0; j<n; ++j) add(rows[i], cols[j]); } @@ -636,6 +585,7 @@ dSMatrix::SSOR (dVector& dst, const double om) const dSMatrixStruct & dSMatrix::get_sparsity_pattern () const { + Assert (cols != 0, ExcMatrixNotInitialized()); return *cols; }; @@ -652,6 +602,7 @@ void dSMatrix::print (ostream &out) const { void dSMatrix::print_formatted (ostream &out, const unsigned int precision) const { + Assert (cols != 0, ExcMatrixNotInitialized()); out.precision (precision); out.setf (ios::scientific, ios::floatfield); // set output format