SIRF  3.5.0
JacobiCG.h
1 /*
2 SyneRBI Synergistic Image Reconstruction Framework (SIRF)
3 Copyright 2020 Rutherford Appleton Laboratory STFC
4 
5 This is software developed for the Collaborative Computational
6 Project in Synergistic Reconstruction for Biomedical Imaging (formerly CCP PETMR)
7 (http://www.ccpsynerbi.ac.uk/).
8 
9 Licensed under the Apache License, Version 2.0 (the "License");
10 you may not use this file except in compliance with the License.
11 You may obtain a copy of the License at
12 http://www.apache.org/licenses/LICENSE-2.0
13 Unless required by applicable law or agreed to in writing, software
14 distributed under the License is distributed on an "AS IS" BASIS,
15 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16 See the License for the specific language governing permissions and
17 limitations under the License.
18 
19 */
20 
21 #pragma once
22 
23 #include <cmath>
24 #include <complex>
25 #include <iostream>
26 #include <memory>
27 
28 #include "sirf/common/Operator.h"
29 
30 /* A simplified implementation of Jacobi-Conjugated Gradients method
31  for computing extreme eigenvalues of real symmetric and Hermitian
32  operators by E. Ovtchinnikov [SIAM Numer. Anal. 46:2567-2619, 2008].
33 */
34 
35 namespace sirf {
36  template<class value_type>
37  class JacobiCG {
38  public:
39  /* This simplified version does fixed number of iterations. */
40  JacobiCG() : nit_(10) {}
41  void set_num_iterations(int nit)
42  {
43  nit_ = nit;
44  }
45 
46  /* Computes the largest eigenvalue of a positive semi-definite operator */
47  template<class vector_type>
48  float largest(Operator<vector_type>& A,
49  vector_type& x /*non-zero initial guess*/, int verb=0)
50  {
51  value_type lmd;
52  value_type a[] = { 0, 0, 0, 0 };
53  value_type mu[2];
54  value_type u[2];
55  value_type v[2];
56  value_type zero = 0;
57  value_type t;
58 
59  std::unique_ptr<vector_type> uptr_y = x.clone();
60  vector_type& y = *uptr_y;
61  std::unique_ptr<vector_type> uptr_z = x.clone();
62  vector_type& z = *uptr_z;
63  std::unique_ptr<vector_type> uptr_w = x.clone();
64  vector_type& w = *uptr_w;
65  std::unique_ptr<vector_type> uptr_Az = x.clone();
66  vector_type& Az = *uptr_Az;
67 
68  float s = x.norm();
69  x.scale(s);
70  std::shared_ptr<vector_type> sptr_Ax = A(x);
71  vector_type& Ax = *sptr_Ax;
72 
73  for (int it = 0; it < nit_; it++) {
74  lmd = Ax.dot(x);
75  y.axpby(1.0, Ax, -lmd, x); // residual
76  if (verb > 0) {
77  if (it == 0)
78  std::cout << '\n' << std::flush;
79  std::cout << "CG iteration " << it
80  << ": largest eigenvalue " << std::real(lmd) << '\n';
81  }
82  if (it) { // conjugate y to the previous search direction
83  t = Az.dot(z);
84  w.axpby(1.0, Az, -t, z);
85  t = w.dot(y) / (lmd - t);
86  y.axpby(1.0, y, t, z);
87  }
88  // normalize y
89  s = y.norm();
90  if (s == 0.0)
91  break; // converged
92  y.scale(s);
93  // orthogonalize y to x
94  t = y.dot(x);
95  y.axpby(1.0, y, -t, x);
96  // normalize y again
97  s = y.norm();
98  if (s == 0.0)
99  break; // converged
100  y.scale(s);
101  // perform Rayleigh-Ritz procedure in span{x,y}
102  std::shared_ptr<vector_type> sptr_Ay = A(y);
103  vector_type& Ay = *sptr_Ay;
104  a[0] = lmd;
105  a[1] = Ay.dot(x);
106  a[2] = x.dot(Ay);
107  a[3] = Ay.dot(y);
108  // compute eigenvalues and eigenvectors of 2x2 matrix a
109  eigh2_(a, mu, u, v);
110  z.axpby(u[0], x, u[1], y);
111  x.axpby(v[0], x, v[1], y);
112  // span{x,y} = span{x,z} => the new x is a linear combination
113  // of the old x and z => we use z as previous search direction
114  // on the next iteration
115  Az.axpby(u[0], Ax, u[1], Ay);
116  Ax.axpby(v[0], Ax, v[1], Ay);
117  lmd = mu[1];
118  s = x.norm();
119  x.scale(s);
120  Ax.scale(s);
121  }
122  return std::real(lmd);
123  }
124  private:
125  int nit_;
126  /* computes eigenvalues and eigenvectors of a 2x2 real symmetric
127  or Hermitian matrix [ [a[0], a[1]], [a[2], a[3]] ]
128  */
129  void eigh2_(const value_type a[4], value_type lmd[2], value_type x[2],
130  value_type y[2]) const
131  {
132  value_type zero = 0;
133  value_type one = 1;
134  value_type c = abs(a[1]);
135  value_type d = abs(a[0] - a[3]);
136  d = d * d + (c + c) * (c + c);
137  value_type s = sqrt(d);
138  value_type t = a[0] + a[3];
139  if (s == zero && t == zero) {
140  lmd[0] = 0.0;
141  lmd[1] = 0.0;
142  x[0] = 1.0;
143  x[1] = 0.0;
144  x[2] = 0.0;
145  x[3] = 1.0;
146  return;
147  }
148  lmd[0] = (t - s) / (one + one);
149  lmd[1] = (t + s) / (one + one);
150  value_type p = a[1];
151  value_type q = lmd[0] - a[0];
152  value_type pc = (p == zero ? 0 : abs(p*p) / p);
153  value_type qc = (q == zero ? 0 : abs(q*q) / q);
154  x[0] = p;
155  x[1] = q;
156  s = sqrt(abs(pc*x[0] + qc*x[1]));
157  x[0] /= s;
158  x[1] /= s;
159  y[0] = -qc / s;
160  y[1] = pc / s;
161  return;
162  }
163  };
164 
165  template<typename vector_type, typename value_type>
166  class Wrapped_sptr {
167  public:
168  Wrapped_sptr(std::shared_ptr<vector_type> sptr) : sptr_(sptr) {}
169  std::shared_ptr<vector_type> sptr() const
170  {
171  return sptr_;
172  }
173  std::unique_ptr<Wrapped_sptr<vector_type, value_type> > clone() const
174  {
175  std::shared_ptr<vector_type> sptr(sptr_->clone());
176  return std::unique_ptr<Wrapped_sptr>(new Wrapped_sptr(sptr));
177  }
178  float norm() const
179  {
180  return sptr_->norm();
181  }
182  void scale(float s)
183  {
184  sptr_->scale(s);
185  }
186  value_type dot(const Wrapped_sptr& y)
187  {
188  value_type s;
189  void* ptr = (void*)&s;
190  sptr_->dot(*y.sptr(), ptr);
191  return s;
192  }
193  void axpby(value_type a, const Wrapped_sptr& x, value_type b, const Wrapped_sptr& y)
194  {
195  void* ptr_a = (void*)&a;
196  void* ptr_b = (void*)&b;
197  sptr_->axpby(ptr_a, *x.sptr(), ptr_b, *y.sptr());
198  }
199  protected:
200  std::shared_ptr<vector_type> sptr_;
201  };
202 }
Definition: JacobiCG.h:37
Definition: Operator.h:7
Definition: JacobiCG.h:166
Abstract data container.
Definition: GeometricalInfo.cpp:141