AX=B Fortran Petsc Code
Tim Stitt
timothy.stitt at ichec.ie
Wed Nov 14 12:04:21 CST 2007
OK...everything is working well now and I am getting the results I
expect. Much appreciated.
Saying that...I am trying to now satisfy the PC FACTOR FILL suggestion
provided by the -info parameter on my sample sparse matrices.
In my case I am getting:
[0] MatLUFactorSymbolic_SeqAIJ(): Reallocs 3 Fill ratio:given 0 needed
2.56568
[0] MatLUFactorSymbolic_SeqAIJ(): Run with -pc_factor_fill 2.56568 or use
[0] MatLUFactorSymbolic_SeqAIJ(): PCFactorSetFill(pc,2.56568);
[0] MatLUFactorSymbolic_SeqAIJ(): for best performance.
So I run my code with ./foo -pc_factor_fill 2.56568
but I continually get
WARNING! There are options you set that were not used!
WARNING! could be spelling mistake, etc!
Option left: name:-pc_factor_fill value: 2.56568
Can someone suggest how I can improve performance with the
pc_factor_fill parameter in my case? As my code stands there is no
MatSetFromOptions() as I set everything explicitly in the code.
Thanks again,
Tim.
Barry Smith wrote:
>
> For sequential codes the index set is as large as the matrix.
>
> In parallel the factor codes do not use the ordering, they do the
> ordering
> internally.
>
> Barry
>
> On Nov 14, 2007, at 10:37 AM, Tim Stitt 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?
>>
>> 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)
>>
>
--
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