1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 | //matrix and vector multiplication when matrix is the upper triangular matrix. inline mxd::t m_mulUv(mxd::t _m, vxd::t _v){ vxd::t v_(_m.rows()); v_.setZero(); i::t idx = 0; while(idx < _m.rows()){ mxd::b lTop = _m.block( 0 , idx, idx + 1 , 1 ); vxd::b yTop = v_.segment(0 , idx +1 ); yTop = v_axpy(_v[idx], lTop, yTop); idx++; } return v_; } //matrix and vector multiplication when matrix is the lower triangular matrix. inline mxd::t m_mulLv(mxd::t _m, vxd::t _v){ vxd::t v_(_m.rows()); v_.setZero(); i::t idx = 0; while(idx < _m.rows()){ mxd::b lBot = _m.block( idx , idx, _m.rows() - idx, 1 ); vxd::b yBot = v_.segment(idx , _v.rows() - idx ); yBot = v_axpy(_v[idx], lBot, yBot); idx++; } return v_; } | cs |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 | //matrix and vector multiplication when matrix is the lower semetric matrix. inline vxd::t v_mulSem(mxd::t _m , vxd::t _v){ i::t idx = 0; vxd::t v_(_m.rows()); v_.setZero(); while(idx < _m.rows()){ mxd::b A00 = _m.block(0 , 0 , idx , idx ); mxd::b a10 = _m.block(idx , 0 , 1 , idx ); mxd::b A20 = _m.block(idx + 1, 0 , _m .rows() - (idx + 1) , idx ); mxd::b a01 = _m.block(0 , idx , idx , 1 ); mxd::b alpha11 = _m.block(idx , idx , 1 , 1 ); mxd::b a21 = _m.block(idx + 1, idx , _m.rows() - (idx + 1) , 1 ); mxd::b A02 = _m.block(0 , idx + 1 , idx , _m.cols() - (idx + 1) ); mxd::b a12 = _m.block(idx , idx + 1 , 1 , _m.cols() - (idx + 1) ); mxd::b A22 = _m.block(idx + 1, idx + 1 , _m.rows() - (idx + 1) , _m.cols() - (idx + 1) ); vxd::b y0 = v_.segment(0 , idx ); vxd::b phi1 = v_.segment(idx , 1 ); vxd::b y2 = v_.segment(idx + 1 , v_.rows() - (idx + 1) ); y0 = v_axpy(_v[idx],a01, y0 ); phi1 = _v[idx] * alpha11 + phi1; y2 = v_axpy(_v[idx],a12.transpose(), y2 ); idx++; } return v_; } | cs |