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