From: pratip on
Dear Group,

Here is a simple code to solve linear systems. It is motivated from
Fullerton University numerics example. Now if someone can suggest how
to parallelize this code.

SORmethod[A0_, B0_, P0_, omega_, max_] :=
Module[{A = N[A0], B = N[B0], i, j, k = 0, n = Length[P0], P = P0,
oldP = P0},
While[ k < max,
Do[P[[i]] = (1 - omega) oldP[[i]] +
omega/A[[i,
i]] (B[[i]] - Sum[A[[i, j]] P[[j]], {j, 1, i - 1}] -
Sum[A[[i, j]] oldP[[j]], {j, 1 + i, n}]) , {i, 1, n}];
oldP = P;
k = k + 1; ];
Return[P] ]

To test the solver.

n = 10;
A = DiagonalMatrix[Table[2., {n}]];
For[i = 1, i <= n - 1, i++,
A[[i, i + 1]] = A[[i + 1, i]] = 1.;];
B = Table[5. - Abs[i - 5], {i, 1, n}];
P = Table[1., {i, n}];
X3 = SORmethod[A, B, P, 1.5, 200] // Chop

Check the above result with Mathematica LinearSolve

LinearSolve[A, B] // N

Straightforward ParallelDo or ParallelSum is not working. One can see
for large n (abt 400-4000) we will significantly get help if we can
parallize the code.

Best regards,

Pratip

From: Jon Harrop on
"pratip" <pratip.chakraborty(a)gmail.com> wrote in message
news:i13tso$7j1$1(a)smc.vnet.net...
> Dear Group,
>
> Here is a simple code to solve linear systems. It is motivated from
> Fullerton University numerics example. Now if someone can suggest how
> to parallelize this code.
>
> SORmethod[A0_, B0_, P0_, omega_, max_] :=
> Module[{A = N[A0], B = N[B0], i, j, k = 0, n = Length[P0], P = P0,
> oldP = P0},
> While[ k < max,
> Do[P[[i]] = (1 - omega) oldP[[i]] +
> omega/A[[i,
> i]] (B[[i]] - Sum[A[[i, j]] P[[j]], {j, 1, i - 1}] -
> Sum[A[[i, j]] oldP[[j]], {j, 1 + i, n}]) , {i, 1, n}];
> oldP = P;
> k = k + 1; ];
> Return[P] ]

Frigo et al. described a cache oblivious divide-and-conquer algorithm that
parallelizes such functions effectively but you'll need modern
infrastructure for load-balanced parallelism on multicores.

I don't really see the point of parallelizing that in Mathematica. The
serial implementation might be useful as an educational exercise but it will
be orders of magnitude slower than necessary and Mathematica cannot express
an efficient parallel solution.

Cheers,
Jon.


 | 
Pages: 1
Prev: Displaying cylinders
Next: FittedModel object type