root / trunk / linbox / linbox / algorithms / gauss.h

Revision 3286, 13.0 kB (checked in by dumas, 3 weeks ago)

default solve now is NOT random, put 0 instead
To get random solve set randomsol parameter to true
Room for improvement in upperTriangularSolve when 0 ...

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1/* -*- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
2
3/* linbox/algorithms/gauss.h
4 * Copyright (C) 1999 Jean-Guillaume Dumas
5 *
6 * Written by Jean-Guillaume Dumas <Jean-Guillaume.Dumas@imag.fr>
7 *
8 * -----------------------------------------------------------
9 * 2003-02-02  Bradford Hovinen  <bghovinen@math.uwaterloo.ca>
10 *
11 * Ported to new matrix archetype; update interface to meet current
12 * standards. Rename gauss_foo as foo and gauss_Uin as upperin
13 *
14 * Move function definitions to gauss.inl
15 * -----------------------------------------------------------
16 *
17 * See COPYING for license information.
18 */
19
20// ========================================================================= //
21// (C) The Linbox Group 1999
22// Calcul de rang par la méthode de Gauss pivot par ligne, sur matrice creuse
23// Time-stamp: <03 Nov 00 19:19:06 Jean-Guillaume.Dumas@imag.fr>
24// ========================================================================= //
25
26#ifndef __GAUSS_H
27#define __GAUSS_H
28
29#include "linbox/util/debug.h"
30#include "linbox/util/commentator.h"
31#include "linbox/field/archetype.h"
32#include "linbox/field/gf2.h"
33#include "linbox/vector/vector-domain.h"
34#include "linbox/matrix/sparse.h"
35#include "linbox/matrix/matrix-domain.h"
36#include "linbox/matrix/archetype.h"
37#include "linbox/solutions/methods.h"
38
39namespace LinBox 
40{
41
42/** \brief Repository of functions for rank by elimination on sparse matrices.
43
44    Several versions allow for adjustment of the pivoting strategy
45    and for choosing in-place elimination or for not modifying the input matrix.
46    Also an LU interface is offered.
47*/
48    template <class _Field>
49    class GaussDomain {
50    public:
51        typedef _Field Field;
52        typedef typename Field::Element Element;
53
54    private:
55        const Field         &_F;
56
57    public:
58
59            /** \brief The field parameter is the domain
60             * over which to perform computations
61             */
62        GaussDomain (const Field &F) : _F (F) {}
63
64            //Copy constructor
65            ///
66        GaussDomain (const GaussDomain &M) : _F (M._F) {}
67
68            /** accessor for the field of computation
69             */
70        const Field &field () const { return _F; }
71
72            /** @name rank
73                Callers of the different rank routines\\
74                -/ The "in" suffix indicates in place computation\\
75                -/ Without Ni, Nj, the Matrix parameter must be a vector of sparse
76                row vectors, NOT storing any zero.\\
77                -/ Calls {@link rankinLinearPivoting rankinLinearPivoting} (by default) or {@link rankinNoReordering rankinNoReordering}
78            */
79            //@{
80            ///     
81        template <class Matrix> unsigned long& rankin(unsigned long &rank,
82                                                      Matrix        &A,
83                                                      SparseEliminationTraits::PivotStrategy   reord = SparseEliminationTraits::PIVOT_LINEAR) const;
84            ///
85        template <class Matrix> unsigned long& rankin(unsigned long &rank,
86                                                      Matrix        &A,
87                                                      unsigned long  Ni,
88                                                      unsigned long  Nj,
89                                                      SparseEliminationTraits::PivotStrategy   reord = SparseEliminationTraits::PIVOT_LINEAR) const;
90            ///       
91        template <class Matrix> unsigned long& rank(unsigned long &rank,
92                                                    const Matrix        &A,
93                                                    SparseEliminationTraits::PivotStrategy   reord = SparseEliminationTraits::PIVOT_LINEAR) const;
94            ///       
95        template <class Matrix> unsigned long& rank(unsigned long &rank,
96                                                    const Matrix        &A,
97                                                    unsigned long  Ni,
98                                                    unsigned long  Nj,
99                                                    SparseEliminationTraits::PivotStrategy   reord = SparseEliminationTraits::PIVOT_LINEAR) const;
100            //@}
101
102            /** @name det
103                Callers of the different determinant routines\\
104                -/ The "in" suffix indicates in place computation\\
105                -/ Without Ni, Nj, the Matrix parameter must be a vector of sparse
106                row vectors, NOT storing any zero.\\
107                -/ Calls {@link LinearPivoting } (by default) or {@link NoReordering}
108            */
109            //@{
110            ///     
111        template <class Matrix> Element& detin(Element &determinant,
112                                               Matrix        &A,
113                                               SparseEliminationTraits::PivotStrategy   reord = SparseEliminationTraits::PIVOT_LINEAR) const;
114            ///
115        template <class Matrix> Element& detin(Element &determinant,
116                                               Matrix        &A,
117                                               unsigned long  Ni,
118                                               unsigned long  Nj,
119                                               SparseEliminationTraits::PivotStrategy   reord = SparseEliminationTraits::PIVOT_LINEAR) const;
120            ///       
121        template <class Matrix> Element& det(Element &determinant,
122                                             const Matrix        &A,
123                                             SparseEliminationTraits::PivotStrategy   reord = SparseEliminationTraits::PIVOT_LINEAR) const;
124            ///       
125        template <class Matrix> Element& det(Element &determinant,
126                                             const Matrix        &A,
127                                             unsigned long  Ni,
128                                             unsigned long  Nj,
129                                             SparseEliminationTraits::PivotStrategy   reord = SparseEliminationTraits::PIVOT_LINEAR) const;
130            //@}
131
132
133            /** \brief Sparse in place Gaussian elimination with reordering to reduce fill-in.
134                pivots are chosen in sparsest column of sparsest row.
135                This runs in linear overhead.
136                It is similar in spirit but different from Markovitz' approach. 
137
138                <pre>
139                Using : SparseFindPivot(..., density) for sparsest column, and
140                eliminate (..., density)
141                </pre>
142
143                The Matrix parameter must meet the LinBox sparse matrix interface.
144                [check details].
145                The computedet indicates whether the algorithm must compute the determionant as it goes
146
147                @ref [Jean-Guillaume Dumas and Gilles Villard,
148                Computing the rank of sparse matrices over finite fields.
149                In Ganzha et~al. CASC'2002, pages 47--62.]
150            */
151        template <class Matrix, class Perm>
152        unsigned long& QLUPin(unsigned long &rank,
153                              Element& determinant,
154                              Perm          &Q,
155                              Matrix        &L,
156                              Matrix        &U,
157                              Perm          &P,
158                              unsigned long Ni, 
159                              unsigned long Nj) const;
160
161        template <class Matrix, class Perm, class Vector1, class Vector2> 
162        Vector1& solve(Vector1& x, unsigned long rank, const Perm& Q, const Matrix& L, const Matrix& U, const Perm& P, const Vector2& b, bool randomsol=false)  const;
163       
164        template <class Matrix, class Perm, class Vector1, class Vector2> 
165        Vector1& solve(Vector1& x, Vector1& w, unsigned long rank, const Perm& Q, const Matrix& L, const Matrix& U, const Perm& P, const Vector2& b)  const;
166       
167
168        template <class Matrix, class Vector1, class Vector2>
169        Vector1& solvein(Vector1        &x,
170                         Matrix         &A,
171                         const Vector2  &b, bool randomsol=false)  const;
172
173
174        template <class Matrix>
175        unsigned long& InPlaceLinearPivoting(unsigned long &rank,
176                                              Element& determinant,
177                                              Matrix        &A,
178                                              unsigned long Ni, 
179                                              unsigned long Nj) const;
180
181
182            /** \brief Sparse Gaussian elimination without reordering.
183
184                Gaussian elimination is done on a copy of the matrix.
185                Using : SparseFindPivot
186                eliminate
187
188                Requirements : SLA is an array of sparse rows
189                WARNING : NOT IN PLACE, THERE IS A COPY.
190                Without reordering (Pivot is first non-zero in row)
191
192            */
193        template <class Matrix>
194        unsigned long& NoReordering (unsigned long &rank, Element& determinant, Matrix &LigneA, unsigned long Ni, unsigned long Nj) const;
195
196            /** \brief Dense in place LU factorization without reordering
197
198                Using : FindPivot and LU
199            */
200        template <class Matrix>
201        unsigned long &LUin (unsigned long &rank, Matrix &A) const;
202
203
204            /** \brief Dense in place Gaussian elimination without reordering
205
206                Using : FindPivot and LU
207            */
208        template <class Matrix>
209        unsigned long &upperin (unsigned long &rank, Matrix &A) const;
210       
211
212
213    protected:
214   
215            //-----------------------------------------
216            // Sparse elimination using a pivot row :
217            // lc <-- lc - lc[k]/lp[0] * lp
218            // D is the number of elements per column
219            //   it is updated and used for reordering
220            // Vector is a vector of Pair (lin_pair.h)
221            //-----------------------------------------
222        template <class Vector, class D>
223        void eliminate (Element             & headpivot,
224                        Vector              &lignecourante,
225                        const Vector        &lignepivot,
226                        const unsigned long indcol,
227                        const long indpermut,
228                        const unsigned long npiv,
229                        D                   &columns) const;
230
231        template <class Vector, class D>
232        void eliminate (Vector              &lignecourante,
233                        const Vector        &lignepivot,
234                        const unsigned long &indcol,
235                        const long &indpermut,
236                        D                   &columns) const;
237
238        template <class Vector>
239        void permute (Vector              &lignecourante,
240                      const unsigned long &indcol,
241                      const long &indpermut) const;
242            //-----------------------------------------
243            // Sparse elimination using a pivot row :
244            // lc <-- lc - lc[k]/lp[0] * lp
245            // No density update
246            // Vector is a vector of Pair (lin_pair.h)
247            //-----------------------------------------
248        template <class Vector>
249        void eliminate (Vector              &lignecourante,
250                        const Vector        &lignepivot,
251                        const unsigned long &indcol,
252                        const long &indpermut) const;
253
254            //-----------------------------------------
255            // Dense elimination using a pivot row :
256            // lc <-- lc - lc[k]/lp[0] * lp
257            // Computing only for k to n (and not 0 to n in LU)
258            //-----------------------------------------
259        template<class Vector>
260        void Upper (Vector       &lignecur,
261                    const Vector &lignepivot,
262                    unsigned long indcol,
263                    long indpermut) const;
264
265            //-----------------------------------------
266            // Dense elimination using a pivot row :
267            // lc <-- lc - lc[k]/lp[0] * lp
268            //-----------------------------------------
269        template <class Vector>
270        void LU (Vector       &lignecur,
271                 const Vector &lignepivot,
272                 unsigned long indcol,
273                 long indpermut) const;
274
275            //------------------------------------------
276            // Looking for a non-zero pivot in a row
277            // Using the column density for reordering
278            // Pivot is chosen as to :
279            // 1. Row density is minimum
280            // 2. Column density is minimum for this row
281            //------------------------------------------
282        template <class Vector, class D>
283        void SparseFindPivot (Vector &lignepivot, unsigned long &indcol, long &indpermut, D &columns, Element& determinant) const;
284       
285            //------------------------------------------
286            // Looking for a non-zero pivot in a row
287            // No reordering
288            //------------------------------------------
289        template <class Vector>
290        void SparseFindPivot (Vector &lignepivot, unsigned long &indcol, long &indpermut, Element& determinant) const;
291       
292            //------------------------------------------
293            // Looking for a non-zero pivot in a row 
294            // Dense search
295            //------------------------------------------
296        template <class Vector>
297        void FindPivot (Vector &lignepivot, unsigned long &k, long &indpermut) const;
298
299    };
300
301} // namespace LinBox
302
303#include "linbox/algorithms/gauss-pivot.inl"
304#include "linbox/algorithms/gauss-elim.inl"
305#include "linbox/algorithms/gauss-solve.inl"
306#include "linbox/algorithms/gauss.inl"
307#include "linbox/algorithms/gauss-rank.inl"
308#include "linbox/algorithms/gauss-det.inl"
309
310#endif // __GAUSS_H
Note: See TracBrowser for help on using the browser.