SpatialOps
MatVecFields.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2014 The University of Utah
3  *
4  * Permission is hereby granted, free of charge, to any person obtaining a copy
5  * of this software and associated documentation files (the "Software"), to
6  * deal in the Software without restriction, including without limitation the
7  * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
8  * sell copies of the Software, and to permit persons to whom the Software is
9  * furnished to do so, subject to the following conditions:
10  *
11  * The above copyright notice and this permission notice shall be included in
12  * all copies or substantial portions of the Software.
13  *
14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
19  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
20  * IN THE SOFTWARE.
21  * ----------------------------------------------------------------------------
22 */
23 
24 #ifndef SpatialOps_MatVecFields_h
25 #define SpatialOps_MatVecFields_h
26 #include <cmath>
27 #include <stdexcept>
28 #include <sstream>
29 #include <spatialops/SpatialOpsConfigure.h>
31 #include <spatialops/structured/SpatialFieldStore.h>
32 
46 namespace SpatialOps
47 {
48 
49 template< typename OpT >
50 class VecScalarOp;
51 template< typename OpT >
52 class MatScalarOp;
53 template< typename OpT >
54 class VecVecOp;
55 template< typename OpT >
56 class MatVecOp;
57 template< typename OpT >
58 class MatMatOp;
59 template< typename OpT >
60 class MatUnaryOp;
61 template< typename FieldT >
62 class EigenDecomposition;
63 template< typename FieldT >
64 class LinearSolve;
65 template< typename FieldT >
66 class MatScalDiagOffset;
67 
79 template< typename FieldT >
80 class FieldVector
81 {
82 private:
83  std::vector<SpatFldPtr<FieldT> > vals_;
84  typedef typename FieldT::value_type ValT;
85 public:
86 
97  FieldVector( const int elements,
98  const MemoryWindow window,
99  const BoundaryCellInfo bc,
100  const GhostData ghosts,
101  const short int devIdx = CPU_INDEX )
102  : vals_( )
103  {
104  int i = elements;
105  vals_.reserve ( i );
106  while( i-- ){
107  vals_.push_back ( SpatialFieldStore::get_from_window<FieldT>( window, bc, ghosts, devIdx ) );
108  }
109  }
110 
117  FieldVector( const std::vector<SpatFldPtr<FieldT> >& vec )
118  : vals_( vec )
119  {
120  }
121 
132  template< typename OpT >
133  void operator = ( const VecScalarOp< OpT >& rhsOp );
134 
145  template< typename OpT >
146  void operator = ( const VecVecOp< OpT >& rhsOp );
147 
158  template< typename OpT >
159  void operator = ( const MatVecOp< OpT >& rhsOp );
160 
172  template< typename OpT >
173  void operator = ( const MatUnaryOp< OpT >& rhsOp );
174 
185  void operator = ( const ValT scalar );
186 
187  FieldT& operator() ( const int j )
188  { return *(vals_[j]); }
189  const FieldT& operator() ( const int j ) const
190  { return *(vals_[j]); }
191  FieldT& at ( const int j )
192  { return *(vals_[j]); }
193  const FieldT& at ( const int j ) const
194  { return *(vals_[j]); }
195 
196  std::vector<SpatFldPtr<FieldT> >& raw_vector(void)
197  { return vals_; }
198  const std::vector<SpatFldPtr<FieldT> >& raw_vector(void) const
199  { return vals_; }
200 
201  int elements() const
202  { return vals_.size(); }
203 
225  void add_device(short int deviceIndex)
226  {
227  int i = vals_.size();
228  while( i-- ){
229  vals_[i]->add_device(deviceIndex);
230  }
231  }
243  void validate_device(short int deviceIndex)
244  {
245  int i = vals_.size();
246  while( i-- ){
247  vals_[i]->validate_device(deviceIndex);
248  }
249  }
260  inline void set_device_as_active( const short int deviceIndex )
261  {
262  int i = vals_.size();
263  while( i-- ){
264  vals_[i]->set_device_as_active(deviceIndex);
265  }
266  }
267 
268 #ifdef __CUDACC__
269 
274  inline bool is_valid(int const deviceIndex) const
275  {
276  int i = vals_.size();
277  bool valid = true;
278  while( i-- ){
279  valid &= vals_[i]->is_valid(deviceIndex);
280  }
281  return valid;
282  }
286  inline int active_device_index(void) const
287  {
288  /*Note: Assume the first field has the same active field as all others*/
289  return vals_[0]->active_device_index();
290  }
291 
297  inline void set_stream( const cudaStream_t& stream )
298  {
299  int i = vals_.size();
300  while( i-- ){
301  vals_[i]->set_stream(stream);
302  }
303  }
304 
311  inline void set_stream( const size_t i, const cudaStream_t& stream )
312  {
313  vals_[i]->set_stream(stream);
314  }
315 
321  inline cudaStream_t const & get_stream( const size_t i ) const
322  {
323  return vals_[i]->get_stream();
324  }
325 
331  inline cudaEvent_t const & get_last_event( const size_t i ) const
332  {
333  return vals_[i]->get_last_event();
334  }
335 #endif /* __CUDACC__ */
336 
337 };
338 
352 template< typename FieldT >
353 class FieldMatrix
354 {
355 private:
356  std::vector<SpatFldPtr<FieldT> > vals_;
357  typedef typename FieldT::value_type ValT;
358 public:
359 
369  FieldMatrix( const int elements,
370  const MemoryWindow window,
371  const BoundaryCellInfo bc,
372  const GhostData ghosts,
373  const short int devIdx = CPU_INDEX )
374  : vals_(),
375  elements_( elements )
376  {
377  int ij = elements * elements;
378  vals_.reserve ( ij );
379  while( ij-- ){
380  vals_.push_back ( SpatialFieldStore::get_from_window<FieldT>( window, bc, ghosts, devIdx ) );
381  }
382  }
383 
390  FieldMatrix( const std::vector<SpatFldPtr<FieldT> >& vec )
391  : vals_( vec ),
392  elements_( std::sqrt( vec.size() ) )
393  {
394 #ifndef NDEBUG
395  double elements = std::sqrt( vec.size() );
396  assert( elements == std::floor( elements ) ); //Ensure that the vector coming in could be square.
397 #endif
398  }
399 
408  void add_to_diagonal(const typename FieldT::value_type c);
409 
419  MatUnaryOp< EigenDecomposition<FieldT> > real_eigenvalues();
428  void eigenvalues(FieldVector<FieldT>& realVector,
429  FieldVector<FieldT>& complexVector);
430 
444  MatVecOp< LinearSolve< FieldT > > solve( const FieldVector< FieldT >& vec );
445 
456  template< typename OpT >
457  void operator = ( const MatScalarOp< OpT >& rhsOp );
458 
469  template< typename OpT >
470  void operator = ( const MatMatOp< OpT >& rhsOp );
471 
478  void operator = ( const ValT scalar );
479 
480  FieldT& operator() ( const int i, const int j )
481  { return *vals_[i * elements_ + j]; }
482 
483  FieldT& operator() ( const int ij )
484  { return *(vals_[ij]); }
485 
486  FieldT& at ( const int ij )
487  { return *(vals_[ij]); }
488 
489  const FieldT& operator() ( const int i, const int j ) const
490  { return *vals_[i * elements_ + j]; }
491 
492  const FieldT& operator() ( const int ij ) const
493  { return *(vals_[ij]); }
494 
495  const FieldT& at ( const int ij ) const
496  { return *(vals_[ij]); }
497 
498  FieldVector<FieldT> row (const int i) const
499  {
500  std::vector<SpatFldPtr<FieldT> > rowVector(vals_.begin() + (i * elements_), vals_.begin() + ((i * elements_) + elements_));
501 
502  return FieldVector<FieldT>(rowVector);
503  }
504 
505  std::vector<SpatFldPtr<FieldT> >& raw_vector(void)
506  { return vals_; }
507  const std::vector<SpatFldPtr<FieldT> >& raw_vector(void) const
508  { return vals_; }
509 
510  int elements() const
511  { return elements_; }
512 
534  void add_device(short int deviceIndex)
535  {
536  int i = vals_.size();
537  while( i-- ){
538  vals_[i]->add_device(deviceIndex);
539  }
540  }
552  void validate_device(short int deviceIndex)
553  {
554  int i = vals_.size();
555  while( i-- ){
556  vals_[i]->validate_device(deviceIndex);
557  }
558  }
569  inline void set_device_as_active( const short int deviceIndex )
570  {
571  int i = vals_.size();
572  while( i-- ){
573  vals_[i]->set_device_as_active(deviceIndex);
574  }
575  }
576 
577 #ifdef __CUDACC__
578 
583  inline bool is_valid(int const deviceIndex) const
584  {
585  int i = vals_.size();
586  bool valid = true;
587  while( i-- ){
588  valid &= vals_[i]->is_valid(deviceIndex);
589  }
590  return valid;
591  }
595  inline int active_device_index(void) const
596  {
597  /*Note: Assume the first field has the same active field as all others*/
598  return vals_[0]->active_device_index();
599  }
600 
606  inline void set_stream( const cudaStream_t& stream )
607  {
608  int i = vals_.size();
609  while( i-- ){
610  vals_[i]->set_stream(stream);
611  }
612  }
613 
620  inline void set_stream( const size_t i, const cudaStream_t& stream )
621  {
622  vals_[i]->set_stream(stream);
623  }
624 
630  inline cudaStream_t const & get_stream( const size_t i ) const
631  {
632  return vals_[i]->get_stream();
633  }
634 
640  inline cudaEvent_t const & get_last_event( const size_t i ) const
641  {
642  return vals_[i]->get_last_event();
643  }
644 #endif /* __CUDACC__ */
645 
646 private:
647  const int elements_;
648 };
649 
650 template< typename FieldT >
651 void FieldVector<FieldT>::operator = ( const ValT scalar )
652 {
653  int i = vals_.size();
654  while( i-- ){
655  *(vals_[i]) <<= scalar;
656  }
657 }
658 
659 template< typename FieldT >
660 void FieldMatrix<FieldT>::operator = ( const ValT scalar )
661 {
662  int i = vals_.size();
663  while( i-- ){
664  *(vals_[i]) <<= scalar;
665  }
666 }
667 
668 template< typename FieldT >
669 template< typename OpT >
670 void FieldVector< FieldT >::operator = ( const VecScalarOp< OpT >& rhsOp ){
671  rhsOp.eval( this );
672 }
673 
674 template< typename FieldT >
675 template< typename OpT >
676 void FieldMatrix< FieldT >::operator = ( const MatScalarOp< OpT >& rhsOp ){
677  rhsOp.eval( this );
678 }
679 
680 template< typename FieldT >
681 template< typename OpT >
682 void FieldVector< FieldT >::operator = ( const VecVecOp< OpT >& rhsOp ){
683  rhsOp.eval( this );
684 }
685 
686 template< typename FieldT >
687 template< typename OpT >
688 void FieldVector< FieldT >::operator = ( const MatVecOp< OpT >& rhsOp ){
689  rhsOp.eval( this );
690 }
691 
692 template< typename FieldT >
693 template< typename OpT >
694 void FieldMatrix< FieldT >::operator = ( const MatMatOp< OpT >& rhsOp ){
695  rhsOp.eval( this );
696 }
697 
698 template< typename FieldT >
699 template< typename OpT >
700 void FieldVector< FieldT>::operator = ( const MatUnaryOp< OpT >& rhsOp ){
701  rhsOp.eval( this );
702 }
703 
704 } // namespace SpatialOps
705 
706 #endif // SpatialOps_MatVecField_h
STL namespace.