32 #ifndef _mpqc_src_lib_math_optimize_conjgrad_h
33 #define _mpqc_src_lib_math_optimize_conjgrad_h
36 #include <math/scmat/abstract.h>
57 double convergence_target = -1.0) {
59 const size_t n = x.nrow() * x.ncol();
60 MPQC_ASSERT(n == (preconditioner.nrow() * preconditioner.ncol()));
78 preconditioner.element_op(findmin_op);
79 preconditioner.element_op(findmax_op);
80 const double precond_min = findmin_op->result().at(0).value;
81 const double precond_max = findmax_op->result().at(0).value;
82 const double cond_number = precond_max / precond_min;
84 if (convergence_target < 0.0) {
85 convergence_target = 1e-15 * cond_number;
88 if (convergence_target < 1e-15 * cond_number)
89 ExEnv::out0() << indent <<
"WARNING: linsolv_conjugate_gradient conv. target (" << convergence_target
90 <<
") may be too low for 64-bit precision" << std::endl;
93 bool converged =
false;
94 const unsigned int max_niter = 500;
96 const unsigned int rhs_size = b.nrow() * b.ncol();
101 XX_i.element_op(multiply_op, preconditioner);
110 ZZ_i.element_op(multiply_op, preconditioner);
117 unsigned int iter = 0;
118 while (not converged) {
121 scalarprod_op->init();
122 RR_i.element_op(scalarprod_op, ZZ_i);
123 const double rz_norm2 = scalarprod_op->result();
126 scalarprod_op->init();
127 PP_i.element_op(scalarprod_op, APP_i);
128 const double pAp_i = scalarprod_op->result();
129 const double alpha_i = rz_norm2 / pAp_i;
133 XX_i.element_op(daxpy_op, PP_i);
137 RR_i.element_op(daxpy_op2, APP_i);
140 RR_i.element_op(norm2_op);
141 const double r_ip1_norm = norm2_op->result() / rhs_size;
142 if (r_ip1_norm < convergence_target) {
149 ZZ_i.element_op(multiply_op, preconditioner);
150 scalarprod_op->init();
151 ZZ_i.element_op(scalarprod_op, RR_i);
152 const double rz_ip1_norm2 = scalarprod_op->result();
154 const double beta_i = rz_ip1_norm2 / rz_norm2;
161 PP_i.element_op(daxpy_op3, ZZ_i);
166 if (iter >= max_niter) {
168 throw MaxIterExceeded(
"linsolv_conjugate_gradient", __FILE__, __LINE__, max_niter);
177 inline size_t size(
const RefSCMatrix& m) {
178 return m.nrow() * m.ncol();
181 inline RefSCMatrix clone(
const RefSCMatrix& m) {
185 inline RefSCMatrix copy(
const RefSCMatrix& m) {
189 inline double minabs_value(
const RefSCMatrix& m) {
190 typedef SCElementFindExtremum<sc::abs_greater<double>, SCMatrixIterationRanges::AllElements> MinOp;
191 Ref<MinOp> findmin_op =
new MinOp(std::numeric_limits<double>::max());
192 m.element_op(findmin_op);
193 return findmin_op->result().at(0).value;
196 inline double maxabs_value(
const RefSCMatrix& m) {
197 typedef SCElementFindExtremum<sc::abs_less<double>, SCMatrixIterationRanges::AllElements> MaxOp;
198 Ref<MaxOp> findmax_op =
new MaxOp;
199 m.element_op(findmax_op);
200 return findmax_op->result().at(0).value;
203 inline void vec_multiply(RefSCMatrix& m1,
const RefSCMatrix& m2) {
204 Ref<SCElementOp2> multiply_op =
new SCElementDestructiveProduct;
205 m1.element_op(multiply_op, m2);
208 inline double dot_product(
const RefSCMatrix& m1,
const RefSCMatrix& m2) {
209 Ref<SCElementScalarProduct> scalarprod_op =
new SCElementScalarProduct;
210 scalarprod_op->init();
211 m1.element_op(scalarprod_op, m2);
212 return scalarprod_op->result();
215 inline void scale(RefSCMatrix& m,
double scaling_factor) {
216 m.scale(scaling_factor);
219 inline void daxpy(RefSCMatrix& y,
double a,
const RefSCMatrix& x) {
220 Ref<SCElementDAXPY> daxpy_op =
new SCElementDAXPY(a);
221 y.element_op(daxpy_op, x);
224 inline void assign(RefSCMatrix& m1,
const RefSCMatrix& m2) {
228 inline double norm2(
const RefSCMatrix& m) {
229 Ref<SCElementKNorm> norm2_op =
new SCElementKNorm(2);
230 m.element_op(norm2_op);
231 return norm2_op->result();
234 inline void print(
const RefSCMatrix& m,
const char* label) {
263 template <
typename D,
typename F>
265 typedef typename D::element_type value_type;
267 value_type operator()(F& a,
270 const D& preconditioner,
271 value_type convergence_target = -1.0) {
273 const size_t n = size(x);
274 MPQC_ASSERT(n == size(preconditioner));
288 const value_type precond_min = minabs_value(preconditioner);
289 const value_type precond_max = maxabs_value(preconditioner);
290 const value_type cond_number = precond_max / precond_min;
293 if (convergence_target < 0.0) {
294 convergence_target = 1e-15 * cond_number;
297 if (convergence_target < 1e-15 * cond_number)
298 std::cout <<
"WARNING: ConjugateGradient convergence target (" << convergence_target
299 <<
") may be too low for 64-bit precision" << std::endl;
302 bool converged =
false;
303 const unsigned int max_niter = 500;
304 value_type rnorm2 = 0.0;
305 const size_t rhs_size = size(b);
310 vec_multiply(XX_i, preconditioner);
319 vec_multiply(ZZ_i, preconditioner);
324 unsigned int iter = 0;
325 while (not converged) {
328 value_type rz_norm2 = dot_product(RR_i, ZZ_i);
331 const value_type pAp_i = dot_product(PP_i, APP_i);
332 const value_type alpha_i = rz_norm2 / pAp_i;
335 daxpy(XX_i, alpha_i, PP_i);
338 daxpy(RR_i, -alpha_i, APP_i);
340 const value_type r_ip1_norm = norm2(RR_i) / rhs_size;
341 if (r_ip1_norm < convergence_target) {
348 vec_multiply(ZZ_i, preconditioner);
350 const value_type rz_ip1_norm2 = dot_product(ZZ_i, RR_i);
352 const value_type beta_i = rz_ip1_norm2 / rz_norm2;
358 daxpy(PP_i, 1.0, ZZ_i);
363 if (iter >= max_niter) {
365 throw std::domain_error(
"ConjugateGradient: max # of iterations exceeded");
377 #endif // end of header guard