| 136 | **********************************************************************/ |
| 137 | |
| 138 | int QRfact(int n, real **h, real *q, int job) |
| 139 | { |
| 140 | real c, s, temp1, temp2, temp3; |
| 141 | int i, j, k, q_ptr, n_minus_1, code=0; |
| 142 | |
| 143 | switch (job) { |
| 144 | case 0: |
| 145 | /* Compute a new factorization of H. */ |
| 146 | code = 0; |
| 147 | for (k=0; k < n; k++) { |
| 148 | |
| 149 | /* Multiply column k by the previous k-1 Givens rotations. */ |
| 150 | for (j=0; j < k-1; j++) { |
| 151 | i = 2*j; |
| 152 | temp1 = h[j][k]; |
| 153 | temp2 = h[j+1][k]; |
| 154 | c = q[i]; |
| 155 | s = q[i+1]; |
| 156 | h[j][k] = c*temp1 - s*temp2; |
| 157 | h[j+1][k] = s*temp1 + c*temp2; |
| 158 | } |
| 159 | |
| 160 | /* Compute the Givens rotation components c and s */ |
| 161 | q_ptr = 2*k; |
| 162 | temp1 = h[k][k]; |
| 163 | temp2 = h[k+1][k]; |
| 164 | if( temp2 == ZERO) { |
| 165 | c = ONE; |
| 166 | s = ZERO; |
| 167 | } else if (ABS(temp2) >= ABS(temp1)) { |
| 168 | temp3 = temp1/temp2; |
| 169 | s = -ONE/RSqrt(ONE+SQR(temp3)); |
| 170 | c = -s*temp3; |
| 171 | } else { |
| 172 | temp3 = temp2/temp1; |
| 173 | c = ONE/RSqrt(ONE+SQR(temp3)); |
| 174 | s = -c*temp3; |
| 175 | } |
| 176 | q[q_ptr] = c; |
| 177 | q[q_ptr+1] = s; |
| 178 | if( (h[k][k] = c*temp1 - s*temp2) == ZERO) code = k+1; |
| 179 | } |
| 180 | break; |
| 181 | |
| 182 | default: |
| 183 | /* Update the factored H to which a new column has been added. */ |
| 184 | n_minus_1 = n - 1; |
| 185 | code = 0; |
| 186 | |
| 187 | /* Multiply the new column by the previous n-1 Givens rotations. */ |
| 188 | for (k=0; k < n_minus_1; k++) { |
| 189 | i = 2*k; |
| 190 | temp1 = h[k][n_minus_1]; |
| 191 | temp2 = h[k+1][n_minus_1]; |
| 192 | c = q[i]; |
| 193 | s = q[i+1]; |
| 194 | h[k][n_minus_1] = c*temp1 - s*temp2; |
| 195 | h[k+1][n_minus_1] = s*temp1 + c*temp2; |