$title Asynchronous Jacobi Methods (JACOBI,SEQ=329)

$onText
This example outlines procedures for implementing various serial and
parallel iterative schemes. For simplicity, a system of linear equations
is selected. This schema extends naturally to other problem types like
nonlinear systems and mixed complementarity problems.

We will implement various ways to solve the problem:

   Gauss Seidel      serial
   Jacobi - parallel sub problems
   Gauss-Seidel      Asynchronous

Bertsekas, D P, and Tsitsiklis, J N, Parallel and distributed
computation: numerical methods. Prentice-Hall, Inc., Upper Saddle
River, NJ, USA, 1989.

Keywords: mixed complementarity problem, Jacobi method, Gauss Seidel method,
          parallel algorithms, dynamic programming, distributed algorithms,
          mathematics
$offText

$ifE %system.gamsversion%<231 $abort Version too old! Need 231 or higher
$eolCom //

$setDDList vars parts iters       // acceptable double dash parameters
$if not set vars  $set vars   50  // number of variables
$if not set parts $set parts   5  // number of partitions
$if not set iters $set iters 100  // max number of iterations
$if not errorfree $abort wrong double dash parameters: --vars=n --parts=n iters=n

Set i 'problem size' / i1*i%vars% /;

Alias (i,j);

Variable x(i);

Equation e(i);

Parameter A(i,j), b(i);

e(i).. sum(j, A(i,j)*x(j)) =e= b(i);

Model lin / e.x /;

b(i)   = 1;
A(i,i) = 1;
A(i,j)$(not sameas(i,j)) = 0.001;

lin.solPrint  = %solPrint.quiet%;       // suppress solution output
lin.solveLink = %solveLink.callModule%; // keep gams memory resident

Set
   iters       'iteration count'           / iter0*iter%iters%     /
   k           'problem partition blocks'  / block_1*block_%parts% /
   active(k,i) 'active vars in partition k'
   fixed(k,i)  'fixed vars in partition k';

Alias (kp,k);

Parameter
   resrep(iters,*) 'summary residual report'
   solrep(i,*)     'summary solution report'
   stats           'summary statistics'
   res(iters)      'max residual'
   h(k)            'handles'
   tol             'convergence tolerance' / 1e-4 /
   iter            'iteration counter'
   curres          'intermediate residual values'
   t1              'temporary timer vars';

active(k,i) = ceil(ord(i)*card(k)/card(i)) = ord(k);
fixed(k,i)  = not active(k,i);

**** solve big problem
t1 = TimeElapsed;

solve lin using mcp;
stats('elapsed','Big Problem') = TimeElapsed - t1;
stats('solves' ,'Big Problem') = 1;

solrep(i,'Big Problem') = x.l(i);

**** Gauss Seidel - all serial
x.l(i)       = 0;
res(iters)   = 0;
res('iter0') = smax(i, abs(b(i)));

t1 = TimeElapsed;
loop(iters$(res(iters) > tol),
   loop(k,
      x.fx(i)$fixed(k,i) = x.l(i);
      solve lin using mcp;
      x.lo(i)$fixed(k,i) = -inf;
      x.up(i)$fixed(k,i) =  inf;
   );
   res(iters+1) = smax(i, abs(b(i) - sum(j, A(i,j)*x.l(j))));
);

stats('elapsed','Gauss Seidel') = TimeElapsed - t1;
stats('solves' ,'Gauss Seidel') = (card(res) - 1)*card(k);

resrep(iters,'Gauss Seidel') = res(iters);
solrep(i,'Gauss Seidel')     = x.l(i);

**** Jacobi - parallel sub problems
lin.solveLink = %solveLink.asyncGrid%; // set grid mode

x.l(i)       = 0;
res(iters)   = 0;
res('iter0') = smax(i, abs(b(i)));

parameter xsave(i);
t1 = TimeElapsed;
loop(iters$(res(iters) > tol),
   loop(k,                            // submitting loop
      x.fx(i)$fixed(k,i) = x.l(i);
      solve lin using mcp;
      h(k) = lin.handle;
      x.lo(i)$fixed(k,i) = -inf;
      x.up(i)$fixed(k,i) =  inf;
   );
   repeat                             // collection loop
      loop(k$handlecollect(h(k)),
         xsave(i)$active(k,i) = x.l(i); // merge in relevant slice of x
         display$handledelete(h(k)) 'trouble deleting handle';
         h(k) = 0;                    // mark problem as solved
      );
      display$sleep(card(h)*0.1) ' sleep a bit';
   until card(h) = 0 or timeelapsed > 10;
   x.l(i) = xsave(i);
   res(iters+1) = smax(i, abs(b(i) - sum(j, A(i,j)*x.l(j))));
);

stats('elapsed','Jacobi') = TimeElapsed - t1;
stats('solves' ,'Jacobi') = (card(res) - 1)*card(k);
resrep(iters,'Jacobi')    = res(iters);
solrep(i,'Jacobi')        = x.l(i);

**** Asynchronous Gauss-Seidel
lin.solveLink = %solveLink.asyncGrid%;   // set grid mode

x.l(i)       = 0;
res(iters)   = 0;
res('iter0') = smax(i, abs(b(i)));
iter = 0;
t1 = TimeElapsed;

loop(k,                          // initial submission loop
   x.fx(i)$fixed(k,i) = x.l(i);
   solve lin using mcp;
   h(k) = lin.handle;
   x.lo(i)$fixed(k,i) = -inf;
   x.up(i)$fixed(k,i) =  inf;
);

xsave(i) = x.l(i);
repeat   // retrieve and submit
   loop(k$handlecollect(h(k)),
      xsave(i)$active(k,i) = x.l(i); // merge in relevant slice of x
      x.l(i) = xsave(i);
      display$handledelete(h(k)) 'trouble deleting handle';
      h(k)   = 0;
      iter   = iter + 1;
      curres = smax(i, abs(b(i) - sum(j, A(i,j)*x.l(j))));
      res(iters)$(ord(iters) = iter + 1) = curres;
      if(curres > tol,
         loop(kp$(h(kp) = 0 and smax(active(kp,i), abs(b(i) - sum(j, A(i,j)*x.l(j)))) > tol),
            x.fx(i)$fixed(kp,i) = x.l(i);
            solve lin using mcp;  // submit new problem
            h(kp) = lin.handle;
            x.lo(i)$fixed(kp,i) = -inf;
            x.up(i)$fixed(kp,i) =  inf;
         );
      );
   );
   display$sleep(card(h)*0.1) 'sleep a bit', curres;
until card(h) = 0 or timeelapsed > 100;

stats('elapsed','Async Gauss') = TimeElapsed - t1;
stats('solves' ,'Async Gauss') = card(res) - 1;

resrep(iters,'Async Gauss') = res(iters);
solrep(i,'Async Gauss')     = x.l(i);

option dispWidth = 15, decimals = 6;

display resrep, solrep, stats;