1 /* =================================================================== 2 * Copyright(C) LinBox 2008 3 * Written by Jean-Guillaume Dumas 4 * Triangular Solve 5 * 6 * ========LICENCE======== 7 * This file is part of the library LinBox. 8 * 9 * LinBox is free software: you can redistribute it and/or modify 10 * it under the terms of the GNU Lesser General Public 11 * License as published by the Free Software Foundation; either 12 * version 2.1 of the License, or (at your option) any later version. 13 * 14 * This library is distributed in the hope that it will be useful, 15 * but WITHOUT ANY WARRANTY; without even the implied warranty of 16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 17 * Lesser General Public License for more details. 18 * 19 * You should have received a copy of the GNU Lesser General Public 20 * License along with this library; if not, write to the Free Software 21 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 22 * ========LICENCE======== 23 *. 24 * Time-stamp: <01 Oct 09 15:38:25 Jean-Guillaume.Dumas@imag.fr> 25 * =================================================================== 26 */ 27 #ifndef __LINBOX_tri_solve_gf2_INL 28 #define __LINBOX_tri_solve_gf2_INL 29 30 #include "linbox/vector/vector-domain.h" 31 #include "linbox/field/gf2.h" 32 33 namespace LinBox 34 { 35 template <class _Matrix, class Vector1, class Vector2> Vector1& upperTriangularSolveBinary(Vector1 & x,const _Matrix & U,const Vector2 & b)36 upperTriangularSolveBinary (Vector1& x, 37 const _Matrix &U, 38 const Vector2& b) 39 { 40 // linbox_check( x.size() == U.coldim() ); 41 // linbox_check( b.size() == U.rowdim() ); 42 typedef _Matrix Matrix; 43 typedef GF2 Field; 44 const GF2 F2; 45 46 commentator().start ("Sparse Elimination Upper Triangular Solve over GF(2)", "utrsmGF2"); 47 48 typename Vector2::const_iterator vec=b.begin(); 49 typename Vector1::iterator res=x.begin(); 50 typename Matrix::const_iterator row=U.begin(); 51 52 // Find last constrained values of x, U and b 53 // for( ; (res != x.end()) && (row != U.rowEnd()); ++res, ++row, ++vec) { } 54 size_t last = x.size(); 55 if( b.size() < last ) last = b.size(); 56 res += (ptrdiff_t)last; 57 row += (ptrdiff_t)last; 58 vec += (ptrdiff_t) last; 59 60 VectorCategories::DenseZeroOneVectorTag DZOtag; 61 VectorCategories::SparseZeroOneVectorTag SZOtag; 62 63 bool consistant = true; 64 for(typename Vector2::const_iterator bcheck=vec; bcheck != b.end(); ++bcheck) { 65 if( ! F2.isZero(*bcheck) ) { 66 consistant = false; 67 break; 68 } 69 } 70 if (consistant) { 71 --vec; --res; --row; 72 73 VectorDomain<Field> VD(F2); 74 for( ; row != U.begin(); --row, --vec, --res) { 75 F2.assign(*res, F2.zero); 76 if (row->size()) { 77 typename Field::Element tmp; 78 VD.dotSpecialized(tmp, x, *row, DZOtag, SZOtag); 79 F2.addin(tmp,*vec); 80 F2.assign(*res,tmp); 81 } 82 else { 83 // Consistency check 84 if( ! F2.isZero(*vec) ) { 85 consistant = false; 86 break; 87 } 88 } 89 } 90 F2.assign(*res, F2.zero); 91 if (row->size()) { 92 typename Field::Element tmp; 93 VD.dotSpecialized(tmp, x, *row, DZOtag, SZOtag); 94 F2.addin(tmp,*vec); 95 F2.assign(*res,tmp); 96 } 97 else { 98 // Consistency check 99 if( ! F2.isZero(*vec) ) consistant = false; 100 } 101 } 102 // if (! consistant) throw LinboxError ("upperTriangularSolveBinary returned INCONSISTENT"); 103 linbox_check( consistant ); 104 105 commentator().stop ("done", NULL, "utrsmGF2"); 106 return x; 107 } 108 109 template <class _Matrix, class Vector1, class Vector2> Vector1& lowerTriangularUnitarySolveBinary(Vector1 & x,const _Matrix & L,const Vector2 & b)110 lowerTriangularUnitarySolveBinary (Vector1& x, 111 const _Matrix &L, 112 const Vector2& b) 113 { 114 linbox_check( b.size() == L.coldim() ); 115 typedef _Matrix Matrix; 116 const GF2 F2; 117 118 commentator().start ("Sparse Elimination Lower Triangular Unitary Solve over GF2", "ltrsmGF2"); 119 120 typename Vector2::const_iterator vec=b.begin(); 121 typename Vector1::iterator res=x.begin(); 122 typename Matrix::const_iterator row=L.begin(); 123 124 VectorCategories::DenseZeroOneVectorTag DZOtag; 125 VectorCategories::SparseZeroOneVectorTag SZOtag; 126 VectorDomain<GF2> VD(F2); 127 for( ; row != L.end(); ++row, ++vec, ++res) { 128 F2.assign(*res, F2.zero); 129 GF2::Element tmp; 130 VD.dotSpecialized(tmp, *row, x, SZOtag, DZOtag); 131 F2.negin(tmp); 132 F2.addin(tmp,*vec); 133 F2.assign(*res,tmp); 134 } 135 136 commentator().stop ("done", NULL, "ltrsmGF2"); 137 return x; 138 } 139 140 } 141 #endif //__LINBOX_tri_solve_gf2_INL 142 143 // Local Variables: 144 // mode: C++ 145 // tab-width: 4 146 // indent-tabs-mode: nil 147 // c-basic-offset: 4 148 // End: 149 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 150