23dfcca431 2011-02-23 kinaba: //------------------------------------------------------------- 23dfcca431 2011-02-23 kinaba: // Linear Equation Solver for Teplitz Matrix 23dfcca431 2011-02-23 kinaba: // O(n^2) 23dfcca431 2011-02-23 kinaba: // 23dfcca431 2011-02-23 kinaba: // Faster solver for matrices satisfying 23dfcca431 2011-02-23 kinaba: // M[i][j] == M[i+1][j+1] 23dfcca431 2011-02-23 kinaba: // for all i, j>0 23dfcca431 2011-02-23 kinaba: // Be careful that this implementation does not support the 23dfcca431 2011-02-23 kinaba: // case M[0][0]=M[1][1]=...=0. (zero-division) 23dfcca431 2011-02-23 kinaba: // 23dfcca431 2011-02-23 kinaba: // TODO: support M[0][0]=M[1][1]=...=0 23dfcca431 2011-02-23 kinaba: // QUICKHACK: 23dfcca431 2011-02-23 kinaba: // find the least degenerated k and 23dfcca431 2011-02-23 kinaba: // calc initial vectors by gauss_jordan... 23dfcca431 2011-02-23 kinaba: // 23dfcca431 2011-02-23 kinaba: // Verified by 23dfcca431 2011-02-23 kinaba: // - (Checked by random data to match with linearEq.cpp) 23dfcca431 2011-02-23 kinaba: //------------------------------------------------------------- 23dfcca431 2011-02-23 kinaba: 23dfcca431 2011-02-23 kinaba: vector<double> solve_teplitz( int n, vector< vector<double> > M, const vector<double>& V ) 23dfcca431 2011-02-23 kinaba: { 23dfcca431 2011-02-23 kinaba: vector<double> f, b, x; 23dfcca431 2011-02-23 kinaba: 23dfcca431 2011-02-23 kinaba: // invariants 23dfcca431 2011-02-23 kinaba: // M|k f|k = (1 0 ... 0) 23dfcca431 2011-02-23 kinaba: // M|k b|k = (0 ... 0 1) 23dfcca431 2011-02-23 kinaba: // M|k x|k = V|k 23dfcca431 2011-02-23 kinaba: f.push_back( 1 / M[0][0] ); 23dfcca431 2011-02-23 kinaba: b.push_back( 1 / M[0][0] ); 23dfcca431 2011-02-23 kinaba: x.push_back( V[0] / M[0][0] ); 23dfcca431 2011-02-23 kinaba: 23dfcca431 2011-02-23 kinaba: for(int k=2; k<=n; ++k) 23dfcca431 2011-02-23 kinaba: { 23dfcca431 2011-02-23 kinaba: // M|k (f|k-1 0) = (1 0 ... 0 ea) 23dfcca431 2011-02-23 kinaba: double ea = 0; 23dfcca431 2011-02-23 kinaba: for(int i=0; i<k-1; ++i) 23dfcca431 2011-02-23 kinaba: ea += M[k-1][i] * f[i]; 23dfcca431 2011-02-23 kinaba: 23dfcca431 2011-02-23 kinaba: // M|k (0 b|k-1) = (eb 0 ... 0 1) 23dfcca431 2011-02-23 kinaba: double eb = 0; 23dfcca431 2011-02-23 kinaba: for(int i=1; i<k; ++i) 23dfcca431 2011-02-23 kinaba: eb += M[0][i] * b[i-1]; 23dfcca431 2011-02-23 kinaba: 23dfcca431 2011-02-23 kinaba: // u(1 ea) + v(eb 1) = (1 0) 23dfcca431 2011-02-23 kinaba: //==> u = 1 / (1-ea.eb), v = -ea /(1-ea.eb) 23dfcca431 2011-02-23 kinaba: double u = 1/(1-ea*eb), v = -ea/(1-ea*eb); 23dfcca431 2011-02-23 kinaba: vector<double> ff; 23dfcca431 2011-02-23 kinaba: for(int i=0; i<k; ++i) 23dfcca431 2011-02-23 kinaba: ff.push_back( u*(i<k-1 ? f[i] : 0) + v*(i>0 ? b[i-1] : 0) ); 23dfcca431 2011-02-23 kinaba: 23dfcca431 2011-02-23 kinaba: // similarly... 23dfcca431 2011-02-23 kinaba: u = -eb/(1-ea*eb), v = 1/(1-ea*eb); 23dfcca431 2011-02-23 kinaba: vector<double> bb; 23dfcca431 2011-02-23 kinaba: for(int i=0; i<k; ++i) 23dfcca431 2011-02-23 kinaba: bb.push_back( u*(i<k-1 ? f[i] : 0) + v*(i>0 ? b[i-1] : 0) ); 23dfcca431 2011-02-23 kinaba: 23dfcca431 2011-02-23 kinaba: // update 23dfcca431 2011-02-23 kinaba: f.swap(ff); 23dfcca431 2011-02-23 kinaba: b.swap(bb); 23dfcca431 2011-02-23 kinaba: 23dfcca431 2011-02-23 kinaba: // M|k (x|k 0) (V|k ec) 23dfcca431 2011-02-23 kinaba: double ec = 0; 23dfcca431 2011-02-23 kinaba: for(int i=0; i<k-1; ++i) 23dfcca431 2011-02-23 kinaba: ec += M[k-1][i] * x[i]; 23dfcca431 2011-02-23 kinaba: x.push_back(0); 23dfcca431 2011-02-23 kinaba: for(int i=0; i<k; ++i) 23dfcca431 2011-02-23 kinaba: x[i] -= b[i] * (ec-V[k-1]); 23dfcca431 2011-02-23 kinaba: } 23dfcca431 2011-02-23 kinaba: 23dfcca431 2011-02-23 kinaba: return x; 23dfcca431 2011-02-23 kinaba: }