/
polydiv.cc
63 lines (48 loc) · 1.36 KB
/
polydiv.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#include "seqspace.hh"
#include "multop.hh"
#include "cgnealg.hh"
#include "terminator.hh"
using namespace RVL;
using namespace RVLAlg;
using namespace RVLUmin;
int main() {
cout <<"Polynomial Division in l2 - computing least squares\n";
cout <<"polynomial approximation to 1/(1+ax) via CG.\n";
double a;
cout <<"Enter coeff a of x: ";
cin >> a;
list<double> alist;
alist.push_back(1.0);
alist.push_back(a);
// cerr<<"op of multiplication by 1+x"<<endl;
PolyMultOp op(alist);
// cerr<<"vector in domain "<<endl;
Vector<double> x(op.getDomain());
// cerr<<"vector in range"<<endl;
Vector<double> b(op.getRange());
// cerr<<"assign to b the coeff list of the const poly 1"<<endl;
list<double> blist;
blist.push_back(1.0);
AssignList initb(blist);
b.eval(initb);
// cerr<<"zero out x"<<endl;
x.zero();
int nsteps;
double tol;
cout<<"Enter limit on number of CG steps: ";
cin >> nsteps;
cout<<"Enter tolerance for norm of residual: ";
cin >> tol;
if (nsteps < 1) nsteps = 50;
if (tol < 100*numeric_limits<double>::epsilon())
tol = 100*numeric_limits<double>::epsilon();
// cerr<<"create CG step"<<endl;
double rnorm;
double nrnorm;
CGNEAlg<double> cg(x,op,b,rnorm,nrnorm,tol,tol,nsteps);
// cerr<<"run"<<endl;
cg.run();
cout<<"\nComputed Solution:\n";
x.write(cout);
// cerr<<"bye-bye\n";
}