AX=B Fortran Petsc Code
Barry Smith
bsmith at mcs.anl.gov
Thu Nov 15 13:43:50 CST 2007
Tim,
There is an field in MatFactorInfo that contains this fill factor
called
fill set it with 2.5 and you should be all set.
Barry
On Nov 14, 2007, at 12:04 PM, Tim Stitt wrote:
> 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