| 748 | } |
| 749 | |
| 750 | template<typename T> static void |
| 751 | SVBkSbImpl_( int m, int n, const T* w, int incw, |
| 752 | const T* u, int ldu, bool uT, |
| 753 | const T* v, int ldv, bool vT, |
| 754 | const T* b, int ldb, int nb, |
| 755 | T* x, int ldx, double* buffer, T eps ) |
| 756 | { |
| 757 | double threshold = 0; |
| 758 | int udelta0 = uT ? ldu : 1, udelta1 = uT ? 1 : ldu; |
| 759 | int vdelta0 = vT ? ldv : 1, vdelta1 = vT ? 1 : ldv; |
| 760 | int i, j, nm = std::min(m, n); |
| 761 | |
| 762 | if( !b ) |
| 763 | nb = m; |
| 764 | |
| 765 | for( i = 0; i < n; i++ ) |
| 766 | for( j = 0; j < nb; j++ ) |
| 767 | x[i*ldx + j] = 0; |
| 768 | |
| 769 | for( i = 0; i < nm; i++ ) |
| 770 | threshold += w[i*incw]; |
| 771 | threshold *= eps; |
| 772 | |
| 773 | // v * inv(w) * uT * b |
| 774 | for( i = 0; i < nm; i++, u += udelta0, v += vdelta0 ) |
| 775 | { |
| 776 | double wi = w[i*incw]; |
| 777 | if( (double)std::abs(wi) <= threshold ) |
| 778 | continue; |
| 779 | wi = 1/wi; |
| 780 | |
| 781 | if( nb == 1 ) |
| 782 | { |
| 783 | double s = 0; |
| 784 | if( b ) |
| 785 | for( j = 0; j < m; j++ ) |
| 786 | s += u[j*udelta1]*b[j*ldb]; |
| 787 | else |
| 788 | s = u[0]; |
| 789 | s *= wi; |
| 790 | |
| 791 | for( j = 0; j < n; j++ ) |
| 792 | x[j*ldx] = (T)(x[j*ldx] + s*v[j*vdelta1]); |
| 793 | } |
| 794 | else |
| 795 | { |
| 796 | if( b ) |
| 797 | { |
| 798 | for( j = 0; j < nb; j++ ) |
| 799 | buffer[j] = 0; |
| 800 | MatrAXPY( m, nb, b, ldb, u, udelta1, buffer, 0 ); |
| 801 | for( j = 0; j < nb; j++ ) |
| 802 | buffer[j] *= wi; |
| 803 | } |
| 804 | else |
| 805 | { |
| 806 | for( j = 0; j < nb; j++ ) |
| 807 | buffer[j] = u[j*udelta1]*wi; |