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