On Fri, May 2, 2014 at 12:53 PM, Xiangdong <epsco...@gmail.com> wrote:
>
> On Thu, May 1, 2014 at 10:21 PM, Barry Smith <bsm...@mcs.anl.gov> wrote:
>
>>
>> On May 1, 2014, at 9:12 PM, Xiangdong <epsco...@gmail.com> wrote:
>>
>> > I came up with a simple example to demonstrate this "eliminating row"
>> behavior. It happens when the solution x to the linearized equation Ax=b is
>> out of the bound set by SNESVISetVariableBounds();
>> >
>> > In the attached example, I use snes to solve a simple function x-b=0.
>> When you run it, it outputs the matrix as 25 rows, while the real Jacobian
>> should be 5*5*2=50 rows. If one changes the lower bound in line 125 to be
>> -inf, it will output 50 rows for the Jacobian. In the first case, the norm
>> given by SNESGetFunctionNorm and SNESGetFunction+VecNorm are also different.
>> >
>> > In solving the nonlinear equations, it is likely that the solution of
>> the linearized equation is out of bound, but then we can reset the
>> out-of-bound solution to be lower or upper bound instead of eliminating the
>> variables (the rows). Any suggestions on doing this in petsc?
>>
>> This is what PETSc is doing. It is using the "active set method".
>> Variables that are at their bounds are “frozen” and then a smaller system
>> is solved (involving just the variables not a that bounds) to get the next
>> search direction. Based on the next search direction some of the variables
>> on the bounds may be unfrozen and other variables may be frozen. There is a
>> huge literature on this topic. See for example our buddies • Nocedal,
>> Jorge; Wright, Stephen J. (2006). Numerical Optimization (2nd ed.). Berlin,
>> New York: Springer-Verlag. ISBN 978-0-387-30303-1..
>>
>> The SNESGetFunctionNorm and SNESGetFunction+VecNorm may return
>> different values with the SNES VI solver. If you care about the function
>> value just use SNESGetFunction() and compute the norm that way. We are
>> eliminating SNESGetFunctionNorm() from PETSc because it is problematic.
>>
>> If you think the SNES VI solver is actually not solving the problem,
>> or giving the wrong answer than please send us the entire simple code and
>> we’ll see if we have introduced any bugs into our solver. But note that the
>> linear system being of different sizes is completely normal for the solver.
>>
>
> Here is an example I do not quite understand. I have a simple function
> F(X) = [x1+x2+100 ; x1-x2; x3+x4; x3-x4+50]. If I solve this F(X)=0 with no
> constraint, the exact solution is [x1=-50; x2=-50; x3=-25; x4=25].
>
> If I specify the constraint as x2>=0 and x4>=0, I expect the solution from
> one iteration of SNES is [-50, 0,-25,25], since the constraint on x2 should
> be active now. However, the petsc outputs the solution [-50, 0, 0, 0].
> Since x3 and x4 does not violate the constraint, why does the solution of
> x3 and x4 change (note that x3 and x3 are decoupled from x1 and x2)? In
> this case, the matrix obtained from KSPGetOperators is only 2-by-2, so two
> variables or constraints are eliminated.
>
This just finds a local solution to the constrained problem, and these need
not be unique.
Matt
> Another thing I noticed is that constraints x2>-1e-7 and x4>-1e-7 gives
> solution [-50,-1e-7,-25,25]; however, constraints x2>-1e-8 and x4>-1e-8
> gives the solution [-50,0,0,0].
>
> Attached please find the simple 130-line code showing this behavior.
> Simply commenting the line 37 to remove the constraints and modifying line
> 92 to change the lower bounds of x2 and x4.
>
> Thanks a lot for your time and help.
>
> Best,
> Xiangdong
>
>
>
>
>>
>>
>> Barry
>>
>>
>> >
>> > Thank you.
>> >
>> > Best,
>> > Xiangdong
>> >
>> > P.S. If we change the lower bound of field u (line 124) to be zero,
>> then the Jacobian matrix is set to be NULL by petsc.
>> >
>> >
>> > On Thu, May 1, 2014 at 3:43 PM, Xiangdong <epsco...@gmail.com> wrote:
>> > Here is the order of functions I called:
>> >
>> > DMDACreate3d();
>> >
>> > SNESCreate();
>> >
>> > SNESSetDM(); (DM with dof=2);
>> >
>> > DMSetApplicationContext();
>> >
>> > DMDASNESSetFunctionLocal();
>> >
>> > SNESVISetVariableBounds();
>> >
>> > DMDASNESetJacobianLocal();
>> >
>> > SNESSetFromOptions();
>> >
>> > SNESSolve();
>> >
>> > SNESGetKSP();
>> > KSPGetSolution();
>> > KSPGetRhs();
>> > KSPGetOperators(); //get operator kspA, kspx, kspb;
>> >
>> > SNESGetFunctionNorm(); ==> get norm fnorma;
>> > SNESGetFunction(); VecNorm(); ==> get norm fnormb;
>> > SNESComputeFunction(); VecNorm(); ==> function evaluation fx at the
>> solution x and get norm fnormc;
>> >
>> > Inside the FormJacobianLocal(), I output the matrix jac and preB;
>> >
>> > I found that fnorma matches the default SNES monitor output "SNES
>> Function norm", but fnormb=fnormc != fnorma. The solution x, the residue fx
>> obtained by snescomputefunction, mat jac and preB are length 50 or
>> 50-by-50, while the kspA, kspx, kspb are 25-by-25 or length 25.
>> >
>> > I checked that kspA=jac(1:2:end,1:2:end) and x(1:2:end)= kspA\kspb;
>> x(2:2:end)=0; It seems that it completely ignores the second degree of
>> freedom (setting it to zero). I saw this for (close to) constant initial
>> guess, while for heterogeneous initial guess, it works fine and the matrix
>> and vector size are correct, and the solution is correct. So this
>> eliminating row behavior seems to be initial guess dependent.
>> >
>> > I saw this even if I use snes_fd, so we can rule out the possibility of
>> wrong Jacobian. For the FormFunctionLocal(), I checked via
>> SNESComputeFunction and it output the correct vector of residue.
>> >
>> > Are the orders of function calls correct?
>> >
>> > Thank you.
>> >
>> > Xiangdong
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> > On Thu, May 1, 2014 at 1:58 PM, Barry Smith <bsm...@mcs.anl.gov> wrote:
>> >
>> > On May 1, 2014, at 10:32 AM, Xiangdong <epsco...@gmail.com> wrote:
>> >
>> > > Under what condition, SNESGetFunctionNorm() will output different
>> results from SENEGetFunction + VecNorm (with NORM_2)?
>> > >
>> > > For most of my test cases, it is the same. However, when I have some
>> special (trivial) initial guess to the SNES problem, I see different norms.
>> >
>> > Please send more details on your “trivial” case where the values are
>> different. It could be that we are not setting the function norm properly
>> on early exit from the solvers.
>> > >
>> > > Another phenomenon I noticed with this is that KSP in SNES squeeze my
>> matrix by eliminating rows. I have a Jacobian supposed to be 50-by-50. When
>> I use KSPGetOperators/rhs/solutions, I found that the operator is 25-by-25,
>> and the rhs and solution is with length 25. Do you have any clue on what
>> triggered this? To my surprise, when I output the Jacobian inside the
>> FormJacobianLocal, it outputs the correct matrix 50-by-50 with correct
>> numerical entries. Why does the operator obtained from KSP is different and
>> got rows eliminated? These rows got eliminated have only one entries per
>> row, but the rhs in that row is not zero. Eliminating these rows would give
>> wrong solutions.
>> >
>> > Hmm, we never squeeze out rows/columns from the Jacobian. The size
>> of the Jacobian set with SNESSetJacobian() should always match that
>> obtained with KSPGetOperators() on the linear system. Please send more
>> details on how you get this. Are you calling the KSPGetOperators() inside a
>> preconditioner where the the preconditioner has chopped up the operator?
>> >
>> > Barry
>> >
>> > >
>> > > Thank you.
>> > >
>> > > Xiangdong
>> > >
>> > >
>> > >
>> > >
>> > >
>> > >
>> > > On Tue, Apr 29, 2014 at 3:12 PM, Matthew Knepley <knep...@gmail.com>
>> wrote:
>> > > On Tue, Apr 29, 2014 at 2:09 PM, Xiangdong <epsco...@gmail.com>
>> wrote:
>> > > It turns out to a be a bug in my FormFunctionLocal(DMDALocalInfo
>> *info,PetscScalar **x,PetscScalar **f,AppCtx *user). I forgot to initialize
>> the array f. Zero the array f solved the problem and gave consistent result.
>> > >
>> > > Just curious, why does not petsc initialize the array f to zero by
>> default inside petsc when passing the f array to FormFunctionLocal?
>> > >
>> > > If you directly set entires, you might not want us to spend the time
>> writing those zeros.
>> > >
>> > > I have another quick question about the array x passed to
>> FormFunctionLocal. If I want to know the which x is evaluated, how can I
>> output x in a vector format? Currently, I created a global vector vecx and
>> a local vector vecx_local, get the array of vecx_local_array, copy the x to
>> vecx_local_array, scatter to global vecx and output vecx. Is there a quick
>> way to restore the array x to a vector and output?
>> > >
>> > > I cannot think of a better way than that.
>> > >
>> > > Matt
>> > >
>> > > Thank you.
>> > >
>> > > Best,
>> > > Xiangdong
>> > >
>> > >
>> > >
>> > > On Mon, Apr 28, 2014 at 10:28 PM, Barry Smith <bsm...@mcs.anl.gov>
>> wrote:
>> > >
>> > > On Apr 28, 2014, at 3:23 PM, Xiangdong <epsco...@gmail.com> wrote:
>> > >
>> > > > Hello everyone,
>> > > >
>> > > > When I run snes program,
>> > >
>> > > ^^^^ what SNES program”?
>> > >
>> > > > it outputs "SNES Function norm 1.23456789e+10". It seems that this
>> norm is different from residue norm (even if solving F(x)=0)
>> > >
>> > > Please send the full output where you see this.
>> > >
>> > > > and also differ from norm of the Jacobian. What is the definition
>> of this "SNES Function Norm”?
>> > >
>> > > The SNES Function Norm as printed by PETSc is suppose to the
>> 2-norm of F(x) - b (where b is usually zero) and this is also the same
>> thing as the “residue norm”
>> > >
>> > > Barry
>> > >
>> > > >
>> > > > Thank you.
>> > > >
>> > > > Best,
>> > > > Xiangdong
>> > >
>> > >
>> > >
>> > >
>> > >
>> > > --
>> > > 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
>> > >
>> >
>> >
>> >
>> > <exdemo.c>
>>
>>
>
--
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
