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

Matthew Knepley knepley at gmail.com
Thu Jun 16 09:07:48 CDT 2011


On Thu, Jun 16, 2011 at 1:43 PM, Alexander Grayver
<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?

   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>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/4851a43b/attachment-0001.htm>


More information about the petsc-users mailing list