source: trunk/linbox/linbox/algorithms/triangular-solve.h @ 4085

Revision 4085, 4.9 KB checked in by bboyer, 3 years ago (diff)

clang

Line 
1/* -*- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
2// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
3/* ===================================================================
4 * Copyright(C) 2008 LinBox
5 * Triangular Solve
6 * See COPYING for license information.
7 * Time-stamp: <16 Jun 10 14:21:18 Jean-Guillaume.Dumas@imag.fr>
8 * ===================================================================
9 */
10#ifndef __LINBOX_triangular_solve_H
11#define __LINBOX_triangular_solve_H
12
13#include "linbox/vector/vector-domain.h"
14
15namespace LinBox
16{
17        template <class _Matrix, class Vector1, class Vector2> Vector1&
18        upperTriangularSolve (Vector1& x,
19                              const _Matrix  &U,
20                              const Vector2& b)
21        {
22                linbox_check( x.size() == U.coldim() );
23                linbox_check( b.size() == U.rowdim() );
24                typedef _Matrix Matrix;
25                typedef typename Matrix::Field Field;
26                const Field& F = U.field();
27
28                commentator.start ("Sparse Elimination Upper Triangular Solve", "utrsm");
29
30                typename Vector2::const_iterator vec=b.begin();
31                typename Vector1::iterator res=x.begin();
32                typename Matrix::ConstRowIterator row=U.rowBegin();
33
34                // Assume U has form (U1, X | 0, 0), where U1 is invertible.
35                // Discover the rank of U so as to use the bottom of x to seed the solution.
36                for(row = U.rowEnd()-1; row >= U.rowBegin() && row->size() == 0; --row);
37                row++; // now points to first zero row of U.
38                res += row - U.rowBegin();
39                vec += row - U.rowBegin();
40
41                bool consistant = true;
42                for(typename Vector2::const_iterator bcheck=vec; bcheck != b.end(); ++bcheck) {
43                        if( ! F.isZero(*bcheck) ) {
44                                consistant = false;
45                                break;
46                        }
47                }
48                if (consistant) {
49                        --vec; --res; --row;
50
51                        VectorDomain<Field> VD(F);
52                        for( ; row != U.rowBegin(); --row, --vec, --res) {
53                                F.init(*res, 0UL);
54                                if (row->size()) {
55                                        typename Field::Element tmp;
56                                        VD.dot(tmp, *row, x);
57                                        F.negin(tmp);
58                                        F.addin(tmp,*vec);
59                                        F.divin(tmp,row->front().second);
60                                        F.assign(*res,tmp);
61                                }
62                                else {
63                                        // Consistency check
64                                        if( ! F.isZero(*vec) ) {
65                                                consistant = false;
66                                                break;
67                                        }
68                                }
69                        }
70
71                        F.init(*res, 0UL);
72                        if (row->size()) {
73                                typename Field::Element tmp;
74                                VD.dot(tmp, *row, x);
75                                F.negin(tmp);
76                                F.addin(tmp,*vec);
77                                F.divin(tmp,row->front().second);
78                                F.assign(*res,tmp);
79                        }
80                        else {
81                                // Consistency check
82                                if( ! F.isZero(*vec) ) consistant = false;
83                        }
84                }
85                if (! consistant) throw LinboxError ("upperTriangularSolve returned INCONSISTENT");
86
87                commentator.stop ("done", NULL, "utrsm");
88                return x;
89        }
90
91        // Suppose x and b are vectors of pairs <index,value>
92        // first rank rows of U are upper triangular and full rank
93        template <class _Matrix, class Vector1, class Vector2> Vector1&
94        upperTriangularSparseSolve (Vector1& x,
95                                    unsigned long rank,
96                                    const _Matrix  &U,
97                                    const Vector2& b)
98        {
99                commentator.start ("SparseElim UpperTriang Sparse Solve", "uSPt");
100
101                x.resize(0);
102
103                if (b.size() != 0) {
104
105                        typedef _Matrix Matrix;
106                        typedef typename Matrix::Field Field;
107                        const Field& F = U.field();
108
109
110
111                        typename Vector2::const_iterator vec=b.begin();
112                        vec += b.size(); --vec;
113
114                        typename Matrix::ConstRowIterator row=U.rowBegin();
115                        row += rank; --row;
116
117                        VectorDomain<Field> VD(F);
118
119
120                        long i=(long)rank;
121                        for(--i; (vec >= b.begin()) && (i>=0); --i,--row) {
122                                if (row->size()) {
123                                        typename Field::Element tmp;
124                                        VD.dot(tmp, *row, x); // x is sparse also
125                                        F.negin(tmp);
126                                        if (static_cast<long>(vec->first) == i) {
127                                                F.addin(tmp,vec->second);
128                                                --vec;
129                                        }
130                                        if (! F.isZero(tmp)) {
131                                                F.divin(tmp,row->front().second);
132                                                x.insert(x.begin(), typename Vector1::value_type((unsigned)i, tmp));
133                                        }
134                                }
135                        }
136                        for(; i>=0; --i,--row) {
137                                if (row->size()) {
138                                        typename Field::Element tmp;
139                                        VD.dot(tmp, *row, x);
140                                        if (! F.isZero(tmp)) {
141                                                F.negin(tmp);
142                                                F.divin(tmp,row->front().second);
143                                                x.insert(x.begin(), typename Vector1::value_type((unsigned)i, tmp));
144                                        }
145                                }
146                        }
147                }
148                //         if (! consistant) throw LinboxError ("upperTriangularSparseSolve returned INCONSISTENT");
149
150                commentator.stop ("done", NULL, "uSPt");
151
152                return x;
153        }
154
155
156        template <class _Matrix, class Vector1, class Vector2> Vector1&
157        lowerTriangularUnitarySolve (Vector1& x,
158                                     const _Matrix  &L,
159                                     const Vector2& b)
160        {
161                linbox_check( b.size() == L.coldim() );
162                typedef _Matrix Matrix;
163                typedef typename Matrix::Field Field;
164                const Field& F = L.field();
165
166                commentator.start ("Sparse Elimination Lower Triangular Unitary Solve", "ltrsm");
167
168                typename Vector2::const_iterator vec=b.begin();
169                typename Vector1::iterator res=x.begin();
170                typename Matrix::ConstRowIterator row=L.rowBegin();
171
172                VectorDomain<Field> VD(F);
173                for( ; row != L.rowEnd(); ++row, ++vec, ++res) {
174                        F.init(*res, 0UL);
175                        typename Field::Element tmp;
176                        VD.dot(tmp, *row, x);
177                        F.negin(tmp);
178                        F.addin(tmp,*vec);
179                        F.assign(*res,tmp);
180                }
181
182                commentator.stop ("done", NULL, "ltrsm");
183                return x;
184        }
185}
186#endif //__LINBOX_triangular_solve_H
187
Note: See TracBrowser for help on using the repository browser.