| 174 | |
| 175 | template <class T> |
| 176 | Matrix<T> MultiVector::T_Orthogonalize (BaseMatrix * ipmat) |
| 177 | { |
| 178 | static Timer t("MultiVector::Orthogonalize"); |
| 179 | RegionTimer reg(t); |
| 180 | auto & mv = *this; |
| 181 | Matrix<T> Rfactor(mv.Size()); |
| 182 | Rfactor = T(0.0); |
| 183 | if (ipmat) |
| 184 | { |
| 185 | auto tmp = mv.RefVec()->CreateVector(); |
| 186 | for (int i = 0; i < mv.Size(); i++) |
| 187 | { |
| 188 | *tmp = *ipmat * *mv[i]; |
| 189 | double norm = sqrt(fabs(InnerProduct<T>(*tmp, *mv[i], true))); |
| 190 | Rfactor(i,i) = norm; |
| 191 | *mv[i] *= 1.0 / norm; |
| 192 | /* |
| 193 | for (int j = i+1; j < mv.Size(); j++) |
| 194 | { |
| 195 | T rij = InnerProduct<T>(*mv[j], *tmp, true) / norm; |
| 196 | Rfactor(i,j) = rij; |
| 197 | *mv[j] -= rij * *mv[i]; |
| 198 | } |
| 199 | */ |
| 200 | IntRange rest(i+1, mv.Size()); |
| 201 | auto mvrest = mv.Range(rest); |
| 202 | Rfactor.Row(i).Range(rest) = (1.0/norm) * mvrest->T_InnerProduct<T> (*tmp, true); |
| 203 | for (auto j : rest) |
| 204 | *mv[j] -= Rfactor(i,j) * *mv[i]; |
| 205 | } |
| 206 | } |
| 207 | else |
| 208 | { |
| 209 | if (mv.Size() == 1) |
| 210 | { |
| 211 | double norm = (*mv[0]).L2Norm(); |
| 212 | Rfactor(0,0) = norm; |
| 213 | *mv[0] *= 1.0 / norm; |
| 214 | } |
| 215 | else |
| 216 | { |
| 217 | int half = mv.Size()/2; |
| 218 | auto r1 = IntRange(0,half); |
| 219 | auto r2 = IntRange(half, mv.Size()); |
| 220 | auto mv1 = mv.Range(r1); |
| 221 | auto mv2 = mv.Range(r2); |
| 222 | Rfactor.Rows(r1).Cols(r1) = mv1->T_Orthogonalize<T>(ipmat); |
| 223 | Matrix<T> ip = Conj(InnerProduct<T> (*mv1, *mv2, true)); |
| 224 | Rfactor.Rows(r2).Cols(r1) = ip; |
| 225 | ip *= -1; |
| 226 | mv2->Add (*mv1, ip); |
| 227 | Rfactor.Rows(r2).Cols(r2) = mv2->T_Orthogonalize<T>(ipmat); |
| 228 | } |
| 229 | } |
| 230 | return Rfactor; |
| 231 | } |
| 232 | |
| 233 | void MultiVector :: Orthogonalize (BaseMatrix * ipmat) |