[petsc-users] [petsc-dev] Number of nonzeros after assembly?

[Barry Smith]

>     for (col=start; col<end;  col++) {

   You are missing all the matrix entries that are not in the "block diagonal" portion of the matrix on each process. This loop needs to start at 0 and go to M-1

  Barry

   
On Aug 22, 2011, at 1:35 PM, Gabriele Jost wrote:

> Hello!
> I have a question about the number of nonzeros in a PETSc matrix (mpiaij format).
> The context is, that I am trying to debug a program which works fine for a certain type of problems, but gives incorrect results for a different type.
> Anyhow, the input matrix has  5222556 nonzeros, this is the smalles data set I can get from the user.  However, if I look at the KSP info I get from PETSc, I see the message:
> 
> KSP Object:
>  type: preonly
>  maximum iterations=5000, initial guess is zero
>  tolerances:  relative=1e-06, absolute=1e-50, divergence=10000
>  left preconditioning
> PC Object:
>  type: lu
>    LU: out-of-place factorization
>      matrix ordering: natural
>    LU: tolerance for zero pivot 1e-12
>    LU: factor fill ratio needed 0
>         Factored matrix follows
>        Matrix Object:
>          type=mpiaij, rows=13608, cols=13608
>          package used to perform factorization: superlu_dist
>          total: nonzeros=0, allocated nonzeros=27216
>            SuperLU_DIST run parameters:
>              Process grid nprow 4 x npcol 4
>              Equilibrate matrix TRUE
>              Matrix input mode 1
>              Replace tiny pivots TRUE
>              Use iterative refinement FALSE
>              Processors in row 4 col partition 4
>              Row permutation LargeDiag
>              Column permutation MMD_AT_PLUS_A
>              Parallel symbolic factorization FALSE
>              Repeated factorization SamePattern
>  linear system matrix = precond matrix:
>  Matrix Object:
>    type=mpiaij, rows=13608, cols=13608
>    total: nonzeros=5227411, allocated nonzeros=5308814
> 
> So, there are more non-zeros than in my matrix. This discrepancy does not show up in the good cases.
> So, I was trying to count the number of non-zeros after the assembly stage. I do this by:
> 
>  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>  if (ierr != 0)  SETERRQ(1,"assembly error");
>  {
>    int p,row,col,start,end,icnt;
>    icnt=0;
>    ierr=MatGetOwnershipRangeColumn(A, &start, &end);
>    for (row=Istarts[mytid]; row<Iends[mytid]; row++) {
>       for (col=start; col<end;  col++) {
>          PetscScalar v=0.;
>          ierr = MatGetValues(A,1,&row,1,&col,&v); CHKERRQ(ierr);
>          if (v!=0.) {
>                icnt++;
>                PetscSynchronizedPrintf(PETSC_COMM_WORLD,"%d %d %lf\n",row, col, v);
>          //PetscSynchronizedPrintf(PETSC_COMM_WORLD,"%d %d %lf\n",row, col,v);
>          }
>        }
>     }
>     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"%s %d %d\n","Processor ",mytid,icnt);
>     PetscSynchronizedFlush(PETSC_COMM_WORLD);
>   }
> 
> If I sum up the nonzeros for each processor, I get 591479... by far not enough. 
> I would greatly appreciate if somebody could advise me on what is going on here. Many thanks in advance for any hints!
> 
> Kind regards,
> Gabriele Jost