AX=B Fortran Petsc Code

Tim Stitt timothy.stitt at
Wed Nov 14 08:13:28 CST 2007

Dear PETSc Users/Developers,

I have the following sequential Fortran PETSc code that I have been 
developing (on and off) based on the kind advice given by members of 
this list, with respect to solving an inverse sparse matrix problem. 
Essentially, the code reads in a square double complex matrix from 
external file of size (order x order) and then proceeds to do a 
MatMatSolve(), where A is the sparse matrix to invert, B is a dense 
identity matrix and X is the resultant dense matrix....hope that makes 

My main problem is that the code stalls on the MatSetValues() for the 
sparse matrix A. With a trivial test matrix of (224 x 224) the program 
terminates successfully (by successfully I mean all instructions 
execute...I am not interested in the validity of X right now). 
Unfortunately, when I move up to a (2352 x 2352) matrix the 
MatSetValues() routine for matrix A is still in progress after 15 
minutes on one processor of our AMD Opteron IBM Cluster. I know that 
people will be screaming "preallocation"...but I have tried to take this 
into account by running a loop over the rows in A and counting the 
non-zero values explicitly prior to creation. I then pass this vector 
into the creation routine for the nnz argument. For the large (2352 x 
2352) problem that seems to be taking forever to most there are 
only 200 elements per row that are non-zero according to the counts.

Can anyone explain why the MatSetValues() routine is taking such a long 
time. Maybe this expected for this specific task...although it seems 
very long?

I did notice that on the trivial (224 x 224) run that I was still 
getting mallocs (approx 2000) for the A assembly when I used the -info 
command line parameter. I thought that it should be 0 if my 
preallocation counts were exact? Does this hint that I am doing 
something wrong. I have checked the code but don't see any obvious 
problems in the logic...not that means anything.

I would be grateful if someone could advise on this matter. Also, if you 
have a few seconds to spare I would be grateful if some experts could 
scan the remaining logic of the code (not in fine detail)  to make sure 
that I am doing all that I need to do to get this calculation 
working...assuming I can resolve the MatSetValues() problem.

Once again many thanks in advance,


 ! Initialise the PETSc MPI Harness
  call PetscInitialize(PETSC_NULL_CHARACTER,error);CHKERRQ(error)
  call MPI_COMM_SIZE(PETSC_COMM_SELF,processes,error);CHKERRQ(error)

  ! Read in Matrix
  read(321) order
  if (ID==0) then
     print *
     print *,processes," Processing Elements being used"
     print *
     print *,"Matrix has order ",order," rows by ",order," columns"
     print *
  end if

  read(321) matrix

  ! Allocate array for nnz

  ! Count number of non-zero elements in each matrix row
  do row=1,order
     do column=1,order
        if (matrix(row,column).ne.(0,0)) count=count+1
     end do
  end do

  ! Declare a PETSc Matrices


  ! Set up zero-based array indexing for use in MatSetValues

  do column=1,order
  end do

  ! Need to transpose values array as row-major arrays are used.

  ! Assemble Matrix A
  call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,error);CHKERRQ(error)
  call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,error);CHKERRQ(error)


  ! Create Index Sets for Factorisation
  call MatFactorInfoInitialize(info,error);CHKERRQ(error)
  call ISSetPermutation(indexSet,error);CHKERRQ(error)
  call MatLUFactorNumeric(A,info,factorMat,error);CHKERRQ(error)

  ! A no-longer needed
  call MatDestroy(A,error);CHKERRQ(error)


  ! Set Diagonal elements in Identity Matrix B
  do row=0,order-1
     call MatSetValue(B,row,row,one,INSERT_VALUES,error);CHKERRQ(error)
  end do

  ! Assemble B
  call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,error);CHKERRQ(error)
  call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,error);CHKERRQ(error)

  ! Assemble X
  call MatAssemblyBegin(X,MAT_FINAL_ASSEMBLY,error);CHKERRQ(error)
  call MatAssemblyEnd(X,MAT_FINAL_ASSEMBLY,error);CHKERRQ(error)

  ! Solve AX=B
  call MatMatSolve(factorMat,B,X,error);CHKERRQ(error)

  ! Deallocate Storage

  call MatDestroy(factorMat,error);CHKERRQ(error)
  call MatDestroy(B,error);CHKERRQ(error)
  call MatDestroy(X,error);CHKERRQ(error)

  call PetscFinalize(error)

Dr. Timothy Stitt <>
HPC Application Consultant - ICHEC (

Dublin Institute for Advanced Studies
5 Merrion Square - Dublin 2 - Ireland

+353-1-6621333 (tel) / +353-1-6621477 (fax)

More information about the petsc-users mailing list