[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

AX=B Fortran Petsc Code

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)