[petsc-users] Matrix dot product

Mon Sep 18 19:31:56 CDT 2017

On Mon, Sep 18, 2017 at 6:16 PM, David Gross <davegwebb10 at gmail.com> wrote:

> Hi Matt,
> I took a deeper look into the dev documentation and the src/mat code.
>

All the source is online, so I think that is a good way to talk about it in
email.

> From what I understand there are multiple places that the code needs to be
> added or referenced in. I created a copy of AXPY in axpy.c
>

Lets look at MatAXPY():

https://bitbucket.org/petsc/petsc/src/a2d7f5eaef7c3b636ab943ac50d4f7c86905df62/src/mat/utils/axpy.c?at=master&fileviewer=file-view-default#axpy.c-22

You just need to make a copy of this function. You want Z = X * Y, right?,
(or maybe Y = X * Y) so its

  MatPointwiseMult(Mat Z,  Mat X, Mat Y)

After the initial checks for compatibility, it dispatches

https://bitbucket.org/petsc/petsc/src/a2d7f5eaef7c3b636ab943ac50d4f7c86905df62/src/mat/utils/axpy.c?at=master&fileviewer=file-view-default#axpy.c-36

Here lets ignore the version optimized for a particular storage format, and
just use the generic version

https://bitbucket.org/petsc/petsc/src/a2d7f5eaef7c3b636ab943ac50d4f7c86905df62/src/mat/utils/axpy.c?at=master&fileviewer=file-view-default#axpy.c-58

Here we take advantage of the built-in ADD_VALUES, but there is no
MULT_VALUES option. Thus, we would need to
get the row of X and Y, compare columns (which come out sorted) and then
insert the result. However, the alternative is
have it error if MatStructure is DIFFERENT_NONZERO_PATTERN, so you can
assume that both rows are identical and
dispense with the column checking.

> and compiled without issue and could even reference it in my code without
> compiler error, but found that it wasn't actually doing anything as best I
> can tell. Am I correct in that one must make implementations for each
> matrix type in mat/impls? Also I assume that it needs to be added to
> matipl.h, include/petscmat.h, and  finclude/petscmat.h .
>
> Is there somewhere that documents the process of adding functions to PETSc
> classes? I couldn't find anything in the developers manual.
>

There is a discussion of this in the tutorials, and in here
http://www.mcs.anl.gov/petsc/developers/index.html there is a discussion of
adding different types of things.

  THanks,

    Matt

> Regards,
> Dave
>
> On Thu, Sep 14, 2017 at 10:07 PM, David Gross <davegwebb10 at gmail.com>
> wrote:
>
>> Hi Matt,
>> Noted for the naming convention.
>>
>> Yes, it is a Newton method (see Pan, V. and Reif, J., Fast and Efficient
>> Parallel Solution of Dense Linear Systems, Computers Math. Applications.
>> Vol. 17, No. 11, pp. 1481-1491, 1989)
>>
>> The dense matrix I have is repeatedly inverted while slowly changing such
>> that the previous inverse is a near perfect guess for the new inverse.
>>
>> Regards,
>> Dave
>>
>> On Thu, Sep 14, 2017 at 2:49 AM, Matthew Knepley <knepley at gmail.com>
>> wrote:
>>
>>> On Wed, Sep 13, 2017 at 6:14 PM, David Gross <davegwebb10 at gmail.com>
>>> wrote:
>>>
>>>> Hi Matt,
>>>> Thank you for getting back to me. Your answer confirms what I thought
>>>> in terms of existing functionality. I think it shouldn't be too hard to
>>>> make a copy of MatAXPY to MatAXY where it performs Xij = A*Xij*Yij (or
>>>> without the A). I could then do the MatNorm of the resulting matrix to get
>>>> what I need.
>>>>
>>>> Is a MatAXY function desirable as a source contribution?
>>>>
>>>
>>> Yes. I would prefer calling it MatPointwiseMult, since you can see it as
>>> VecPointwiseMult on a Vec obtained
>>> from forgetting the linear operator structure of the matrix (forgetful
>>> functor).
>>>
>>>
>>>> I am hoping to use PETSc for performing basic vector and matrix
>>>> operations on dense matrices and 1D vectors. The main uses are matmult,
>>>> matmatmult and matrix additions and scaling. The application is for
>>>> implementing a parallel version on an existing Pan-Reif matrix inversion
>>>> algorithm.
>>>>
>>>
>>> Is this Newton's method on the vector space of matrices?
>>>
>>>   Thanks,
>>>
>>>     Matt
>>>
>>>
>>>> The choice of using PETSc is mostly due to us already using it in the
>>>> same program to solve sparse matrices (with MUMPS) with the goal of
>>>> avoiding adding yet another package (ex ScaLAPACK/PBLAS) into the code even
>>>> if other packages may be more directly oriented towards my application.
>>>>
>>>> Regards,
>>>> Dave
>>>>
>>>>
>>>>
>>>> On Mon, Sep 11, 2017 at 2:19 AM, Matthew Knepley <knepley at gmail.com>
>>>> wrote:
>>>>
>>>>> On Sun, Sep 10, 2017 at 5:51 PM, David Gross <davegwebb10 at gmail.com>
>>>>> wrote:
>>>>>
>>>>>> Hello,
>>>>>> I was wondering if there was a matrix equivalent to the vecDot
>>>>>> function (Frobenius inner product)? As far as I can tell the closest thing
>>>>>> is MatNorm with NORM_FROBENIUS, but obviously this is acting on only one
>>>>>> matrix.
>>>>>>
>>>>>> If there is not a built in function, what is the best way to compute
>>>>>> this? I am working fortran90.
>>>>>>
>>>>>
>>>>> We do not have this. However, it would be trivial to add since we have
>>>>>
>>>>>   http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpage
>>>>> s/Mat/MatAXPY.html
>>>>>
>>>>> since you just replace + with * in our code. You could argue that we
>>>>> should have written for
>>>>> a general ring, but C makes this cumbersome. Do you think you could
>>>>> make the change?
>>>>>
>>>>> What are you using this for?
>>>>>
>>>>>   Thanks,
>>>>>
>>>>>      Matt
>>>>>
>>>>>
>>>>>> Regards,
>>>>>> Dave
>>>>>>
>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> What most experimenters take for granted before they begin their
>>>>> experiments is infinitely more interesting than any results to which their
>>>>> experiments lead.
>>>>> -- Norbert Wiener
>>>>>
>>>>> http://www.caam.rice.edu/~mk51/
>>>>>
>>>>
>>>>
>>>
>>>
>>> --
>>> What most experimenters take for granted before they begin their
>>> experiments is infinitely more interesting than any results to which their
>>> experiments lead.
>>> -- Norbert Wiener
>>>
>>> http://www.caam.rice.edu/~mk51/
>>>
>>
>>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

http://www.caam.rice.edu/~mk51/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170918/b89b02e6/attachment-0001.html>