Cyclops Tensor Framework
parallel arithmetic on multidimensional arrays
jacobi.cxx
Go to the documentation of this file.
1 
8 #include <ctf.hpp>
9 using namespace CTF;
10 
11 // compute a single Jacobi iteration to get new x, elementwise: x_i <== d_i*(b_i-sum_j R_ij*x_j)
12 // solves Ax=b where R_ij=A_ij for i!=j, while R_ii=0, and d_i=1/A_ii
14  x["i"] = -R["ij"]*x["j"];
15  x["i"] += b["i"];
16  x["i"] *= d["i"];
17 }
18 
19 int jacobi(int n,
20  World & dw){
21 
22  Matrix<> spA(n, n, SP, dw, "spA");
23  Matrix<> dnA(n, n, dw, "dnA");
24  Vector<> b(n, dw);
25  Vector<> c1(n, dw);
26  Vector<> c2(n, dw);
27  Vector<> res(n, dw);
28 
29  srand48(dw.rank);
30  b.fill_random(0.0,1.0);
31  c1.fill_random(0.0,1.0);
32  c2["i"] = c1["i"];
33 
34  //make diagonally dominant matrix
35  dnA.fill_random(0.0,1.0);
36  spA["ij"] += dnA["ij"];
37  //sparsify
38  spA.sparsify(.5);
39  spA["ii"] += 2.*n;
40  dnA["ij"] = spA["ij"];
41 
42  Vector<> d(n, dw);
43  d["i"] = spA["ii"];
44  Transform<> inv([](double & d){ d=1./d; });
45  inv(d["i"]);
46 
47  Matrix<> spR(n, n, SP, dw, "spR");
48  Matrix<> dnR(n, n, dw, "dnR");
49  spR["ij"] = spA["ij"];
50  dnR["ij"] = dnA["ij"];
51  spR["ii"] = 0;
52  dnR["ii"] = 0;
53 
54 /* spR.print();
55  dnR.print(); */
56 
57  //do up to 100 iterations
58  double res_norm;
59  int iter;
60  for (iter=0; iter<100; iter++){
61  jacobi_iter(dnR, b, d, c1);
62 
63  res["i"] = b["i"];
64  res["i"] -= dnA["ij"]*c1["j"];
65 
66  res_norm = res.norm2();
67  if (res_norm < 1.E-4) break;
68  }
69 #ifndef TEST_SUITE
70  if (dw.rank == 0)
71  printf("Completed %d iterations of Jacobi with dense matrix, residual F-norm is %E\n", iter, res_norm);
72 #endif
73 
74  for (iter=0; iter<100; iter++){
75  jacobi_iter(spR, b, d, c2);
76 
77  res["i"] = b["i"];
78  res["i"] -= spA["ij"]*c2["j"];
79 
80  res_norm = res.norm2();
81  if (res_norm < 1.E-4) break;
82  }
83 #ifndef TEST_SUITE
84  if (dw.rank == 0)
85  printf("Completed %d iterations of Jacobi with sparse matrix, residual F-norm is %E\n", iter, res_norm);
86 #endif
87 
88  c2["i"] -= c1["i"];
89 
90  bool pass = c2.norm2() <= 1.E-6;
91 
92  if (dw.rank == 0){
93  if (pass)
94  printf("{ Jacobi x[\"i\"] = (1./A[\"ii\"])*(b[\"j\"] - (A[\"ij\"]-A[\"ii\"])*x[\"j\"]) with sparse A } passed \n");
95  else
96  printf("{ Jacobi x[\"i\"] = (1./A[\"ii\"])*(b[\"j\"] - (A[\"ij\"]-A[\"ii\"])*x[\"j\"]) with sparse A } failed \n");
97  }
98  return pass;
99 }
100 
101 
102 #ifndef TEST_SUITE
103 char* getCmdOption(char ** begin,
104  char ** end,
105  const std::string & option){
106  char ** itr = std::find(begin, end, option);
107  if (itr != end && ++itr != end){
108  return *itr;
109  }
110  return 0;
111 }
112 
113 
114 int main(int argc, char ** argv){
115  int rank, np, n, pass;
116  int const in_num = argc;
117  char ** input_str = argv;
118 
119  MPI_Init(&argc, &argv);
120  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
121  MPI_Comm_size(MPI_COMM_WORLD, &np);
122 
123  if (getCmdOption(input_str, input_str+in_num, "-n")){
124  n = atoi(getCmdOption(input_str, input_str+in_num, "-n"));
125  if (n < 0) n = 7;
126  } else n = 7;
127 
128 
129  {
130  World dw(argc, argv);
131 
132  if (rank == 0){
133  printf("Running Jacobi method on random %d-by-%d sparse matrix\n",n,n);
134  }
135  pass = jacobi(n, dw);
136  assert(pass);
137  }
138 
139  MPI_Finalize();
140  return 0;
141 }
147 #endif
Matrix class which encapsulates a 2D tensor.
Definition: matrix.h:18
int jacobi(int n, World &dw)
Definition: jacobi.cxx:19
def rank(self)
Definition: core.pyx:312
Definition: common.h:37
Vector class which encapsulates a 1D tensor.
Definition: vector.h:14
an instance of the CTF library (world) on a MPI communicator
Definition: world.h:19
string
Definition: core.pyx:456
dtype norm2()
computes the frobenius norm of the tensor (needs sqrt()!)
Definition: tensor.h:811
void fill_random(dtype rmin, dtype rmax)
fills local unique tensor elements to random values in the range [min,max] works only for dtype in {f...
Definition: tensor.cxx:928
int rank
rank of local processor
Definition: world.h:24
int main(int argc, char **argv)
Definition: jacobi.cxx:114
void sparsify()
reduce tensor to sparse format, storing only nonzero data, or data above a specified threshold...
Definition: tensor.cxx:449
Definition: apsp.cxx:17
char * getCmdOption(char **begin, char **end, const std::string &option)
Definition: jacobi.cxx:103
void jacobi_iter(Matrix<> &R, Vector<> &b, Vector<> &d, Vector<> &x)
Definition: jacobi.cxx:13
def np(self)
Definition: core.pyx:315