AX=B Fortran Petsc Code

[Tim Stitt]

Thanks Barry.

Barry Smith wrote:
>
>   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)
>>
>


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