| 228 | **********************************************************************/ |
| 229 | |
| 230 | int QRsol(int n, real **h, real *q, real *b) |
| 231 | { |
| 232 | real c, s, temp1, temp2; |
| 233 | int i, k, q_ptr, code=0; |
| 234 | |
| 235 | /* Compute Q*b. */ |
| 236 | |
| 237 | for (k=0; k < n; k++) { |
| 238 | q_ptr = 2*k; |
| 239 | c = q[q_ptr]; |
| 240 | s = q[q_ptr+1]; |
| 241 | temp1 = b[k]; |
| 242 | temp2 = b[k+1]; |
| 243 | b[k] = c*temp1 - s*temp2; |
| 244 | b[k+1] = s*temp1 + c*temp2; |
| 245 | } |
| 246 | |
| 247 | /* Solve R*x = Q*b. */ |
| 248 | |
| 249 | for (k=n-1; k >= 0; k--) { |
| 250 | if (h[k][k] == ZERO) { |
| 251 | code = k + 1; |
| 252 | break; |
| 253 | } |
| 254 | b[k] /= h[k][k]; |
| 255 | for (i=0; i < k; i++) b[i] -= b[k]*h[i][k]; |
| 256 | } |
| 257 | |
| 258 | return (code); |
| 259 | } |
| 260 | |
| 261 | } |
| 262 | |