flower/diagonal-storage.cc

   1 /*
   2   diagonal-storage.cc -- implement Diagonal_storage
   3
   4   source file of the Flower Library
   5
   6   (c) 1997 Han-Wen Nienhuys <hanwen@stack.nl>
   7 */
   8
   9
  10 #include "diagonal-storage.hh"
  11
  12
  13 #ifdef INLINE
  14 #undef INLINE
  15 #endif
  16
  17 #define INLINE inline
  18
  19 #include "full-storage.icc"
  20
  21 int
  22 Diagonal_storage::dim()const
  23 {
  24     return band_.rows();
  25 }
  26
  27 Diagonal_storage::Diagonal_storage()
  28 {
  29 }
  30
  31 int
  32 Diagonal_storage::rows() const
  33 {
  34     return band_.rows();
  35 }
  36
  37 int
  38 Diagonal_storage::cols() const
  39 {
  40     return band_.rows();
  41 }
  42
  43 int
  44 Diagonal_storage::band_size_i()const
  45 {
  46     return (band_.cols()-1)/2;
  47 }
  48
  49 void
  50 Diagonal_storage::set_band_size(int s)
  51 {
  52     assert( s>=0);
  53     Full_storage f(dim(), 2*s+1);
  54     for (int i=0; i < dim(); i++) {
  55         int k=-s;
  56         for ( ;  k < -band_size_i(); k++)
  57             f.elem(i,k + s) = 0.0;
  58         for ( ; k <=  band_size_i()&& k<=s ; k++)
  59             f.elem(i, k + s) = band_.elem(i,k+ band_size_i());
  60         for( ; k <= s; k++)
  61             f.elem(i, k + s) =0.0;
  62     }
  63
  64     band_ = f;
  65 }
  66
  67
  68
  69 /*
  70   any takers?
  71   */
  72 void
  73 Diagonal_storage::insert_row(int )
  74 {
  75     assert(false);
  76 }
  77
  78 void
  79 Diagonal_storage::delete_row(int )
  80 {
  81     assert(false);
  82 }
  83
  84 void
  85 Diagonal_storage::resize(int,int)
  86 {
  87 }
  88
  89 void
  90 Diagonal_storage::resize(int)
  91 {
  92 }
  93
  94 void
  95 Diagonal_storage::delete_column(int )
  96 {
  97     assert(false);
  98 }
  99
 100 Diagonal_storage::~Diagonal_storage()
 101 {
 102 }
 103
 104
 105 bool
 106 Diagonal_storage::band_elt_b(int i,int j)const
 107 {
 108     return abs(i-j) <= band_size_i();
 109 }
 110
 111 void
 112 Diagonal_storage::assert_valid(int i,int j )const
 113 {
 114     assert( band_elt_b(i,j) );
 115     assert( i >=0 && j >=0 && i < dim() && j < dim());
 116 }
 117
 118
 119 void
 120 Diagonal_storage::resize_dim(int d)
 121 {
 122     Full_storage f(d, 2*band_size_i()+1);
 123     for (int i=0; i < d && i < dim(); i++) {
 124         for ( int k=0;  k < 2*band_size_i(); k++)
 125             f.elem(i,k) = elem(i,k);
 126     }
 127
 128     band_ = f;
 129 }
 130
 131
 132
 133 bool
 134 Diagonal_storage::mult_ok(int i,int )const
 135 {
 136     return i < dim();
 137 }
 138
 139 void
 140 Diagonal_storage::mult_next(int &i, int &j)const
 141 {
 142     j++;
 143     if ( j < i - band_size_i() )
 144         j = i- band_size_i();
 145     if ( j > i + band_size_i() || j >= dim() ) {
 146         i++;
 147         j = 0 >? i - band_size_i();
 148     }
 149 }
 150
 151 bool
 152 Diagonal_storage::trans_ok(int ,int j)const
 153 {
 154     return j < dim();
 155 }
 156
 157 void
 158 Diagonal_storage::trans_next(int &i, int& j)const
 159 {
 160     i++;
 161     if ( i < j - band_size_i())
 162         i = j-band_size_i();
 163
 164     if ( i >= dim() || i > j + band_size_i() ) {
 165         j++;
 166         i = 0 >? j - band_size_i();
 167     }
 168 }
 169
 170 static Real nul_entry=0.0;
 171
 172 Real
 173 Diagonal_storage::elem(int i, int j)const
 174 {
 175     if (abs ( i-j ) > band_size_i())
 176         return 0;
 177     else
 178         return band_.elem(i, j - i +band_size_i());
 179 }
 180
 181 Real &
 182 Diagonal_storage::elem(int i, int j)
 183 {
 184     /*
 185       if this fails, the previous call fucked up
 186       */
 187     assert(!nul_entry);
 188
 189     if (abs ( i-j ) > band_size_i())
 190         return nul_entry;
 191     else
 192         return band_.elem(i, j - i + band_size_i());
 193 }
 194
 195 /*
 196   Hairy routine to try to save some fp-ops
 197  */
 198
 199 bool
 200 Diagonal_storage::try_right_multiply(Matrix_storage*dest,
 201                                      const Matrix_storage*right)const
 202 {
 203     if ( right->name() != Diagonal_storage::static_name() )
 204         return false;
 205
 206     const Diagonal_storage*  right_diag = (Diagonal_storage const*)right;
 207     int band2 = right_diag->band_size_i();
 208     int n = dim();
 209     /*
 210       should check if dest is a Diagonal_storage of sufficient size too.
 211      */
 212     for (int i=0;  i < n; i++) {
 213         for (int j = 0; j < n; j++) {
 214             int startk = i - band_size_i() >? 0 >? j - band2;
 215             int stopk = i + band_size_i() <? n-1 <? j  + band2;
 216             int relk =  startk + band_size_i() -i;
 217             Real sum =0.0;
 218             for ( int k = startk; k <= stopk; k++)
 219                 sum += band_.elem(i, relk++) * right_diag->elem(k, j);
 220             dest->elem(i, j) = sum;
 221
 222         }
 223     }
 224     return true;
 225 }
 226
 227 IMPLEMENT_IS_TYPE_B1(Diagonal_storage, Matrix_storage);
 228
 229
 230 Diagonal_storage::Diagonal_storage(Matrix_storage*stor_l, int band_i)
 231 {
 232     set_band_size(band_i);
 233     resize_dim(stor_l->dim());
 234
 235     for ( int i=0,j=0; mult_ok(i,j); mult_next(i,j))
 236         band_.elem(i, j + band_i -i ) = stor_l->elem(i,j);
 237 }
 238
 239 void
 240 Diagonal_storage::OK() const
 241 {
 242     band_.OK();
 243 }
 244
 245 IMPLEMENT_VIRTUAL_COPY_CONS(Diagonal_storage, Matrix_storage);