| 513 | |
| 514 | |
| 515 | void MultTriangularUR3 (SliceMatrix<double> X, BareSliceMatrix<double> T) |
| 516 | { |
| 517 | size_t n = X.Width(); |
| 518 | size_t m = X.Height(); |
| 519 | |
| 520 | /* |
| 521 | for (size_t i = n; i-->0; ) |
| 522 | { |
| 523 | X.Col(i) *= T(i,i); |
| 524 | for (size_t j = 0; j < i; j++) |
| 525 | X.Col(i) += T(j,i) * X.Col(j); |
| 526 | } |
| 527 | */ |
| 528 | |
| 529 | constexpr size_t BS = 2*SIMD<double>::Size(); |
| 530 | alignas (SIMD<double>) double memb[BS*128]; |
| 531 | size_t i = n; |
| 532 | |
| 533 | for ( ; i >= BS; i -= BS) |
| 534 | { |
| 535 | // copy b |
| 536 | for (size_t j = 0; j < i; j++) |
| 537 | for (size_t k = 0; k < BS; k++) |
| 538 | memb[BS*j+k] = T(j,i+k-BS); |
| 539 | for (int k = 0; k < BS; k++) |
| 540 | for (int j = 0; j < k; j++) |
| 541 | memb[BS*(i-BS+k)+j] = 0; |
| 542 | |
| 543 | /* |
| 544 | FlatMatrix<> subb(i, BS, &memb[0]); |
| 545 | for (size_t j = 0; j < m; j++) |
| 546 | { |
| 547 | Vec<BS> c = Trans(subb) * X.Row(j).Range(0,i); |
| 548 | X.Row(j).Range(i-BS,i) = c; |
| 549 | } |
| 550 | */ |
| 551 | size_t j = 0; |
| 552 | for ( ; j+4 <= m; j+=4) |
| 553 | MatKernelMultAB<4,2,SET> (i, &X(j,0), X.Dist(), &memb[0], BS, &X(j,i-BS), X.Dist()); |
| 554 | for ( ; j+1 <= m; j+=1) |
| 555 | MatKernelMultAB<1,2,SET> (i, &X(j,0), X.Dist(), &memb[0], BS, &X(j,i-BS), X.Dist()); |
| 556 | } |
| 557 | |
| 558 | for ( ; i >= 1; i--) |
| 559 | { |
| 560 | X.Col(i-1) *= T(i-1, i-1); |
| 561 | if (i > 1) |
| 562 | X.Col(i-1) += X.Cols(0, i-1) * T.Col(i-1).Range(0, i-1); |
| 563 | } |
| 564 | } |
| 565 | |
| 566 | void MultTriangularUR2 (SliceMatrix<double> X, BareSliceMatrix<double> T) |
| 567 | { |