Cyclops Tensor Framework
parallel arithmetic on multidimensional arrays
qr.cxx
Go to the documentation of this file.
1 
8 #include <ctf.hpp>
9 #include "conj.h"
10 using namespace CTF;
11 
12 
13 template <typename dtype>
15  int m,
16  int n,
17  World & dw){
18 
19  // Perform QR
20  Matrix<dtype> Q,R;
21  A.qr(Q,R);
22 
23  // Test orthogonality
24  Matrix<dtype> E(n,n,dw);
25 
26  E["ii"] = 1.;
27 
28  E["ij"] -= Q["ki"]*conj<dtype>(Q)["kj"];
29 
30  bool pass_orthogonality = true;
31 
32  double nrm;
33  E.norm2(nrm);
34  if (nrm > m*n*1.E-6){
35  pass_orthogonality = false;
36  }
37 
38  A["ij"] -= Q["ik"]*R["kj"];
39 
40  bool pass_residual = true;
41  A.norm2(nrm);
42  if (nrm > m*n*n*1.E-6){
43  pass_residual = false;
44  }
45 
46 #ifndef TEST_SUITE
47  if (dw.rank == 0){
48  printf("QR orthogonality check returned %d, residual check %d\n", pass_orthogonality, pass_residual);
49  }
50 #endif
51  return pass_residual & pass_orthogonality;
52 }
53 
54 bool test_qr(int m, int n, World dw){
55  bool pass = true;
56  Matrix<float> A(m,n,dw);
57  Matrix<float> AA(m,n,dw);
58  A.fill_random(0.,1.);
59  AA.fill_random(0.,1.);
60  pass = pass & qr<float>(A,m,n,dw);
61 
62  Matrix<double> B(m,n,dw);
63  Matrix<double> BB(m,n,dw);
64  B.fill_random(0.,1.);
65  BB.fill_random(0.,1.);
66  pass = pass & qr<double>(B,m,n,dw);
67 
68  Matrix<std::complex<float>> cA(m,n,dw);
69  cA["ij"] = Function<float,float,std::complex<float>>([](float a, float b){ return std::complex<float>(a,b); })(A["ij"],AA["ij"]);
70  pass = pass & qr<std::complex<float>>(cA,m,n,dw);
71 
72 
73  Matrix<std::complex<double>> cB(m,n,dw);
74  cB["ij"] = Function<double,double,std::complex<double>>([](double a, double b){ return std::complex<double>(a,b); })(B["ij"],BB["ij"]);
75  pass = pass & qr<std::complex<double>>(cB,m,n,dw);
76  if (dw.rank == 0){
77  if (pass){
78  printf("{ A = QR and Q^TQ = I } passed\n");
79  } else {
80  printf("{ A = QR and Q^TQ = I } failed\n");
81  }
82  }
83 
84 
85  return pass;
86 }
87 
88 #ifndef TEST_SUITE
89 char* getCmdOption(char ** begin,
90  char ** end,
91  const std::string & option){
92  char ** itr = std::find(begin, end, option);
93  if (itr != end && ++itr != end){
94  return *itr;
95  }
96  return 0;
97 }
98 
99 
100 int main(int argc, char ** argv){
101  int rank, np, m, n, pass;
102  int const in_num = argc;
103  char ** input_str = argv;
104 
105  MPI_Init(&argc, &argv);
106  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
107  MPI_Comm_size(MPI_COMM_WORLD, &np);
108 
109  if (getCmdOption(input_str, input_str+in_num, "-m")){
110  m = atoi(getCmdOption(input_str, input_str+in_num, "-m"));
111  if (m < 0) m = 13;
112  } else m = 13;
113 
114 
115  if (getCmdOption(input_str, input_str+in_num, "-n")){
116  n = atoi(getCmdOption(input_str, input_str+in_num, "-n"));
117  if (n < 0) n = 7;
118  } else n = 7;
119 
120 
121  {
122  World dw(argc, argv);
123 
124  if (rank == 0){
125  printf("Testing %d-by-%d QR factorization\n", m, n);
126  }
127  pass = test_qr(m, n, dw);
128  assert(pass);
129  }
130 
131  MPI_Finalize();
132  return 0;
133 }
139 #endif
bool qr(Matrix< dtype > A, int m, int n, World &dw)
Definition: qr.cxx:14
void qr(Matrix< dtype > &Q, Matrix< dtype > &R)
Definition: matrix.cxx:425
Matrix class which encapsulates a 2D tensor.
Definition: matrix.h:18
def rank(self)
Definition: core.pyx:312
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: qr.cxx:100
Definition: apsp.cxx:17
char * getCmdOption(char **begin, char **end, const std::string &option)
Definition: qr.cxx:89
bool test_qr(int m, int n, World dw)
Definition: qr.cxx:54
def np(self)
Definition: core.pyx:315