<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Wed, May 13, 2015 at 4:59 PM, Justin Chang <span dir="ltr"><<a href="mailto:jychang48@gmail.com" target="_blank">jychang48@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>Okay that makes sense. Reason for asking this is because I am solving the following convex optimization problem for a linear diffusion problem:<br><br></div><div>min 1/2 x^T * A * x + x^T * g<br></div><div>s.t. 0 <= x <= 1<br><br></div><div>Using Tao's BLMVM solver, the objective function and gradient vector reads as follows:<br><br></div><div>f = 0.5*x^T*A*x + x^T*(r - A*x_0)<br></div><div>g = A*(x - x_0) + r<br><br></div><div>Where A, x_0, and r are the jacobian, initial guess, and residual respectively. To compute these, I have TaoSetObjectiveAndGradientRoutine() invoke the following function:<br><br>#undef __FUNCT__<br>#define __FUNCT__ "FormFunctionGradient"<br>PetscErrorCode FormFunctionGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx)<br>{<br> AppCtx *user = (AppCtx*)ctx;<br> PetscScalar xtHx;<br> PetscErrorCode ierr;<br><br> PetscFunctionBegin;<br> ierr = MatMult(user->A,x,g);CHKERRQ(ierr);<br> ierr = VecDot(x,g,&xtHx);CHKERRQ(ierr);<br> ierr = VecDot(x,user->f,f);CHKERRQ(ierr);<br> *f += 0.5*xtHx;<br> ierr = VecAXPY(g,1.0,user->f);CHKERRQ(ierr);<br> PetscFunctionReturn(0);<br>}<br><br></div><div>where user->f = (r - A*x_0) The timing results from -log_summary seem to indicate that MatMult(user->A,x,g) is a bottleneck within TaoSolve() so I was wondering if I could somehow make the computation of Vec g matrix-free considering that user->A is constant. The problem proliferates if I decide to use TRON. Or is this unavoidable (Jason any thoughts?)<br></div></div></blockquote><div> </div><div>If you have a first order discretization, its still better to form the matrix since it uses fewer flops on repeated application.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div></div><div>Thanks,<br></div><div>Justin<br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Wed, May 13, 2015 at 3:15 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><span>On Wed, May 13, 2015 at 2:54 PM, Justin Chang <span dir="ltr"><<a href="mailto:jychang48@gmail.com" target="_blank">jychang48@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div><div>Hello everyone,<br><br></div>From what I read online, one could use a matrix-free method if using the conjugate gradient method for a symmetric and positive definite problem. Some of the PETSc SNES examples like ex12 compute the Jacobian A explicitly, but I was wondering if it's possible to convert this problem into matrix free form. That is, to directly assemble the product of A*x into vector form without first forming A and then invoking MatMult().<br></div></div></div></blockquote><div><br></div></span><div>Yes, Jed advocates this all the time, but there is some infrastructure that you want for this</div><div>to be really effective. For quadratic and higher orders, it makes sense to compute the action</div><div>of your operator MF. This is done by defining something that looks a lot like your residual</div><div>evaluation (element vec in and out) but computes the linearization. However, you still usually</div><div>need a low order explicit matrix for preconditioning. I have not spent time hooking that all together.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div></div>Thanks,<br></div>Justin<span><font color="#888888"><br></font></span></div><span><font color="#888888">
</font></span></blockquote></div><span><font color="#888888"><br><br clear="all"><span class="HOEnZb"><font color="#888888"><div><br></div>-- <br><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</font></span></font></span></div></div>
</blockquote></div><br></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</div></div>