[petsc-users] Getting access to matrix rows without and setting values simultaneously

[Alexander Grayver]

On 16.06.2011 16:07, Matthew Knepley wrote:
> On Thu, Jun 16, 2011 at 1:43 PM, Alexander Grayver 
> <agrayver at gfz-potsdam.de <mailto:agrayver at gfz-potsdam.de>> wrote:
>
>     >> Could you explain more about the task you're trying to do.
>
>     Well, I can try. That is pretty specific task, I wouldn't like to
>     go deep into details.
>
>     Let's say we have a discretized model of some physical parameter m
>     (say acoustic velocity). Number of model parameters is N. We need
>     to take M measurements d within model (say time traveling of
>     acoustic wave) based on M different sets of receiver/source positions.
>     In our case N >> M (e.g. 10^7 >> 10^3).
>     We have operator F (which might linear or not) that defines
>     relationship between m and d:
>     d=F(m)
>
>     The operator F is just a system of equations actually. I have no
>     problem now to solve it for any m.
>
>     What I need now is to compute the variation of this operator:
>
>     A = \frac{\partial d_i}{\partial m_j}, for all i=1..M, j=1..N
>
>     After some maths this calculation could be reformulated as a
>     triple product:
>
>     A_i = C*F(m)^-1*v
>
>     Where C some sparse matrix.
>
>     Once the A is computed I want to solve another problem:
>
>     (A'*A)b=A'*r
>
>     Which is a Gauss-Newton system now,
>     A'*A is truncated Hessian,
>     A'*r is gradient,
>     r = d - d_obs, where d_obs is the real observed data.
>     b -- model change which have to applied to original model m in
>     order to explain your observed data better
>
>     The latter problem I want to solve using petsc matrix-free
>     formulation and some iterative solver (haven't decided yet which,
>     could you advice one?).
>
>     But the POINT here is that the latter problem must be solved not
>     in the original m-space, but in transformed x-space. For that we
>     need another A:
>
>     A_tr = \frac{\partial d_i}{\partial x_j}
>
>     However, using chain rule you can represent A_tr in terms of A:
>
>     A = \frac{\partial d_i}{\partial m_j} * \frac{\partial
>     m_j}{\partial x_j}
>
>     \frac{\partial m_j}{\partial x_j} - is the exactly transformation
>     I want to apply to matrix A. It's a simple scalar expression, but
>     has to be applied to each element of A.
>
>     Did it help? :)
>
>
> Do you ever actually use A? If not, why not just build the transformed 
> operator directly?

Two reasons:
1. It's wrong from the architectural point view. At the point where I 
compute A I prefer know nothing about transformation, just working with 
original physical parameters is much more natural.
2. This original A might be used in the future.

I really suggested that would be easier to apply this transformation to 
the whole matrix. If it's impossible than I have to find another way. 
Just want to be sure of that.

>
>    Matt
>
>     Regards,
>     Alexander
>
>
>
>     On 16.06.2011 14:55, Jed Brown wrote:
>>     On Thu, Jun 16, 2011 at 14:48, Alexander Grayver
>>     <agrayver at gfz-potsdam.de <mailto:agrayver at gfz-potsdam.de>> wrote:
>>
>>         by the grid I meant that this matrix is 2d array, not a real
>>         grid of some physical parameters or whatever and it's also
>>         not a linear operator itself.
>>
>>
>>     Could you explain more about the task you're trying to do.
>>     Pseudocode or Matlab for the whole process would be useful. There
>>     might be a better way to do this without using Mat to store
>>     something that is not a linear operator.
>
>
>
>
> -- 
> 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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110616/68a8bba8/attachment.htm>