|
Blender
V2.59
|
00001 00005 /* 00006 00007 * 00008 * Template Numerical Toolkit (TNT): Linear Algebra Module 00009 * 00010 * Mathematical and Computational Sciences Division 00011 * National Institute of Technology, 00012 * Gaithersburg, MD USA 00013 * 00014 * 00015 * This software was developed at the National Institute of Standards and 00016 * Technology (NIST) by employees of the Federal Government in the course 00017 * of their official duties. Pursuant to title 17 Section 105 of the 00018 * United States Code, this software is not subject to copyright protection 00019 * and is in the public domain. The Template Numerical Toolkit (TNT) is 00020 * an experimental system. NIST assumes no responsibility whatsoever for 00021 * its use by other parties, and makes no guarantees, expressed or implied, 00022 * about its quality, reliability, or any other characteristic. 00023 * 00024 * BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE 00025 * see http://math.nist.gov/tnt for latest updates. 00026 * 00027 */ 00028 00029 00030 // 2D Regions for arrays and matrices 00031 00032 #ifndef REGION2D_H 00033 #define REGION2D_H 00034 00035 #include "index.h" 00036 #include <iostream> 00037 #include <cassert> 00038 00039 namespace TNT 00040 { 00041 00042 template <class Array2D> 00043 class const_Region2D; 00044 00045 00046 template <class Array2D> 00047 class Region2D 00048 { 00049 protected: 00050 00051 Array2D & A_; 00052 Subscript offset_[2]; // 0-offset internally 00053 Subscript dim_[2]; 00054 00055 public: 00056 typedef typename Array2D::value_type T; 00057 typedef Subscript size_type; 00058 typedef T value_type; 00059 typedef T element_type; 00060 typedef T* pointer; 00061 typedef T* iterator; 00062 typedef T& reference; 00063 typedef const T* const_iterator; 00064 typedef const T& const_reference; 00065 00066 Array2D & array() { return A_; } 00067 const Array2D & array() const { return A_; } 00068 Subscript lbound() const { return A_.lbound(); } 00069 Subscript num_rows() const { return dim_[0]; } 00070 Subscript num_cols() const { return dim_[1]; } 00071 Subscript offset(Subscript i) const // 1-offset 00072 { 00073 #ifdef TNT_BOUNDS_CHECK 00074 assert( A_.lbound() <= i); 00075 assert( i<= dim_[0] + A_.lbound()-1); 00076 #endif 00077 return offset_[i-A_.lbound()]; 00078 } 00079 00080 Subscript dim(Subscript i) const 00081 { 00082 #ifdef TNT_BOUNDS_CHECK 00083 assert( A_.lbound() <= i); 00084 assert( i<= dim_[0] + A_.lbound()-1); 00085 #endif 00086 return dim_[i-A_.lbound()]; 00087 } 00088 00089 00090 00091 Region2D(Array2D &A, Subscript i1, Subscript i2, Subscript j1, 00092 Subscript j2) : A_(A) 00093 { 00094 #ifdef TNT_BOUNDS_CHECK 00095 assert( i1 <= i2 ); 00096 assert( j1 <= j2); 00097 assert( A.lbound() <= i1); 00098 assert( i2<= A.dim(A.lbound()) + A.lbound()-1); 00099 assert( A.lbound() <= j1); 00100 assert( j2<= A.dim(A.lbound()+1) + A.lbound()-1 ); 00101 #endif 00102 00103 00104 offset_[0] = i1-A.lbound(); 00105 offset_[1] = j1-A.lbound(); 00106 dim_[0] = i2-i1+1; 00107 dim_[1] = j2-j1+1; 00108 } 00109 00110 Region2D(Array2D &A, const Index1D &I, const Index1D &J) : A_(A) 00111 { 00112 #ifdef TNT_BOUNDS_CHECK 00113 assert( I.lbound() <= I.ubound() ); 00114 assert( J.lbound() <= J.ubound() ); 00115 assert( A.lbound() <= I.lbound()); 00116 assert( I.ubound()<= A.dim(A.lbound()) + A.lbound()-1); 00117 assert( A.lbound() <= J.lbound()); 00118 assert( J.ubound() <= A.dim(A.lbound()+1) + A.lbound()-1 ); 00119 #endif 00120 00121 offset_[0] = I.lbound()-A.lbound(); 00122 offset_[1] = J.lbound()-A.lbound(); 00123 dim_[0] = I.ubound() - I.lbound() + 1; 00124 dim_[1] = J.ubound() - J.lbound() + 1; 00125 } 00126 00127 Region2D(Region2D<Array2D> &A, Subscript i1, Subscript i2, 00128 Subscript j1, Subscript j2) : A_(A.A_) 00129 { 00130 #ifdef TNT_BOUNDS_CHECK 00131 assert( i1 <= i2 ); 00132 assert( j1 <= j2); 00133 assert( A.lbound() <= i1); 00134 assert( i2<= A.dim(A.lbound()) + A.lbound()-1); 00135 assert( A.lbound() <= j1); 00136 assert( j2<= A.dim(A.lbound()+1) + A.lbound()-1 ); 00137 #endif 00138 offset_[0] = (i1 - A.lbound()) + A.offset_[0]; 00139 offset_[1] = (j1 - A.lbound()) + A.offset_[1]; 00140 dim_[0] = i2-i1 + 1; 00141 dim_[1] = j2-j1+1; 00142 } 00143 00144 Region2D<Array2D> operator()(Subscript i1, Subscript i2, 00145 Subscript j1, Subscript j2) 00146 { 00147 #ifdef TNT_BOUNDS_CHECK 00148 assert( i1 <= i2 ); 00149 assert( j1 <= j2); 00150 assert( A_.lbound() <= i1); 00151 assert( i2<= dim_[0] + A_.lbound()-1); 00152 assert( A_.lbound() <= j1); 00153 assert( j2<= dim_[1] + A_.lbound()-1 ); 00154 #endif 00155 return Region2D<Array2D>(A_, 00156 i1+offset_[0], offset_[0] + i2, 00157 j1+offset_[1], offset_[1] + j2); 00158 } 00159 00160 00161 Region2D<Array2D> operator()(const Index1D &I, 00162 const Index1D &J) 00163 { 00164 #ifdef TNT_BOUNDS_CHECK 00165 assert( I.lbound() <= I.ubound() ); 00166 assert( J.lbound() <= J.ubound() ); 00167 assert( A_.lbound() <= I.lbound()); 00168 assert( I.ubound()<= dim_[0] + A_.lbound()-1); 00169 assert( A_.lbound() <= J.lbound()); 00170 assert( J.ubound() <= dim_[1] + A_.lbound()-1 ); 00171 #endif 00172 00173 return Region2D<Array2D>(A_, I.lbound()+offset_[0], 00174 offset_[0] + I.ubound(), offset_[1]+J.lbound(), 00175 offset_[1] + J.ubound()); 00176 } 00177 00178 inline T & operator()(Subscript i, Subscript j) 00179 { 00180 #ifdef TNT_BOUNDS_CHECK 00181 assert( A_.lbound() <= i); 00182 assert( i<= dim_[0] + A_.lbound()-1); 00183 assert( A_.lbound() <= j); 00184 assert( j<= dim_[1] + A_.lbound()-1 ); 00185 #endif 00186 return A_(i+offset_[0], j+offset_[1]); 00187 } 00188 00189 inline const T & operator() (Subscript i, Subscript j) const 00190 { 00191 #ifdef TNT_BOUNDS_CHECK 00192 assert( A_.lbound() <= i); 00193 assert( i<= dim_[0] + A_.lbound()-1); 00194 assert( A_.lbound() <= j); 00195 assert( j<= dim_[1] + A_.lbound()-1 ); 00196 #endif 00197 return A_(i+offset_[0], j+offset_[1]); 00198 } 00199 00200 00201 Region2D<Array2D> & operator=(const Region2D<Array2D> &R) 00202 { 00203 Subscript M = num_rows(); 00204 Subscript N = num_cols(); 00205 00206 // make sure both sides conform 00207 assert(M == R.num_rows()); 00208 assert(N == R.num_cols()); 00209 00210 00211 Subscript start = R.lbound(); 00212 Subscript Mend = start + M - 1; 00213 Subscript Nend = start + N - 1; 00214 for (Subscript i=start; i<=Mend; i++) 00215 for (Subscript j=start; j<=Nend; j++) 00216 (*this)(i,j) = R(i,j); 00217 00218 return *this; 00219 } 00220 00221 Region2D<Array2D> & operator=(const const_Region2D<Array2D> &R) 00222 { 00223 Subscript M = num_rows(); 00224 Subscript N = num_cols(); 00225 00226 // make sure both sides conform 00227 assert(M == R.num_rows()); 00228 assert(N == R.num_cols()); 00229 00230 00231 Subscript start = R.lbound(); 00232 Subscript Mend = start + M - 1; 00233 Subscript Nend = start + N - 1; 00234 for (Subscript i=start; i<=Mend; i++) 00235 for (Subscript j=start; j<=Nend; j++) 00236 (*this)(i,j) = R(i,j); 00237 00238 return *this; 00239 } 00240 00241 Region2D<Array2D> & operator=(const Array2D &R) 00242 { 00243 Subscript M = num_rows(); 00244 Subscript N = num_cols(); 00245 00246 // make sure both sides conform 00247 assert(M == R.num_rows()); 00248 assert(N == R.num_cols()); 00249 00250 00251 Subscript start = R.lbound(); 00252 Subscript Mend = start + M - 1; 00253 Subscript Nend = start + N - 1; 00254 for (Subscript i=start; i<=Mend; i++) 00255 for (Subscript j=start; j<=Nend; j++) 00256 (*this)(i,j) = R(i,j); 00257 00258 return *this; 00259 } 00260 00261 Region2D<Array2D> & operator=(const T &scalar) 00262 { 00263 Subscript start = lbound(); 00264 Subscript Mend = lbound() + num_rows() - 1; 00265 Subscript Nend = lbound() + num_cols() - 1; 00266 00267 00268 for (Subscript i=start; i<=Mend; i++) 00269 for (Subscript j=start; j<=Nend; j++) 00270 (*this)(i,j) = scalar; 00271 00272 return *this; 00273 } 00274 00275 }; 00276 00277 //************************ 00278 00279 template <class Array2D> 00280 class const_Region2D 00281 { 00282 protected: 00283 00284 const Array2D & A_; 00285 Subscript offset_[2]; // 0-offset internally 00286 Subscript dim_[2]; 00287 00288 public: 00289 typedef typename Array2D::value_type T; 00290 typedef T value_type; 00291 typedef T element_type; 00292 typedef const T* const_iterator; 00293 typedef const T& const_reference; 00294 00295 const Array2D & array() const { return A_; } 00296 Subscript lbound() const { return A_.lbound(); } 00297 Subscript num_rows() const { return dim_[0]; } 00298 Subscript num_cols() const { return dim_[1]; } 00299 Subscript offset(Subscript i) const // 1-offset 00300 { 00301 #ifdef TNT_BOUNDS_CHECK 00302 assert( TNT_BASE_OFFSET <= i); 00303 assert( i<= dim_[0] + TNT_BASE_OFFSET-1); 00304 #endif 00305 return offset_[i-TNT_BASE_OFFSET]; 00306 } 00307 00308 Subscript dim(Subscript i) const 00309 { 00310 #ifdef TNT_BOUNDS_CHECK 00311 assert( TNT_BASE_OFFSET <= i); 00312 assert( i<= dim_[0] + TNT_BASE_OFFSET-1); 00313 #endif 00314 return dim_[i-TNT_BASE_OFFSET]; 00315 } 00316 00317 00318 const_Region2D(const Array2D &A, Subscript i1, Subscript i2, 00319 Subscript j1, Subscript j2) : A_(A) 00320 { 00321 #ifdef TNT_BOUNDS_CHECK 00322 assert( i1 <= i2 ); 00323 assert( j1 <= j2); 00324 assert( TNT_BASE_OFFSET <= i1); 00325 assert( i2<= A.dim(TNT_BASE_OFFSET) + TNT_BASE_OFFSET-1); 00326 assert( TNT_BASE_OFFSET <= j1); 00327 assert( j2<= A.dim(TNT_BASE_OFFSET+1) + TNT_BASE_OFFSET-1 ); 00328 #endif 00329 00330 offset_[0] = i1-TNT_BASE_OFFSET; 00331 offset_[1] = j1-TNT_BASE_OFFSET; 00332 dim_[0] = i2-i1+1; 00333 dim_[1] = j2-j1+1; 00334 } 00335 00336 const_Region2D(const Array2D &A, const Index1D &I, const Index1D &J) 00337 : A_(A) 00338 { 00339 #ifdef TNT_BOUNDS_CHECK 00340 assert( I.lbound() <= I.ubound() ); 00341 assert( J.lbound() <= J.ubound() ); 00342 assert( TNT_BASE_OFFSET <= I.lbound()); 00343 assert( I.ubound()<= A.dim(TNT_BASE_OFFSET) + TNT_BASE_OFFSET-1); 00344 assert( TNT_BASE_OFFSET <= J.lbound()); 00345 assert( J.ubound() <= A.dim(TNT_BASE_OFFSET+1) + TNT_BASE_OFFSET-1 ); 00346 #endif 00347 00348 offset_[0] = I.lbound()-TNT_BASE_OFFSET; 00349 offset_[1] = J.lbound()-TNT_BASE_OFFSET; 00350 dim_[0] = I.ubound() - I.lbound() + 1; 00351 dim_[1] = J.ubound() - J.lbound() + 1; 00352 } 00353 00354 00355 const_Region2D(const_Region2D<Array2D> &A, Subscript i1, 00356 Subscript i2, 00357 Subscript j1, Subscript j2) : A_(A.A_) 00358 { 00359 #ifdef TNT_BOUNDS_CHECK 00360 assert( i1 <= i2 ); 00361 assert( j1 <= j2); 00362 assert( TNT_BASE_OFFSET <= i1); 00363 assert( i2<= A.dim(TNT_BASE_OFFSET) + TNT_BASE_OFFSET-1); 00364 assert( TNT_BASE_OFFSET <= j1); 00365 assert( j2<= A.dim(TNT_BASE_OFFSET+1) + TNT_BASE_OFFSET-1 ); 00366 #endif 00367 offset_[0] = (i1 - TNT_BASE_OFFSET) + A.offset_[0]; 00368 offset_[1] = (j1 - TNT_BASE_OFFSET) + A.offset_[1]; 00369 dim_[0] = i2-i1 + 1; 00370 dim_[1] = j2-j1+1; 00371 } 00372 00373 const_Region2D<Array2D> operator()(Subscript i1, Subscript i2, 00374 Subscript j1, Subscript j2) 00375 { 00376 #ifdef TNT_BOUNDS_CHECK 00377 assert( i1 <= i2 ); 00378 assert( j1 <= j2); 00379 assert( TNT_BASE_OFFSET <= i1); 00380 assert( i2<= dim_[0] + TNT_BASE_OFFSET-1); 00381 assert( TNT_BASE_OFFSET <= j1); 00382 assert( j2<= dim_[0] + TNT_BASE_OFFSET-1 ); 00383 #endif 00384 return const_Region2D<Array2D>(A_, 00385 i1+offset_[0], offset_[0] + i2, 00386 j1+offset_[1], offset_[1] + j2); 00387 } 00388 00389 00390 const_Region2D<Array2D> operator()(const Index1D &I, 00391 const Index1D &J) 00392 { 00393 #ifdef TNT_BOUNDS_CHECK 00394 assert( I.lbound() <= I.ubound() ); 00395 assert( J.lbound() <= J.ubound() ); 00396 assert( TNT_BASE_OFFSET <= I.lbound()); 00397 assert( I.ubound()<= dim_[0] + TNT_BASE_OFFSET-1); 00398 assert( TNT_BASE_OFFSET <= J.lbound()); 00399 assert( J.ubound() <= dim_[1] + TNT_BASE_OFFSET-1 ); 00400 #endif 00401 00402 return const_Region2D<Array2D>(A_, I.lbound()+offset_[0], 00403 offset_[0] + I.ubound(), offset_[1]+J.lbound(), 00404 offset_[1] + J.ubound()); 00405 } 00406 00407 00408 inline const T & operator() (Subscript i, Subscript j) const 00409 { 00410 #ifdef TNT_BOUNDS_CHECK 00411 assert( TNT_BASE_OFFSET <= i); 00412 assert( i<= dim_[0] + TNT_BASE_OFFSET-1); 00413 assert( TNT_BASE_OFFSET <= j); 00414 assert( j<= dim_[1] + TNT_BASE_OFFSET-1 ); 00415 #endif 00416 return A_(i+offset_[0], j+offset_[1]); 00417 } 00418 00419 }; 00420 00421 00422 // ************** std::ostream algorithms ******************************* 00423 00424 template <class Array2D> 00425 std::ostream& operator<<(std::ostream &s, const const_Region2D<Array2D> &A) 00426 { 00427 Subscript start = A.lbound(); 00428 Subscript Mend=A.lbound()+ A.num_rows() - 1; 00429 Subscript Nend=A.lbound() + A.num_cols() - 1; 00430 00431 00432 s << A.num_rows() << " " << A.num_cols() << "\n"; 00433 for (Subscript i=start; i<=Mend; i++) 00434 { 00435 for (Subscript j=start; j<=Nend; j++) 00436 { 00437 s << A(i,j) << " "; 00438 } 00439 s << "\n"; 00440 } 00441 00442 00443 return s; 00444 } 00445 00446 template <class Array2D> 00447 std::ostream& operator<<(std::ostream &s, const Region2D<Array2D> &A) 00448 { 00449 Subscript start = A.lbound(); 00450 Subscript Mend=A.lbound()+ A.num_rows() - 1; 00451 Subscript Nend=A.lbound() + A.num_cols() - 1; 00452 00453 00454 s << A.num_rows() << " " << A.num_cols() << "\n"; 00455 for (Subscript i=start; i<=Mend; i++) 00456 { 00457 for (Subscript j=start; j<=Nend; j++) 00458 { 00459 s << A(i,j) << " "; 00460 } 00461 s << "\n"; 00462 } 00463 00464 00465 return s; 00466 00467 } 00468 00469 } // namespace TNT 00470 00471 #endif // REGION2D_H 00472