ITL - the correct procedure for solving Ax=b

Hi all!

I got an interesting observation, which I would call a feature rather than a bug. 

In all examples of ITL library for solving equation mA*vx=vb the following listing is suggested (with change SSOR to other solution approaches)


int max_it = 50;

itl::SSOR<Matrix> precond(mA);
itl::noisy_iteration<double> iter(vb, max_it, 0, 1.0e-6);
itl::cgs(mA, vx, vb, precond(), iter);



BUT ! working with vector vb, which has tiny components (say all 1.0e-9, it was my case), all you get after running the procedure is

vx == vb (!!)

(surely, matrix mA is not an identity :] )

So, the problem is in guessing a "good" value of the last parameter for iter() procedure. The way out is simple: one needs to put not *absolute* but *relative* noise, so the correct variant of listing will be


int max_it = 50;

itl::SSOR<Matrix> precond(mA);
// last parameter of iter() is in *relative* to |vb| values
itl::noisy_iteration<double> iter(vb, max_it, 0, 1.0e-6 * itl::two_norm(vb) ); 
itl::cgs(mA, vx, vb, precond(), iter);


You can substitute the   itl::two_norm(vb)   with any value, giving you the order of magnitude of  vb (any norm)

Hope many ppl find it useful.

Dimitri Treebushny
PhD student
UCEWP, Ukraine
www.ucewp.kiev.ua
-- 
___________________________________________________
OperaMail free e-mail - http://www.operamail.com
OperaMail Premium - 28MB, POP3, more! US$29.99/year

Powered by Outblaze
_______________________________________________
This list is archived at http://www.osl.iu.edu/MailArchives/mtl-devel/