[petsc-users] PETSc set off-diagonal matrix pre-allocation
Barry Smith
bsmith at mcs.anl.gov
Tue Feb 14 20:05:59 CST 2017
1) on a single process all the columns are in the diagonal block and none in the off diagonal block (note that the column partitioning corresponding to a partitioning of the vector in the product A *x
2) on two processes you guessed wrong how many columns PETSc would put on the first process. You guessed it would put the first three on the first process and the last three on the second process.
a) it cannot do that, every column has to been owned by some process (in the vector x above) so it cannot be 3 and 3, it has to be 4 and 3 or 3 and 4.
b) PETSc puts the "extra" columns on the earlier processes not the later processes.
For rectangular matrices you really cannot get away with using "PETSC_DECIDE" for local columns when using preallocation
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ex3.c
Type: application/octet-stream
Size: 1678 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170214/dbdf7c02/attachment.obj>
-------------- next part --------------
, since it may decide something different than what you assume.
I've attached the parallel code that behaves as expected.
Barry
> On Feb 14, 2017, at 1:06 PM, Andrew Ho <andrewh0 at uw.edu> wrote:
>
> The problem isn't only with 1 process, but it seems to be with all non-square matrices?
>
> For example, here I have a two process program which tries to set a 6x7 matrix to (each process owns 3 rows):
>
> 0 0 0 1 1 1 1
> 0 0 0 1 1 1 1
> 0 0 0 1 1 1 1
> 0 0 0 1 1 1 1
> 0 0 0 1 1 1 1
> 0 0 0 1 1 1 1
>
> #include <petsc.h>
> #include <mpi.h>
>
> int main(int argc, char** argv)
> {
> PetscErrorCode err;
> err = PetscInitialize(&argc, &argv, NULL, "help");
> CHKERRQ(err);
>
> int rank, size;
> MPI_Comm_rank(MPI_COMM_WORLD, &rank);
> MPI_Comm_size(MPI_COMM_WORLD, &size);
> if(size != 2)
> {
> printf("must run with 2 processes");
> MPI_Abort(MPI_COMM_WORLD, -1);
> }
>
> // create a sparse AIJ matrix distributed across MPI
> Mat A;
> err = MatCreate(MPI_COMM_WORLD, &A);
> CHKERRQ(err);
> err = MatSetType(A, MATMPIAIJ);
> CHKERRQ(err);
> // setup pre-allocation for matrix space
> {
> err =
> MatSetSizes(A, 3, PETSC_DECIDE, 6, 7);
> CHKERRQ(err);
> if(rank == 0)
> {
> PetscInt d_nnz[] = {0, 0, 0};
> PetscInt o_nnz[] = {4, 4, 4};
> err = MatMPIAIJSetPreallocation(A, 0, d_nnz, 0, o_nnz);
> CHKERRQ(err);
> }
> else
> {
> PetscInt d_nnz[] = {3, 3, 3};
> PetscInt o_nnz[] = {1, 1, 1};
> err = MatMPIAIJSetPreallocation(A, 0, d_nnz, 0, o_nnz);
> CHKERRQ(err);
> }
> }
> err = MatSetUp(A);
> CHKERRQ(err);
>
> // set values inside the matrix
> for (PetscInt row = 0; row < 3; ++row)
> {
> for (PetscInt col = 3; col < 7; ++col)
> {
> err = MatSetValue(A, 3 * rank + row, col, 1, INSERT_VALUES);
> CHKERRQ(err);
> }
> }
>
> err = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
> CHKERRQ(err);
> err = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
> CHKERRQ(err);
>
> err = MatView(A, PETSC_VIEWER_STDOUT_WORLD);
> CHKERRQ(err);
>
> // free memory
> err = MatDestroy(&A);
> CHKERRQ(err);
>
> // cleanup any internal PETSc data at end of program
> err = PetscFinalize();
> CHKERRQ(err);
> }
>
>
> On Tue, Feb 14, 2017 at 5:26 AM, Matthew Knepley <knepley at gmail.com> wrote:
> On Tue, Feb 14, 2017 at 6:41 AM, Andrew Ho <andrewh0 at uw.edu> wrote:
> I have a 3x6 matrix, and I want to set it to (just as an example):
>
> 0 0 0 1 1 1
> 0 0 0 1 1 1
> 0 0 0 1 1 1
>
> From what I can tell, according to the documentation the MPIAIJ sparsity of this matrix is:
>
> I believe this is an inconsistency in PETSc for rectangular matrices. We divide matrices into diagonal and
> off-diagonal parts mainly to separate communication from computation. Thus on 1 proc, we never divide
> them, even if the matrix is rectangular. Therefore we misinterpret your preallocation.
>
> Barry, do you think we want to try and change this?
>
> Matt
>
> d_nnz = [0, 0, 0]
> o_nnz = [3, 3, 3]
>
> However, when I do this, I get the following error:
>
> [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
> [0]PETSC ERROR: Argument out of range
> [0]PETSC ERROR: New nonzero at (0,3) caused a malloc
> Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off this check
> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.7.5, unknown
> [0]PETSC ERROR: ./ex3 on a arch-linux2-c-opt Tue Feb 14 04:23:56 2017
> [0]PETSC ERROR: Configure options --with-debugging=0 --COPTFLAGS="-O3 -march=native" --CXXOPTFLAGS="-O3 -march=native" --FOPTFLAGS="-O3 -march=native"
> [0]PETSC ERROR: #1 MatSetValues_MPIAIJ() line 582 in petsc/src/mat/impls/aij/mpi/mpiaij.c
> [0]PETSC ERROR: #2 MatSetValues() line 1190 in petsc/src/mat/interface/matrix.c
> [0]PETSC ERROR: #3 main() line 36 in ex3.c
> [0]PETSC ERROR: No PETSc Option Table entries
> [0]PETSC ERROR: ----------------End of Error Message -------send entire error message to petsc-maint at mcs.anl.gov----------
>
> Here's a working test code:
>
> #include <petsc.h>
> #include <mpi.h>
> int main(int argc, char** argv)
> {
> PetscErrorCode err;
> err = PetscInitialize(&argc, &argv, NULL, "help");
> CHKERRQ(err);
> // create a sparse AIJ matrix distributed across MPI
> PetscInt global_width = 6;
> PetscInt global_height = 3;
> Mat A;
> err = MatCreate(MPI_COMM_WORLD, &A);
> CHKERRQ(err);
> err = MatSetType(A, MATMPIAIJ);
> CHKERRQ(err);
> // setup pre-allocation for matrix space
> {
> err =
> MatSetSizes(A, global_height, PETSC_DECIDE, global_height, global_width);
> CHKERRQ(err);
> PetscInt d_nnz[] = {0, 0, 0};
> PetscInt o_nnz[] = {3, 3, 3};
> err = MatMPIAIJSetPreallocation(A, 0, d_nnz, 0, o_nnz);
> CHKERRQ(err);
> }
> err = MatSetUp(A);
> CHKERRQ(err);
> // set values inside the matrix
> for (PetscInt row = 0; row < global_height; ++row)
> {
> for (PetscInt col = global_height; col < global_width; ++col)
> {
> err = MatSetValue(A, row, col, 1, INSERT_VALUES);
> CHKERRQ(err);
> }
> }
> err = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
> CHKERRQ(err);
> err = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
> CHKERRQ(err);
> err = MatView(A, PETSC_VIEWER_STDOUT_WORLD);
> CHKERRQ(err);
> // free memory
> err = MatDestroy(&A);
> CHKERRQ(err);
> // cleanup any internal PETSc data at end of program
> err = PetscFinalize();
> CHKERRQ(err);
> }
>
> Am I mis-understanding what the d_nnz and o_nnz parameter are supposed to mean?
> --
> Andrew Ho
>
>
>
> --
> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
> -- Norbert Wiener
>
>
>
> --
> Andrew Ho
More information about the petsc-users
mailing list