AX=B Fortran Petsc Code
Tim Stitt
timothy.stitt at ichec.ie
Wed Nov 14 10:37:43 CST 2007
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?
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)
More information about the petsc-users
mailing list