Segfault 5.9.1 on large matrix operations

hi daniel,

about mnewton ---

> my guess is that solve_the_matrix_equation_for_deltas can be written
> much more efficiently than the method in mnewton (which calculates the
> symbolic inverse, pre-multiplies by the symbolic equations,
> substitutes the guesses and tries to evaluate the whole thing to a
> floating point vector.

mnewton at present calls determinant and adjoint,
which  is likely to be a relatively slow way to solve equations ...

here is a revised version of mnewton which calls
linsolve (gaussian elimination) to solve F(x) + DF(x) h = 0.
it seems to work as intended on the examples in mnewton.dem.
i'll probably commit this after cutting out the debug output.

for what it's worth,
robert dodier


NEWTONEPSILON:10.0^(-fpprec/2)$
NEWTONMAXITER:50$

mnewton(FuncList, VarList, GuessList):=
  block([nfunc, Solutions, Increments, solved:false, h,
         i, j, k, keepfloat:true, ratprint:false],

    /* Depuration of the function arguments */
    GuessList:float(GuessList),
    nfunc:length(FuncList),
    if length(VarList) # nfunc then (
      print("mnewton: incorrect number of variable names (", nfunc,
            "functions but", length(VarList), "variable names)."),
      return(false)
    ),
    if length(GuessList) # nfunc then (
      print("mnewton: incorrect number of approach values (", nfunc,
            "variables but", length(GuessList), "approximation values)."),
      return(false)
    ),
    apply(kill, VarList),

    DF: zeromatrix (nfunc, nfunc),
    h : makelist (h[i], i, 1, nfunc),
    for i:1 thru nfunc do
      for j:1 thru nfunc do
        DF[i][j] : diff (FuncList[i], VarList[j]),  /* Jacobian matrix */

    GuessList: float (GuessList),
print ("DEBUG: initial GuessList:", GuessList),

    /* Solution of the functions system */
    for k:1 thru NEWTONMAXITER do (
      Solutions:map("=", VarList, GuessList),

      /* Solve 0 = F(x) + DF(x) h for h.
       */
      block
        ([DFx: sublis (Solutions, DF),
        Fx: sublis (Solutions, FuncList)],
        Increments: linsolve (transpose (DFx . h + Fx)[1], h)),

      GuessList:GuessList + map (rhs, Increments),
      GuessList: float (GuessList),
print ("DEBUG: Increments:", Increments, "GuessList:", GuessList),

      solved:true,
      for i:1 thru nfunc do
        solved: solved and abs (rhs (Increments[i])) < NEWTONEPSILON,
      if solved then return(0)
    ),

    /* Return of solution or error */
    if solved = false then (
      print("mnewton: the process doesn't converge or it converges too
slowly."),
      return([])
    ),
    Solutions:map("=", VarList, GuessList),
print ("DEBUG: sublis (Solutions, FuncList) =>", sublis (Solutions, FuncList)),
    return([Solutions])
  )$