[petsc-dev] Number of nonzeros after assembly?

[Gabriele Jost]

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