| 95 | } |
| 96 | |
| 97 | private int pcg(double[] g, double[] M, double[] s, double[] r) { |
| 98 | int i, inc = 1; |
| 99 | int n = fun_obj.get_nr_variable(); |
| 100 | double one = 1; |
| 101 | double[] d = new double[n]; |
| 102 | double[] Hd = new double[n]; |
| 103 | double zTr, znewTrnew, alpha, beta, cgtol, dHd; |
| 104 | double[] z = new double[n]; |
| 105 | double Q = 0, newQ, Qdiff; |
| 106 | |
| 107 | for (i = 0; i < n; i++) { |
| 108 | s[i] = 0; |
| 109 | r[i] = -g[i]; |
| 110 | z[i] = r[i] / M[i]; |
| 111 | d[i] = z[i]; |
| 112 | } |
| 113 | |
| 114 | zTr = Blas.ddot_(n, z, inc, r, inc); |
| 115 | double gMinv_norm = Math.sqrt(zTr); |
| 116 | cgtol = Math.min(eps_cg, Math.sqrt(gMinv_norm)); |
| 117 | int cg_iter = 0; |
| 118 | int max_cg_iter = Math.max(n, 5); |
| 119 | |
| 120 | while (cg_iter < max_cg_iter) { |
| 121 | cg_iter++; |
| 122 | |
| 123 | fun_obj.Hv(d, Hd); |
| 124 | dHd = Blas.ddot_(n, d, inc, Hd, inc); |
| 125 | // avoid 0/0 in getting alpha |
| 126 | if (dHd <= 1.0e-16) |
| 127 | break; |
| 128 | |
| 129 | alpha = zTr / dHd; |
| 130 | Blas.daxpy_(n, alpha, d, inc, s, inc); |
| 131 | alpha = -alpha; |
| 132 | Blas.daxpy_(n, alpha, Hd, inc, r, inc); |
| 133 | |
| 134 | // Using quadratic approximation as CG stopping criterion |
| 135 | newQ = -0.5 * (Blas.ddot_(n, s, inc, r, inc) - Blas.ddot_(n, s, inc, g, inc)); |
| 136 | Qdiff = newQ - Q; |
| 137 | if (newQ <= 0 && Qdiff <= 0) { |
| 138 | if (cg_iter * Qdiff >= cgtol * newQ) |
| 139 | break; |
| 140 | } else { |
| 141 | info("WARNING: quadratic approximation > 0 or increasing in CG%n"); |
| 142 | break; |
| 143 | } |
| 144 | Q = newQ; |
| 145 | |
| 146 | for (i = 0; i < n; i++) |
| 147 | z[i] = r[i] / M[i]; |
| 148 | znewTrnew = Blas.ddot_(n, z, inc, r, inc); |
| 149 | beta = znewTrnew / zTr; |
| 150 | Blas.dscal_(n, beta, d, inc); |
| 151 | Blas.daxpy_(n, one, z, inc, d, inc); |
| 152 | zTr = znewTrnew; |
| 153 | } |
| 154 | |