# Question:Towards a compiled Newton Raphson

## Question:Towards a compiled Newton Raphson

Acer,

Thanks. This request from me comes from the need that Maple does not have a real DAE solver. (It converts ODEs to DAEs and then solves, which fails for nonlinear and most non-trivial DAEs).

Your question is valid. For situations which can be compiled (i.e, no Heaviside or functions that can't be compiled) at hardware floats, I would like to have a Newton Raphson solver. Some parameters (time values) change when you integrate DAEs, which means that there is a need to solve the same set of nonlinear equations again and again for minor changes in parameters.

While LinearSolve is fast, evaluating F and Jac for a large system of equations can kill the speed. I am attaching an example for Newton wherein F, Jac and a linear solver is compiled. But the Newton procedure is not. If we can compile the Newton procedure also, it will be much more faster.

Also, having a compiled code with parametric form will help in optimization.

 > restart:
 > Digits:=15;
 (1)
 > N:=4;
 (2)
 > f:=array([seq(evalf(-2*x[i]+x[i]^2+0.99/i^2),i=1..N)]);
 (3)
 > fsolve({f[1],f[2],f[3],f[4]},{x[1]=1.2,x[2]=1,x[3]=0,x[4]=0});
 (4)
 >
 > eqsA := [seq(F[i]=f[i],i=1..N)]:
 > irform2 := StatSeq(seq(Assign(op(i)),i=eqsA)):
 > prccons:= codegen[intrep2maple](Proc(Parameters(x::Array,F::Array),irform2)):
 > f3:=Compiler:-Compile(prccons):
 > x0:=Vector([seq(1/2,i=1..N)],datatype=float[8]);
 (5)
 > f0:=Vector([seq(0,i=1..N)],datatype=float[8]);
 (6)
 > eqsJ := [seq(seq(Jac[i,j]=evalf(diff(f[i],x[j])),j=1..N),i=1..N)];
 (7)

 > irformJ := StatSeq(seq(Assign(op(i)),i=eqsJ)):
 > prcconsJ:= codegen[intrep2maple](Proc(Parameters(x::Array,Jac::Array),irformJ)):
 > J3:=Compiler:-Compile(prcconsJ):
 > j0:=Matrix(1..N,1..N,datatype=float[8]):

Following linear sovler algorithm is basiclly from dsolve/numeric/SC

.

 > 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; end proc:
 > s3c:=Compiler:-Compile(s3):
 > x0:=Vector(N,[seq(0.5,i=1..N)],datatype=float[8]);
 (8)
 > f0;
 (9)
 > j0:=Matrix(1..N,1..N,datatype=float[8]);
 (10)
 (11)
 > Newton:=proc(x0,f0,j0,N) global f3,s3c; local xnew,err,xold,ip,i::integer; ip:=Vector(N,[seq(1,i=1..N)],datatype=integer[4]); err:=10; xold:=x0; while err>1e-8 do f3(xold,f0); J3(xold,j0); s3c(N,-j0,ip,f0); xnew:=xold+f0; err:=max(abs(f0))/N; xold:=xnew; end: xnew; end proc;
 >
 (12)
 > Compiler:-Compile(Newton);
 Error, (in CodeGeneration:-IssueError) cannot analyze non-integer range boundary args[4]
 > t11:=time():Newton(x0,f0,j0,N);time()-t11;
 (13)