jacobi.gms : Asynchronous Jacobi Methods

Description

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

Reference

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

Large Model of Type : MCP

Category : GAMS Model library

Main file : jacobi.gms

``````\$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.holdFixed = 1;                      // treat fixed vars as constants

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

x.l(i)       = 0;
res(iters)   = 0;
res('iter0') = smax(i, abs(b(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)),
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;
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,'Jacopi')        = x.l(i);

**** Asynchronous Gauss-Seidel

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

repeat   // retriev and submit
loop(k\$handlecollect(h(k)),
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;
``````