Advanced Search

Question: Towards a compiled Newton Raphson

Question:Towards a compiled Newton Raphson

Dr. Venkat Subramanian 406
May 11 2015
2

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;

Digits := 15

 (1)
> 

N:=4;

N := 4

 (2)
> 

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

f := Vector[row](4, {(1) = -2.*x[1]+x[1]^2+.99, (2) = -2.*x[2]+x[2]^2+.247500000000000, (3) = -2.*x[3]+x[3]^2+.110000000000000, (4) = -2.*x[4]+x[4]^2+0.618750000000000e-1})

 (3)
> 

fsolve({f[1],f[2],f[3],f[4]},{x[1]=1.2,x[2]=1,x[3]=0,x[4]=0});

{x[1] = .900000000000000, x[2] = .132532421355126, x[3] = 1.94339811320566, x[4] = 0.314314686094742e-1}

 (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]);

x0 := Vector(4, {(1) = .500000000000000, (2) = .500000000000000, (3) = .500000000000000, (4) = .500000000000000})

 (5)
> 

f0:=Vector([seq(0,i=1..N)],datatype=float[8]);

f0 := Vector(4, {(1) = 0., (2) = 0., (3) = 0., (4) = 0.})

 (6)
> 

eqsJ := [seq(seq(Jac[i,j]=evalf(diff(f[i],x[j])),j=1..N),i=1..N)];

eqsJ := [Jac[1, 1] = -2.+2.*x[1], Jac[1, 2] = 0., Jac[1, 3] = 0., Jac[1, 4] = 0., Jac[2, 1] = 0., Jac[2, 2] = -2.+2.*x[2], Jac[2, 3] = 0., Jac[2, 4] = 0., Jac[3, 1] = 0., Jac[3, 2] = 0., Jac[3, 3] = -2.+2.*x[3], Jac[3, 4] = 0., Jac[4, 1] = 0., Jac[4, 2] = 0., Jac[4, 3] = 0., Jac[4, 4] = -2.+2.*x[4]]

 (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]);

x0 := Vector(4, {(1) = .500000000000000, (2) = .500000000000000, (3) = .500000000000000, (4) = .500000000000000})

 (8)
> 

f0;

Vector(4, {(1) = 0., (2) = 0., (3) = 0., (4) = 0.})

 (9)
> 

j0:=Matrix(1..N,1..N,datatype=float[8]);

j0 := Matrix(4, 4, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (2, 1) = 0., (2, 2) = 0., (2, 3) = 0., (2, 4) = 0., (3, 1) = 0., (3, 2) = 0., (3, 3) = 0., (3, 4) = 0., (4, 1) = 0., (4, 2) = 0., (4, 3) = 0., (4, 4) = 0.})

 (10)

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;
> 

 

Newton := proc (x0, f0, j0, N) local xnew, err, xold, ip, i::integer; global f3, s3c; ip := Vector(N, [seq(1, i = 1 .. N)], datatype = integer[4]); err := 10; xold := x0; while 0.1e-7 < err do f3(xold, f0); J3(xold, j0); s3c(N, -j0, ip, f0); xnew := xold+f0; err := max(abs(f0))/N; xold := xnew end do; 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;

Vector(4, {(1) = .900000000000000, (2) = .132532421355126, (3) = 0.566018867943396e-1, (4) = 0.314314686094742e-1})

0.

 (13)


Download NewtonRaphsoncompileattempt.mws
Please Wait...