/* ===================================================================
 * Copyright(C) LinBox 2008
 * Written by Jean-Guillaume Dumas
 * Triangular Solve
 *
 * ========LICENCE========
 * This file is part of the library LinBox.
 *
 * LinBox is free software: you can redistribute it and/or modify
 * it under the terms of the  GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 * ========LICENCE========
 *.
 * Time-stamp: <01 Oct 09 15:38:25 Jean-Guillaume.Dumas@imag.fr>
 * ===================================================================
 */
#ifndef __LINBOX_tri_solve_gf2_INL
#define __LINBOX_tri_solve_gf2_INL

#include "linbox/vector/vector-domain.h"
#include "linbox/field/gf2.h"

namespace LinBox
{
	template <class _Matrix, class Vector1, class Vector2> Vector1&
	upperTriangularSolveBinary (Vector1& x,
				    const _Matrix  &U,
				    const Vector2& b)
	{
		//         linbox_check( x.size() == U.coldim() );
		//         linbox_check( b.size() == U.rowdim() );
		typedef _Matrix Matrix;
		typedef GF2 Field;
		const GF2 F2;

		commentator().start ("Sparse Elimination Upper Triangular Solve over GF(2)", "utrsmGF2");

		typename Vector2::const_iterator vec=b.begin();
		typename Vector1::iterator res=x.begin();
		typename Matrix::const_iterator row=U.begin();

		// Find last constrained values of x, U and b
		//         for( ; (res != x.end()) && (row != U.rowEnd()); ++res, ++row, ++vec) { }
		size_t last = x.size();
		if( b.size() < last ) last = b.size();
		res += (ptrdiff_t)last;
		row += (ptrdiff_t)last;
		vec += (ptrdiff_t) last;

		VectorCategories::DenseZeroOneVectorTag  DZOtag;
		VectorCategories::SparseZeroOneVectorTag SZOtag;

		bool consistant = true;
		for(typename Vector2::const_iterator bcheck=vec; bcheck != b.end(); ++bcheck) {
			if( ! F2.isZero(*bcheck) ) {
				consistant = false;
				break;
			}
		}
		if (consistant) {
			--vec; --res; --row;

			VectorDomain<Field> VD(F2);
			for( ; row != U.begin(); --row, --vec, --res) {
				F2.assign(*res, F2.zero);
				if (row->size()) {
					typename Field::Element tmp;
					VD.dotSpecialized(tmp, x, *row, DZOtag, SZOtag);
					F2.addin(tmp,*vec);
					F2.assign(*res,tmp);
				}
				else {
					// Consistency check
					if( ! F2.isZero(*vec) ) {
						consistant = false;
						break;
					}
				}
			}
			F2.assign(*res, F2.zero);
			if (row->size()) {
				typename Field::Element tmp;
				VD.dotSpecialized(tmp, x, *row, DZOtag, SZOtag);
				F2.addin(tmp,*vec);
				F2.assign(*res,tmp);
			}
			else {
				// Consistency check
				if( ! F2.isZero(*vec) ) consistant = false;
			}
		}
		//         if (! consistant) throw LinboxError ("upperTriangularSolveBinary returned INCONSISTENT");
		linbox_check( consistant );

		commentator().stop ("done", NULL, "utrsmGF2");
		return x;
	}

	template <class _Matrix, class Vector1, class Vector2> Vector1&
	lowerTriangularUnitarySolveBinary (Vector1& x,
					   const _Matrix  &L,
					   const Vector2& b)
	{
		linbox_check( b.size() == L.coldim() );
		typedef _Matrix Matrix;
		const GF2 F2;

		commentator().start ("Sparse Elimination Lower Triangular Unitary Solve over GF2", "ltrsmGF2");

		typename Vector2::const_iterator vec=b.begin();
		typename Vector1::iterator res=x.begin();
		typename Matrix::const_iterator row=L.begin();

		VectorCategories::DenseZeroOneVectorTag  DZOtag;
		VectorCategories::SparseZeroOneVectorTag SZOtag;
		VectorDomain<GF2> VD(F2);
		for( ; row != L.end(); ++row, ++vec, ++res) {
			F2.assign(*res, F2.zero);
			GF2::Element tmp;
			VD.dotSpecialized(tmp, *row, x, SZOtag, DZOtag);
			F2.negin(tmp);
			F2.addin(tmp,*vec);
			F2.assign(*res,tmp);
		}

		commentator().stop ("done", NULL, "ltrsmGF2");
		return x;
	}

}
#endif //__LINBOX_tri_solve_gf2_INL

// Local Variables:
// mode: C++
// tab-width: 4
// indent-tabs-mode: nil
// c-basic-offset: 4
// End:
// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s