Sparse Matrix Inversion using PETSc

[Dr. Timothy Stitt]

Perfect...thanks Barry.

Barry Smith wrote:
>   Tim,
>
>    Ok, this is a bit more reasonable :-)
>
>    The proceedure I discussed previously is the same. Just
> use for the right hand side matrix a dense vector with
> M columns whose entries are the first M columns of the identity,
> then solve with it.
>
>    Barry
>
>
> On Thu, 16 Aug 2007, Dr. Timothy Stitt wrote:
>
>   
>> Barry,
>>
>> If it helps I was speaking to some of the project members and they mention
>> that they actually only need the first M columns of the inverse. The equations
>> can also be rewritten so that they only require the last M columns also or
>> first/last M rows. I believe the retarded Green's Function forms an integral.
>> For the integral in the real plane (which is repeated for thousands of energy
>> points and hence requires thousands of inversions) they need only the first
>> few columns as described above.
>>
>> I would be grateful if you could suggest a possible approach based on this
>> extra information. I much appreciate your comments already.
>>
>> Thanks,
>>
>> Tim.
>>
>> Barry Smith wrote:
>>     
>>>   Tim,
>>>
>>>    Do you know what the G_r is then used for?
>>>
>>>    Thanks
>>>
>>>      Barry
>>>
>>>
>>> On Wed, 15 Aug 2007, Dr. Timothy Stitt wrote:
>>>
>>>   
>>>       
>>>> Barry,
>>>>
>>>> The group I am working with are calculating what they call retarded
>>>> Green's
>>>> Functions of the form:
>>>>
>>>> G_r=(E-H)^(-1)
>>>>
>>>> where (E-H) is a matrix. Apparently they say there is no way to avoid this
>>>> calculation.
>>>>
>>>> Tim.
>>>>
>>>> Barry Smith wrote:
>>>>     
>>>>         
>>>>>   Tim,
>>>>>
>>>>>     A dense matrix with 100,000 rows and columns requires 80 gigabytes
>>>>> to
>>>>> store
>>>>> the result. 1,000,000 rows and columns requires 8,000 gigabytes.    With
>>>>> this range of sizes I wouldn't even consider iterative solvers and
>>>>> would only use direct solvers. I would divide MPI_COMM_WORLD into N
>>>>> subcommunicators of size n each and have each subcommunicator work on a
>>>>> collection of columns of the inverse. For matrix sizes of 5,000 to say
>>>>> 20,000?
>>>>> I'd make n be 1 and just use PETSc's native LU solver and use
>>>>> MatMatSolve()
>>>>> and not use KSP at all. For larger matrices you may be able to use an n
>>>>> of 2
>>>>> to possibly as large as 8?
>>>>>
>>>>>    Now my numerical analysis training :-) requires me to state the
>>>>> following.
>>>>> It is completely insane to compute the EXPLICIT inverse of large sparse
>>>>> matrices
>>>>> since they are dense. Please tell me what the inverses are used for and
>>>>> perhaps
>>>>> we can come up with an approach the does not require computing them.
>>>>>
>>>>>    Barry
>>>>>
>>>>>
>>>>>
>>>>> On Wed, 15 Aug 2007, Dr. Timothy Stitt wrote:
>>>>>
>>>>>         
>>>>>           
>>>>>> Firstly, many thanks to everyone who has replied with information. It
>>>>>> has
>>>>>> been
>>>>>> very useful indeed. Much appreciated.
>>>>>>
>>>>>> Barry, in this case the sparse matrices would be of ~ order 5000x5000.
>>>>>> They
>>>>>> could grow in size but this is the sample matrices I am working with
>>>>>> right
>>>>>> now. We would love a scalable approach so we can deal with more
>>>>>> interesting
>>>>>> problems and hence larger sparse matrices. Hope that helps.
>>>>>>
>>>>>> Many thanks again.
>>>>>>
>>>>>> Tim.
>>>>>>
>>>>>> Barry Smith wrote:
>>>>>>             
>>>>>>             
>>>>>>>   Tim,
>>>>>>>
>>>>>>>     How large are you matrices?
>>>>>>>
>>>>>>>    Barry
>>>>>>>
>>>>>>>
>>>>>>> On Wed, 15 Aug 2007, Dr. Timothy Stitt wrote:
>>>>>>>
>>>>>>>                   
>>>>>>>               
>>>>>>>> Hi all,
>>>>>>>>
>>>>>>>> I am currently investigating the best way to perform the inversion
>>>>>>>> of
>>>>>>>> a
>>>>>>>> large
>>>>>>>> sparse matrix and came upon the idea of using PETSc as a framework
>>>>>>>> for
>>>>>>>> testing
>>>>>>>> various strategies from direct to iterative methods on my sample
>>>>>>>> matrices.
>>>>>>>> In
>>>>>>>> this setup for an NxN sparse matrix A I would have N rhs's
>>>>>>>> representing
>>>>>>>> the
>>>>>>>> Identity matrix and then solve for X. I wanted to experiment with
>>>>>>>> both
>>>>>>>> parallel and serial strategies ranging from LU Decomposition using
>>>>>>>> SuperLU,
>>>>>>>> MUMPS etc. to iterative methods using GMRES etc. Am I right in
>>>>>>>> thinking
>>>>>>>> that
>>>>>>>> all this can be done in PETSc by setting up a core framework and
>>>>>>>> then
>>>>>>>> varying
>>>>>>>> the solver methods etc?
>>>>>>>>
>>>>>>>> I have looked over the sample KSP Solver codes although they only
>>>>>>>> seem
>>>>>>>> to
>>>>>>>> suggest single vectors for x and b. Can this be changed to accept
>>>>>>>> multiple
>>>>>>>> vectors? Can anyone suggest a sample code that maybe demonstrates
>>>>>>>> the
>>>>>>>> sort
>>>>>>>> of
>>>>>>>> thing I want to achieve...if it is in fact possible.
>>>>>>>>
>>>>>>>> Thanks in advance for any assistance given,
>>>>>>>>
>>>>>>>> Regards,
>>>>>>>>
>>>>>>>> Tim.
>>>>>>>>
>>>>>>>>
>>>>>>>>                         
>>>>>>>>                 
>>>>>>>                   
>>>>>>>               
>>>>>>             
>>>>>>             
>>>>>         
>>>>>           
>>>>     
>>>>         
>>>   
>>>       
>>
>>     
>
>   


-- 
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) / +353-874195427 (mobile)