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     {
  56       int k=-s;
  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());
  61       for (; k <= s; k++)
  62         f.elem (i, k + s) =0.0;
  63     }
  64
  65   band_ = f;
  66 }
  67
  68
  69
  70 /*
  71   any takers?
  72   */
  73 void
  74 Diagonal_storage::insert_row (int)
  75 {
  76   assert (false);
  77 }
  78
  79 void
  80 Diagonal_storage::delete_row (int)
  81 {
  82   assert (false);
  83 }
  84
  85 void
  86 Diagonal_storage::resize (int,int)
  87 {
  88 }
  89
  90 void
  91 Diagonal_storage::resize (int)
  92 {
  93 }
  94
  95 void
  96 Diagonal_storage::delete_column (int)
  97 {
  98   assert (false);
  99 }
 100
 101 Diagonal_storage::~Diagonal_storage()
 102 {
 103 }
 104
 105
 106 bool
 107 Diagonal_storage::band_elt_b (int i,int j) const
 108 {
 109   return abs (i-j) <= band_size_i();
 110 }
 111
 112 void
 113 Diagonal_storage::assert_valid (int i,int j) const
 114 {
 115   assert (band_elt_b (i,j));
 116   assert (i >=0 && j >=0 && i < dim() && j < dim ());
 117 }
 118
 119
 120 void
 121 Diagonal_storage::resize_dim (int d)
 122 {
 123   Full_storage f (d, 2*band_size_i()+1);
 124   for (int i=0; i < d && i < dim(); i++)
 125     {
 126       for (int k=0;  k < 2*band_size_i(); k++)
 127         f.elem (i,k) = elem (i,k);
 128     }
 129
 130   band_ = f;
 131 }
 132
 133
 134
 135 bool
 136 Diagonal_storage::mult_ok (int i,int) const
 137 {
 138   return i < dim();
 139 }
 140
 141 void
 142 Diagonal_storage::mult_next (int &i, int &j) const
 143 {
 144   j++;
 145   if (j < i - band_size_i())
 146     j = i- band_size_i();
 147   if (j > i + band_size_i() || j >= dim ())
 148     {
 149       i++;
 150       j = 0 >? i - band_size_i();
 151     }
 152 }
 153
 154 bool
 155 Diagonal_storage::trans_ok (int ,int j) const
 156 {
 157   return j < dim();
 158 }
 159
 160 void
 161 Diagonal_storage::trans_next (int &i, int& j) const
 162 {
 163   i++;
 164   if (i < j - band_size_i())
 165     i = j-band_size_i();
 166
 167   if (i >= dim() || i > j + band_size_i ())
 168     {
 169       j++;
 170       i = 0 >? j - band_size_i();
 171     }
 172 }
 173
 174 static Real nul_entry=0.0;
 175
 176 Real
 177 Diagonal_storage::elem (int i, int j) const
 178 {
 179   if (abs (i-j) > band_size_i())
 180     return 0;
 181   else
 182     return band_.elem (i, j - i +band_size_i());
 183 }
 184
 185 Real &
 186 Diagonal_storage::elem (int i, int j)
 187 {
 188   /*
 189     if this fails, the previous call fucked up
 190     */
 191   assert (!nul_entry);
 192
 193   if (abs (i-j) > band_size_i())
 194     return nul_entry;
 195   else
 196     return band_.elem (i, j - i + band_size_i());
 197 }
 198
 199 /*
 200   Hairy routine to try to save some fp-ops
 201  */
 202
 203 bool
 204 Diagonal_storage::try_right_multiply (Matrix_storage*dest,
 205                                       const Matrix_storage*right) const
 206 {
 207   if (right->name() != Diagonal_storage::static_name ())
 208     return false;
 209
 210   const Diagonal_storage*  right_diag = (Diagonal_storage const*)right;
 211   int band2 = right_diag->band_size_i();
 212   int n = dim();
 213   /*
 214     should check if dest is a Diagonal_storage of sufficient size too.
 215     */
 216   for (int i=0;  i < n; i++)
 217     {
 218       for (int j = 0; j < n; j++)
 219         {
 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;
 223           Real sum =0.0;
 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;
 227
 228         }
 229     }
 230   return true;
 231 }
 232
 233 IMPLEMENT_IS_TYPE_B1(Diagonal_storage, Matrix_storage);
 234
 235
 236 Diagonal_storage::Diagonal_storage (Matrix_storage*stor_l, int band_i)
 237 {
 238   set_band_size (band_i);
 239   resize_dim (stor_l->dim());
 240
 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);
 243 }
 244
 245 void
 246 Diagonal_storage::OK() const
 247 {
 248   band_.OK();
 249 }
 250
 251 IMPLEMENT_VIRTUAL_COPY_CONS(Diagonal_storage, Matrix_storage);