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