AX=B Fortran Petsc Code

[Matthew Knepley]

On Nov 14, 2007 10:37 AM, Tim Stitt <timothy.stitt at ichec.ie> wrote:
> Can I just ask a question about MatLUFactorSymbolic() in this context?
> What sizes should the 'row' and 'col' index sets be? Should they span
> all global rows/columns in A?

Yes, the matrix is permuted as a whole.

   Matt

> Matthew Knepley wrote:
> > You appear to be setting every value in the sparse matrix. We do not
> > throw out 0 values (since sometimes they are necessary for structural
> > reasons). Thus you are allocating a ton of times. You need to remove
> > the 0 values before calling MatSetValues (and their associated
> > column entires as well).
> >
> >   Matt
> >
> > On Nov 14, 2007 8:13 AM, Tim Stitt <timothy.stitt at ichec.ie> wrote:
> >
> >> 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
> >> sense.
> >>
> >> 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 set...at 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,
> >>
> >> Tim.
> >>
> >>  ! Initialise the PETSc MPI Harness
> >>   call PetscInitialize(PETSC_NULL_CHARACTER,error);CHKERRQ(error)
> >>
> >>   call MPI_COMM_SIZE(PETSC_COMM_SELF,processes,error);CHKERRQ(error)
> >>   call MPI_COMM_RANK(PETSC_COMM_SELF,ID,error);CHKERRQ(error)
> >>
> >>   ! Read in Matrix
> >>   open(321,file='Hamiltonian.bin',form='unformatted')
> >>   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
> >>
> >>   allocate(matrix(order,order))
> >>   read(321) matrix
> >>   close(321)
> >>
> >>   ! Allocate array for nnz
> >>   allocate(numberZero(order))
> >>
> >>   ! Count number of non-zero elements in each matrix row
> >>   do row=1,order
> >>      count=0
> >>      do column=1,order
> >>         if (matrix(row,column).ne.(0,0)) count=count+1
> >>      end do
> >>      numberZero(row)=count
> >>   end do
> >>
> >>   ! Declare a PETSc Matrices
> >>
> >>   call
> >> MatCreateSeqAIJ(PETSC_COMM_SELF,order,order,PETSC_NULL_INTEGER,numberZero,A,error);CHKERRQ(error)
> >>   call
> >> MatCreateSeqAIJ(PETSC_COMM_SELF,order,order,0,PETSC_NULL_INTEGER,factorMat,error);CHKERRQ(error)
> >>   call
> >> MatCreateSeqDense(PETSC_COMM_SELF,order,order,PETSC_NULL_SCALAR,X,error);CHKERRQ(error)
> >>   call
> >> MatCreateSeqDense(PETSC_COMM_SELF,order,order,PETSC_NULL_SCALAR,B,error);CHKERRQ(error)
> >>
> >>   ! Set up zero-based array indexing for use in MatSetValues
> >>   allocate(columnIndices(order))
> >>
> >>   do column=1,order
> >>      columnIndices(column)=column-1
> >>   end do
> >>
> >>   ! Need to transpose values array as row-major arrays are used.
> >>   call
> >> MatSetValues(A,order,columnIndices,order,columnIndices,transpose(matrix),INSERT_VALUES,error);CHKERRQ(error)
> >>
> >>   ! Assemble Matrix A
> >>   call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,error);CHKERRQ(error)
> >>   call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,error);CHKERRQ(error)
> >>
> >>   deallocate(matrix)
> >>
> >>   ! Create Index Sets for Factorisation
> >>   call
> >> ISCreateGeneral(PETSC_COMM_SELF,order,columnIndices,indexSet,error);CHKERRQ(error)
> >>   call MatFactorInfoInitialize(info,error);CHKERRQ(error)
> >>   call ISSetPermutation(indexSet,error);CHKERRQ(error)
> >>   call
> >> MatLUFactorSymbolic(A,indexSet,indexSet,info,factorMat,error);CHKERRQ(error)
> >>   call MatLUFactorNumeric(A,info,factorMat,error);CHKERRQ(error)
> >>
> >>   ! A no-longer needed
> >>   call MatDestroy(A,error);CHKERRQ(error)
> >>
> >>   one=(1,0)
> >>
> >>   ! 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
> >>   deallocate(columnIndices)
> >>
> >>   call MatDestroy(factorMat,error);CHKERRQ(error)
> >>   call MatDestroy(B,error);CHKERRQ(error)
> >>   call MatDestroy(X,error);CHKERRQ(error)
> >>
> >>   call PetscFinalize(error)
> >>
> >> --
> >> Dr. Timothy Stitt <timothy_dot_stitt_at_ichec.ie>
> >> HPC Application Consultant - ICHEC (www.ichec.ie)
> >>
> >> Dublin Institute for Advanced Studies
> >> 5 Merrion Square - Dublin 2 - Ireland
> >>
> >> +353-1-6621333 (tel) / +353-1-6621477 (fax)
> >>
> >>
> >>
> >
> >
> >
> >
>
>
> --
> Dr. Timothy Stitt <timothy_dot_stitt_at_ichec.ie>
> HPC Application Consultant - ICHEC (www.ichec.ie)
>
> Dublin Institute for Advanced Studies
> 5 Merrion Square - Dublin 2 - Ireland
>
> +353-1-6621333 (tel) / +353-1-6621477 (fax)
>
>



-- 
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