2 diagonal-storage.cc -- implement Diagonal_storage
4 source file of the Flower Library
6 (c) 1997 Han-Wen Nienhuys <hanwen@stack.nl>
10 #include "diagonal-storage.hh"
19 #include "full-storage.icc"
22 Diagonal_storage::dim() const
27 Diagonal_storage::Diagonal_storage()
32 Diagonal_storage::rows() const
38 Diagonal_storage::cols() const
44 Diagonal_storage::band_size_i() const
46 return (band_.cols()-1)/2;
50 Diagonal_storage::set_band_size (int s)
53 Full_storage f (dim(), 2*s+1);
54 for (int i=0; i < dim(); i++)
57 for (; k < -band_size_i(); k++)
58 f.elem (i,k + s) = 0.0;
59 for (; k <= band_size_i()&& k<=s ; k++)
60 f.elem (i, k + s) = band_.elem (i,k+ band_size_i());
62 f.elem (i, k + s) =0.0;
74 Diagonal_storage::insert_row (int)
80 Diagonal_storage::delete_row (int)
86 Diagonal_storage::resize (int,int)
91 Diagonal_storage::resize (int)
96 Diagonal_storage::delete_column (int)
101 Diagonal_storage::~Diagonal_storage()
107 Diagonal_storage::band_elt_b (int i,int j) const
109 return abs (i-j) <= band_size_i();
113 Diagonal_storage::assert_valid (int i,int j) const
115 assert (band_elt_b (i,j));
116 assert (i >=0 && j >=0 && i < dim() && j < dim ());
121 Diagonal_storage::resize_dim (int d)
123 Full_storage f (d, 2*band_size_i()+1);
124 for (int i=0; i < d && i < dim(); i++)
126 for (int k=0; k < 2*band_size_i(); k++)
127 f.elem (i,k) = elem (i,k);
136 Diagonal_storage::mult_ok (int i,int) const
142 Diagonal_storage::mult_next (int &i, int &j) const
145 if (j < i - band_size_i())
146 j = i- band_size_i();
147 if (j > i + band_size_i() || j >= dim ())
150 j = 0 >? i - band_size_i();
155 Diagonal_storage::trans_ok (int ,int j) const
161 Diagonal_storage::trans_next (int &i, int& j) const
164 if (i < j - band_size_i())
167 if (i >= dim() || i > j + band_size_i ())
170 i = 0 >? j - band_size_i();
174 static Real nul_entry=0.0;
177 Diagonal_storage::elem (int i, int j) const
179 if (abs (i-j) > band_size_i())
182 return band_.elem (i, j - i +band_size_i());
186 Diagonal_storage::elem (int i, int j)
189 if this fails, the previous call fucked up
193 if (abs (i-j) > band_size_i())
196 return band_.elem (i, j - i + band_size_i());
200 Hairy routine to try to save some fp-ops
204 Diagonal_storage::try_right_multiply (Matrix_storage*dest,
205 const Matrix_storage*right) const
207 if (right->name() != Diagonal_storage::static_name ())
210 const Diagonal_storage* right_diag = (Diagonal_storage const*)right;
211 int band2 = right_diag->band_size_i();
214 should check if dest is a Diagonal_storage of sufficient size too.
216 for (int i=0; i < n; i++)
218 for (int j = 0; j < n; j++)
220 int startk = i - band_size_i() >? 0 >? j - band2;
221 int stopk = i + band_size_i() <? n-1 <? j + band2;
222 int relk = startk + band_size_i() -i;
224 for (int k = startk; k <= stopk; k++)
225 sum += band_.elem (i, relk++) * right_diag->elem (k, j);
226 dest->elem (i, j) = sum;
233 IMPLEMENT_IS_TYPE_B1(Diagonal_storage, Matrix_storage);
236 Diagonal_storage::Diagonal_storage (Matrix_storage*stor_l, int band_i)
238 set_band_size (band_i);
239 resize_dim (stor_l->dim());
241 for (int i=0,j=0; mult_ok (i,j); mult_next (i,j))
242 band_.elem (i, j + band_i -i) = stor_l->elem (i,j);
246 Diagonal_storage::OK() const
251 IMPLEMENT_VIRTUAL_COPY_CONS(Diagonal_storage, Matrix_storage);