| 82 | **********************************************************************/ |
| 83 | |
| 84 | int ClassicalGS(N_Vector *v, real **h, int k, int p, real *new_vk_norm, |
| 85 | N_Vector temp, real *s) |
| 86 | { |
| 87 | int i, k_minus_1, i0; |
| 88 | real vk_norm; |
| 89 | |
| 90 | k_minus_1 = k - 1; |
| 91 | |
| 92 | /* Perform Classical Gram-Schmidt */ |
| 93 | |
| 94 | vk_norm = RSqrt(N_VDotProd(v[k], v[k])); |
| 95 | |
| 96 | i0 = MAX(k-p, 0); |
| 97 | for (i=i0; i < k; i++) { |
| 98 | h[i][k_minus_1] = N_VDotProd(v[i], v[k]); |
| 99 | } |
| 100 | |
| 101 | for (i=i0; i < k; i++) { |
| 102 | N_VLinearSum(ONE, v[k], -h[i][k_minus_1], v[i], v[k]); |
| 103 | } |
| 104 | |
| 105 | /* Compute the norm of the new vector at v[k]. */ |
| 106 | |
| 107 | *new_vk_norm = RSqrt(N_VDotProd(v[k], v[k])); |
| 108 | |
| 109 | /* Reorthogonalize if necessary */ |
| 110 | |
| 111 | if ((FACTOR * (*new_vk_norm)) < vk_norm) { |
| 112 | |
| 113 | for (i=i0; i < k; i++) { |
| 114 | s[i] = N_VDotProd(v[i], v[k]); |
| 115 | } |
| 116 | |
| 117 | if (i0 < k) { |
| 118 | N_VScale(s[i0], v[i0], temp); |
| 119 | h[i0][k_minus_1] += s[i0]; |
| 120 | } |
| 121 | for (i=i0+1; i < k; i++) { |
| 122 | N_VLinearSum(s[i], v[i], ONE, temp, temp); |
| 123 | h[i][k_minus_1] += s[i]; |
| 124 | } |
| 125 | N_VLinearSum(ONE, v[k], -ONE, temp, v[k]); |
| 126 | |
| 127 | *new_vk_norm = RSqrt(N_VDotProd(v[k],v[k])); |
| 128 | } |
| 129 | |
| 130 | return(0); |
| 131 | } |
| 132 | |
| 133 | /*************** QRfact ********************************************** |
| 134 | This implementation of QRfact is a slight modification of a previous |
no test coverage detected