restart:

Digits := 15:

N := 50;

N := 50

f := array([seq(evalf(-2*x[i]+x[i]^2+0.99/i^2),i=1..N)]):

CodeTools:-Usage(
  fsolve({seq(f[i],i=1..N)},{seq(x[i]=0.5,i=1..N)})  ):
type(eval(%,1),set(name=numeric));

memory used=43.53MiB, alloc change=64.00MiB, cpu time=2.59s, real time=2.60s, gc time=15.60ms
true

# s3 is not recreated for each new problem, so does not contribute
# to problem specific timing. Even compilation to s3c would be done
# at most once, upon package load say.
#
s3 := proc(n::posint,
           A::Array( datatype = float[8] ),
           ip::Array( datatype = integer[4] ),
           b::Array( datatype = float[8] ) )
   local i::integer, j::integer, k::integer, m::integer, t::float;
      ip[n] := 1;
      for k to n-1 do
        m := k;
        for i from k+1 to n do
          if abs(A[m,k]) < abs(A[i,k]) then
            m := i
           end if
         end do;
        ip[k] := m;
        if m <> k then
          ip[n] := -ip[n]
         end if;
       t := A[m,k];
       A[m,k] := A[k,k];
       A[k,k] := t;
       if t = 0 then
         ip[n] := 0;
         #return ip
         end if;
       for i from k+1 to n do
         A[i,k] := -A[i,k]/t
         end do;
       for j from k+1 to n do
         t := A[m,j];
         A[m,j] := A[k,j];
         A[k,j] := t;
         if t <> 0 then
           for i from k+1 to n do
             A[i,j] := A[i,j]+A[i,k]*t
             end do
           end if
         end do
       end do;
     if A[n,n] = 0 then
       ip[n] := 0
       end if;
     #ip[n]
      if ip[n] = 0 then
        error "The matrix A is numerically singular"
       end if;
      if 1 < n then
        for k to n-1 do
          m := ip[k];
          t := b[m];
          b[m] := b[k];
          b[k] := t;
          for i from k+1 to n do
           b[i] := b[i]+A[i,k]*t
           end do
         end do;
       for k from n by -1 to 2 do
         b[k] := b[k]/A[k,k];
         t := -b[k];
        for i to k-1 do
           b[i] := b[i]+A[i,k]*t
           end do
         end do
       end if;
     b[1] := b[1]/A[1,1];
     #return 0.0; # Better handling of explicit returns needed when inlining
end proc:

# Newtondummy is not recreated for each new problem. But its
# instantiation is done for each problem, and that is accounted
# for in the base timing below.
#
Newtondummy := proc(x::Array(datatype=float[8]),
               f0::Array(datatype=float[8]),
               j0::Array(datatype=float[8]),
               N::integer,
               ip::Array(datatype=integer[4]),
               maxiter::integer[4],
               tol::float[8])
   local ferr::float[8], temp::float[8],
         i::integer, k::integer, iter::integer;
      ferr := 10.0*tol;
      iter := 0;
      while ferr > tol and iter < maxiter do
        _DUMMY1; # f3(x, f0);
        _DUMMY2; # J3(x, j0);
        j0[1,1] := j0[1,1]; # else j0 not recognized as 2-dim?
        _DUMMY3; # s3c(N, j0, ip, f0);
        for i from 1 to N do
           x[i] := x[i] - f0[i];
        end do;
        ferr := abs(f0[1]);
        for i from 2 to N do
           temp := abs(f0[i]);
           if temp > ferr then
              ferr := temp;
           end if;
        end do;
        ferr := ferr;
        iter := iter + 1;
     end do:
     return ferr;
end proc:

# The following are all the problem specific assembly and compilation
# tasks, timed together. (They could be part of an appliable module,
# eventually).
#
(st,str) := time(),time[real]():

eqsA := [seq(F[i]=evalf(f[i]),i=1..N)]:
irform2 := StatSeq(seq(Assign(op(i)),i=eqsA)):
prccons:= codegen[intrep2maple](Proc(Parameters(x::Array,F::Array),irform2)):
eqsJ := remove(type,[seq(seq(Jac[i,j]=evalf(diff(f[i],x[j])),j=1..N),i=1..N)],name=0.0):
irformJ := StatSeq(seq(Assign(op(i)),i=eqsJ)):
prcconsJ:= codegen[intrep2maple](Proc(Parameters(x::Array,Jac::Array),irformJ)):

(basetime,basetimer) := time()-st,time[real]()-str;

basetime, basetimer := 0.31e-1, 0.38e-1

# I need to document what I've done here.
# (I ought to consider writing a more general purpose "function-call Inliner".)
#
(st,str) := time(),time[real]():

UNewtondummy := ToInert(eval(Newtondummy)):
UNdlocals := [op([2,..],UNewtondummy)]:
Uprccons := ToInert(eval(prccons)):
UprcconsJ := ToInert(eval(prcconsJ)):
Us3 := ToInert(eval(s3)):
Us3locals := [op([2,..],Us3)]:
newlocals := [seq(`if`(t::'specfunc(anything,_Inert_DCOLON)',
                  _Inert_DCOLON(_Inert_NAME(cat(op([1,1],t),"_n1")),op([2..],t)),NULL),
                  t=Us3locals)]:
Ls3 := [seq(_Inert_LOCAL(i)=_Inert_LOCAL(i+nops(op(2,UNewtondummy))),i=1..nops(Us3locals))]:
T1 := subs(_Inert_NAME("_DUMMY1")=op([5,..],Uprccons), # param sequence 1st 2nd entries match
           _Inert_NAME("_DUMMY2")=op(subs([_Inert_PARAM(2)=_Inert_PARAM(3)],op([5],UprcconsJ))),
           _Inert_NAME("_DUMMY3")=op(subs([_Inert_PARAM(1)=_Inert_PARAM(4),
                                           _Inert_PARAM(2)=_Inert_PARAM(3),
                                           _Inert_PARAM(3)=_Inert_PARAM(5),
                                           _Inert_PARAM(4)=_Inert_PARAM(2),
                                           op(Ls3)],op([5],Us3))),
           subsop(2=_Inert_LOCALSEQ(op([2,..],UNewtondummy),op(newlocals)),UNewtondummy)):

# This is like `Newton`, but with the code from prccons,prcconsJ,s3 inlined.
InlinedNewton := FromInert(T1):

InlinedNewtonc := Compiler:-Compile(InlinedNewton):

(fusedcompiletime,fusedcompiletimer) := time()-st,time[real]()-str;

fusedcompiletime, fusedcompiletimer := .983, 1.526

# And now, essentially what was done previously. Ie, all of J3, f3, and s3c are
# Compiled. And so is `Newton`, except now it is accessing those three lexically
# while before they were explicitly declared global in `Newton`.
#
(st,str) := time(),time[real]():

f3 := Compiler:-Compile(prccons):
J3 := Compiler:-Compile(prcconsJ):
s3c := Compiler:-Compile(s3):
Newton := subs([_DUMMY1='f3(x, f0)',
                _DUMMY2='J3(x, j0)',
                _DUMMY3='s3c(N, j0, ip, f0)'],eval(Newtondummy)):
Newtonc := Compiler:-Compile(Newton):

(separatedcompiletime,separatedcompiletimer) := time()-st,time[real]()-str;

Warning, the function names {J3, f3, s3c} are not recognized in the target language

separatedcompiletime, separatedcompiletimer := .780, 1.306

# These need only be created for the first problem, and could be
# efficiently grown (but not replaced) for new problems with larger size N.
#
x0 := Array(1..N,datatype=float[8]):
f0 := Array(1..N,datatype=float[8]):
j0 := Array(1..N,1..N,datatype=float[8]):
ip := Array(1..N,datatype=integer[4]):

repeats := 1000:

ArrayTools:-Fill(0.0,f0):
ArrayTools:-Fill(1,ip):
(t11,t11r) := time(),time[real]():
for ii from 1 to repeats do
  ArrayTools:-Fill(0.5,x0):
  ArrayTools:-Fill(0.0,j0):
  InlinedNewtonc(x0,f0,j0,N,ip,100,HFloat(2e-6)); # singly compiled
end do:
(tottime,totrtime) := time()-t11,time[real]()-t11r:
evalf[4](tottime/repeats)*'seconds',evalf[4](totrtime/repeats)*'seconds[real]';

ArrayTools:-Fill(0.5,x0):
ArrayTools:-Fill(0.0,j0):
InlinedNewtonc(x0,f0,j0,N,ip,100,HFloat(2e-6)); # singly compiled

0.8420e-3*seconds, 0.8500e-3*seconds[real]
0.463564004320658215e-6

ArrayTools:-Fill(0.0,f0):
ArrayTools:-Fill(1,ip):
(t11,t11r) := time(),time[real]():
for ii from 1 to repeats do
  ArrayTools:-Fill(0.5,x0):
  ArrayTools:-Fill(0.0,j0):
  Newtonc(x0,f0,j0,N,ip,100,HFloat(2e-6)); # separetly compiled
end do:
(tottime,totrtime) := time()-t11,time[real]()-t11r:
evalf[4](tottime/repeats)*'seconds',evalf[4](totrtime/repeats)*'seconds[real]';

ArrayTools:-Fill(0.5,x0):
ArrayTools:-Fill(0.0,j0):
Newtonc(x0,f0,j0,N,ip,100,HFloat(2e-6)); # separately compiled

0.1077e-2*seconds, 0.1087e-2*seconds[real]
0.463564004320658215e-6