PetscMatTools.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 #include "PetscMatTools.hpp"
00030 #include <cassert>
00031 
00032 
00034 // Implementation
00036 
00037 void PetscMatTools::SetElement(Mat matrix, PetscInt row, PetscInt col, double value)
00038 {
00039     PetscInt lo, hi;
00040     GetOwnershipRange(matrix, lo, hi);
00041 
00042     if (row >= lo && row < hi)
00043     {
00044         MatSetValue(matrix, row, col, value, INSERT_VALUES);
00045     }
00046 }
00047 
00048 void PetscMatTools::AddToElement(Mat matrix, PetscInt row, PetscInt col, double value)
00049 {
00050     PetscInt lo, hi;
00051     GetOwnershipRange(matrix, lo, hi);
00052 
00053     if (row >= lo && row < hi)
00054     {
00055         MatSetValue(matrix, row, col, value, ADD_VALUES);
00056     }
00057 }
00058 
00059 void PetscMatTools::AssembleFinal(Mat matrix)
00060 {
00061     MatAssemblyBegin(matrix, MAT_FINAL_ASSEMBLY);
00062     MatAssemblyEnd(matrix, MAT_FINAL_ASSEMBLY);
00063 }
00064 
00065 void PetscMatTools::AssembleIntermediate(Mat matrix)
00066 {
00067     MatAssemblyBegin(matrix, MAT_FLUSH_ASSEMBLY);
00068     MatAssemblyEnd(matrix, MAT_FLUSH_ASSEMBLY);
00069 }
00070 
00071 void PetscMatTools::Display(Mat matrix)
00072 {
00073     MatView(matrix,PETSC_VIEWER_STDOUT_WORLD);
00074 }
00075 
00076 void PetscMatTools::SetRow(Mat matrix, PetscInt row, double value)
00077 {
00078     PetscInt lo, hi;
00079     GetOwnershipRange(matrix, lo, hi);
00080 
00081     if (row >= lo && row < hi)
00082     {
00083         PetscInt rows, cols;
00084         MatGetSize(matrix, &rows, &cols);
00085         for (PetscInt i=0; i<cols; i++)
00086         {
00087             SetElement(matrix, row, i, value);
00088         }
00089     }
00090 }
00091 
00092 void PetscMatTools::ZeroRowsWithValueOnDiagonal(Mat matrix, std::vector<unsigned>& rRows, double diagonalValue)
00093 {
00094     AssembleFinal(matrix);
00095 
00096     // Important! Petsc by default will destroy the sparsity structure for this row and deallocate memory
00097     // when the row is zeroed, and if there is a next timestep, the memory will have to reallocated
00098     // when assembly to done again. This can kill performance. The following makes sure the zeroed rows
00099     // are kept.
00100 #if (PETSC_VERSION_MAJOR == 3) //PETSc 3.x.x
00101  #if (PETSC_VERSION_MINOR == 0)
00102     MatSetOption(matrix, MAT_KEEP_ZEROED_ROWS, PETSC_TRUE);
00103  #else
00104     MatSetOption(matrix, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
00105  #endif
00106 #else
00107     MatSetOption(matrix, MAT_KEEP_ZEROED_ROWS);
00108 #endif
00109 
00110     PetscInt* rows = new PetscInt[rRows.size()];
00111     for(unsigned i=0; i<rRows.size(); i++)
00112     {
00113         rows[i] = rRows[i];
00114     }
00115 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
00116     IS is;
00117     ISCreateGeneral(PETSC_COMM_WORLD, rRows.size(), rows, &is);
00118     MatZeroRows(matrix, is, &diagonalValue);
00119     ISDestroy(is);
00120     /*
00121      * 
00122 [2]PETSC ERROR: MatMissingDiagonal_SeqAIJ() line 1011 in src/mat/impls/aij/seq/aij.c
00123 [2]PETSC ERROR: Petsc has generated inconsistent data!
00124 [2]PETSC ERROR: Matrix is missing diagonal number 15!
00125 [2]PETSC ERROR: MatILUFactorSymbolic_SeqAIJ() line 906 in src/mat/impls/aij/seq/aijfact.c
00126      * 
00127      * 
00128     //There appears to be a problem with MatZeroRows not setting diagonals correctly
00129     //While we are supporting PETSc 2.2, we have to do this the slow way
00130 
00131     AssembleFinal(matrix);
00132     PetscInt lo, hi;
00133     GetOwnershipRange(matrix, lo, hi);
00134     PetscInt size=GetSize(matrix);
00136     for(unsigned index=0; index<rRows.size(); index++)
00137     {
00138         PetscInt row = rRows[index];
00139         if (row >= lo && row < hi)
00140         {
00141             std::vector<unsigned> non_zero_cols;
00142             //This row is local, so zero it.
00143             for (PetscInt column = 0; column < size; column++)
00144             {
00145                 if (GetElement(matrix, row, column) != 0.0)
00146                 {
00147                     non_zero_cols.push_back(column);
00148                 }
00149             }
00150             //Separate "gets" from "sets" so that we don't have to keep going into "assembled" mode
00151             for (unsigned i=0; i<non_zero_cols.size();i++)
00152             {
00153                 SetElement(matrix, row, non_zero_cols[i], 0.0);
00154             }
00155             //Set the diagonal
00156             SetElement(matrix, row, row, diagonalValue);
00157         }
00158         //Everyone communicate after row is finished
00159         AssembleFinal(matrix);
00160     }
00161     */     
00162 #else
00163     MatZeroRows(matrix, rRows.size(), rows, diagonalValue);
00164 #endif
00165     delete [] rows;
00166 }
00167 
00168 
00169 void PetscMatTools::ZeroRowsAndColumnsWithValueOnDiagonal(Mat matrix, std::vector<unsigned>& rRowColIndices, double diagonalValue)
00170 {
00171     AssembleFinal(matrix);
00172 
00173     PetscInt lo, hi;
00174     GetOwnershipRange(matrix, lo, hi);
00175     std::vector<unsigned>* p_nonzero_rows_per_column = new std::vector<unsigned>[rRowColIndices.size()];
00176 
00177     // for each column: collect all the row indices corresponding to a non-zero entry
00178     // We do all the columns at once, before doing the zeroing, as otherwise
00179     // a MatAssemblyBegin() & MatAssemblyEnd() would have to be called
00180     // after every MatSetValues and before the below GetMatrixElement()
00181     for(unsigned index=0; index<rRowColIndices.size(); index++)
00182     {
00183         unsigned column = rRowColIndices[index];
00184 
00185         // determine which rows in this column are non-zero (and
00186         // therefore need to be zeroed)
00187         for (PetscInt row = lo; row < hi; row++)
00188         {
00189             if (GetElement(matrix, row, column) != 0.0)
00190             {
00191                 p_nonzero_rows_per_column[index].push_back(row);
00192             }
00193         }
00194     }
00195 
00196     // Now zero each column in turn
00197     for(unsigned index=0; index<rRowColIndices.size(); index++)
00198     {
00199         // set those rows to be zero by calling MatSetValues
00200         unsigned size = p_nonzero_rows_per_column[index].size();
00201         PetscInt* rows = new PetscInt[size];
00202         PetscInt cols[1];
00203         double* zeros = new double[size];
00204 
00205         cols[0] = rRowColIndices[index];
00206 
00207         for (unsigned i=0; i<size; i++)
00208         {
00209             rows[i] = p_nonzero_rows_per_column[index][i];
00210             zeros[i] = 0.0;
00211         }
00212 
00213         MatSetValues(matrix, size, rows, 1, cols, zeros, INSERT_VALUES);
00214         delete [] rows;
00215         delete [] zeros;
00216     }
00217     delete[] p_nonzero_rows_per_column;
00218 
00219     // now zero the rows and add the diagonal entries
00220     ZeroRowsWithValueOnDiagonal(matrix, rRowColIndices, diagonalValue);
00221 }
00222 
00223 void PetscMatTools::ZeroColumn(Mat matrix, PetscInt col)
00224 {
00225     AssembleFinal(matrix);
00226 
00227     PetscInt lo, hi;
00228     GetOwnershipRange(matrix, lo, hi);
00229 
00230     // determine which rows in this column are non-zero (and
00231     // therefore need to be zeroed)
00232     std::vector<unsigned> nonzero_rows;
00233     for (PetscInt row = lo; row < hi; row++)
00234     {
00235         if (GetElement(matrix, row, col) != 0.0)
00236         {
00237             nonzero_rows.push_back(row);
00238         }
00239     }
00240 
00241     // set those rows to be zero by calling MatSetValues
00242     unsigned size = nonzero_rows.size();
00243     PetscInt* rows = new PetscInt[size];
00244     PetscInt cols[1];
00245     double* zeros = new double[size];
00246 
00247     cols[0] = col;
00248 
00249     for (unsigned i=0; i<size; i++)
00250     {
00251         rows[i] = nonzero_rows[i];
00252         zeros[i] = 0.0;
00253     }
00254 
00255     MatSetValues(matrix, size, rows, 1, cols, zeros, INSERT_VALUES);
00256     delete [] rows;
00257     delete [] zeros;
00258 }
00259 
00260 void PetscMatTools::Zero(Mat matrix)
00261 {
00262     MatZeroEntries(matrix);
00263 }
00264 
00265 unsigned PetscMatTools::GetSize(Mat matrix) 
00266 {
00267     PetscInt rows, cols;
00268     
00269     MatGetSize(matrix, &rows, &cols);
00270     assert(rows == cols);
00271     return (unsigned) rows;
00272 }
00273 
00274 void PetscMatTools::GetOwnershipRange(Mat matrix, PetscInt& lo, PetscInt& hi)
00275 {
00276     MatGetOwnershipRange(matrix, &lo, &hi);
00277 }
00278 
00279 double PetscMatTools::GetElement(Mat matrix, PetscInt row, PetscInt col)
00280 {
00281     PetscInt lo, hi;
00282     GetOwnershipRange(matrix, lo, hi);
00283 
00284     assert(lo <= row && row < hi);
00285     PetscInt row_as_array[1];
00286     row_as_array[0] = row;
00287     PetscInt col_as_array[1];
00288     col_as_array[0] = col;
00289 
00290     double ret_array[1];
00291 
00292     MatGetValues(matrix, 1, row_as_array, 1, col_as_array, ret_array);
00293 
00294     return ret_array[0];
00295 }
00296 
00297 
00298 void PetscMatTools::SetOption(Mat matrix, MatOption option)
00299 {
00300 #if (PETSC_VERSION_MAJOR == 3) //PETSc 3.x.x
00301     MatSetOption(matrix, option, PETSC_TRUE);
00302 #else
00303     MatSetOption(matrix, option);
00304 #endif
00305 }
00306 
00307 
00308 
00309 Vec PetscMatTools::GetMatrixRowDistributed(Mat matrix, unsigned rowIndex)
00310 {
00311     //  We need to make sure that lhs_ith_row doesn't ignore off processor entries when assembling,
00312     //  otherwise the VecSetValues call a few lines below will not work as expected.
00313 
00314     PetscInt lo, hi;
00315     PetscMatTools::GetOwnershipRange(matrix, lo, hi);
00316     unsigned size = PetscMatTools::GetSize(matrix);
00317 
00318     Vec mat_ith_row = PetscTools::CreateVec(size, hi-lo, false);
00319 
00320     PetscInt num_entries;
00321     const PetscInt* column_indices;
00322     const PetscScalar* values;
00323 
00324     bool am_row_owner = (PetscInt)rowIndex >= lo && (PetscInt)rowIndex < hi;
00325 
00326     // Am I the owner of the row? If so get the non-zero entries and add them lhs_ith_row.
00327     // In parallel, VecAssembly{Begin,End} will send values to the rest of processors.
00328     if (am_row_owner)
00329     {
00330         MatGetRow(matrix, rowIndex, &num_entries, &column_indices, &values);
00331         VecSetValues(mat_ith_row, num_entries, column_indices, values, INSERT_VALUES);
00332     }
00333 
00334     VecAssemblyBegin(mat_ith_row);
00335     VecAssemblyEnd(mat_ith_row);
00336 
00337     if (am_row_owner)
00338     {
00339         MatRestoreRow(matrix, rowIndex, &num_entries, &column_indices, &values);
00340     }
00341 
00342     return mat_ith_row;
00343 }

Generated on Tue May 31 14:31:42 2011 for Chaste by  doxygen 1.5.5