MOOCHO (Single Doxygen Collection) Version of the Day
DenseLinAlgPack_DMatrixClass.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 //
00042 
00043 #include <iomanip>
00044 
00045 #include "DenseLinAlgPack_DMatrixClass.hpp"
00046 
00047 namespace DenseLinAlgPack {
00048 
00049 // ////////////////////////////////////////////////////////////////////////////////
00050 // DMatrixSlice
00051 
00052 DVectorSlice DMatrixSlice::p_diag(difference_type k) const {
00053   if(k > 0) {
00054     validate_col_subscript(k+1);
00055     // upper diagonal (k > 0)
00056     return DVectorSlice( const_cast<value_type*>(col_ptr(1)) + k*max_rows()
00057       , cols()-k > rows() ? rows() : cols()-k, max_rows()+1 );
00058   }
00059   // lower diagonal (k < 0) or center diagonal (k = 0)
00060   validate_row_subscript(-k+1);
00061   return DVectorSlice( const_cast<value_type*>(col_ptr(1)) - k
00062     , rows()+k > cols() ? cols() : rows()+k, max_rows()+1 );
00063 }
00064 
00065 EOverLap DMatrixSlice::overlap(const DMatrixSlice& gms) const
00066 {
00067   typedef DMatrixSlice::size_type size_type;
00068   
00069   const DVectorSlice::value_type
00070     *raw_ptr1 = col_ptr(1),
00071     *raw_ptr2 = gms.col_ptr(1);
00072 
00073   if( !raw_ptr1 || !raw_ptr2 )
00074     return NO_OVERLAP;    // If one of the views is unbound there can't be any overlap
00075 
00076   DVectorSlice::size_type
00077     max_rows1 = max_rows(),
00078     max_rows2 = gms.max_rows(),
00079     rows1 = rows(),
00080     rows2 = gms.rows(),
00081     cols1 = cols(),
00082     cols2 = gms.cols();
00083 
00084   // Establish a frame of reference where raw_ptr1 < raw_ptr2
00085   if(raw_ptr1 > raw_ptr2) {
00086     std::swap(raw_ptr1,raw_ptr2);
00087     std::swap(max_rows1,max_rows2);
00088     std::swap(rows1,rows2);
00089     std::swap(cols1,cols2);
00090   }
00091 
00092   if( raw_ptr2 > (raw_ptr1 + (cols1 - 1) * max_rows1 + (rows1 - 1)) ) {
00093     return NO_OVERLAP; // can't be any overlap
00094   }
00095 
00096   DVectorSlice::size_type
00097     start1 = 0,
00098     start2 = raw_ptr2 - raw_ptr1;
00099 
00100   if(start1 == start2 && max_rows1 == max_rows2 && rows1 == rows2 && cols1 == cols2)
00101     return SAME_MEM;
00102   if(start1 + (rows1 - 1) + (cols1 - 1) * max_rows1 < start2)
00103     return NO_OVERLAP;  // start2 is past the last element in m1 so no overlap.
00104   // There may be some overlap so determine if start2 lays in the region of m1.
00105   // Determine the number of elements in v that start2 is past the start of a
00106   // column of m1.  If start2 was the first element in one of m1's cols
00107   // row_i would be 1, and if it was just before for first element of the next
00108   // column of m1 then row_i would equal to max_rows1.
00109   size_type row_i = (start2 - start1 + 1) % max_rows1;
00110   if(row_i <= rows1)
00111     return SOME_OVERLAP; // start2 is in a column of m1 so there is some overlap
00112   // Determine how many rows in the original matrix are below the last row in m1
00113   size_type lower_rows = max_rows1 - (start1 % max_rows1 + rows1);
00114   if(row_i < rows1 + lower_rows)
00115     return NO_OVERLAP; // m2 lays below m1 in the original matrix
00116   // If you get here start2 lays above m1 in original matix so if m2 has enough
00117   // rows then the lower rows of m2 will overlap the upper rows of m1.
00118   if(row_i + rows2 - 1 <= max_rows1)
00119     return NO_OVERLAP; // m2 completely lays above m1
00120   return SOME_OVERLAP; // Some lower rows of m2 extend into m1 
00121 }
00122 
00123 #ifdef LINALGPACK_CHECK_RANGE
00124 void DMatrixSlice::validate_row_subscript(size_type i) const
00125 {
00126   if( i > rows() || !i )
00127     throw std::out_of_range( "DMatrixSlice::validate_row_subscript(i) :"
00128                   "row index i is out of bounds"        );
00129 }
00130 #endif
00131 
00132 #ifdef LINALGPACK_CHECK_RANGE
00133 void DMatrixSlice::validate_col_subscript(size_type j) const
00134 {
00135   if( j > cols() || !j )
00136     throw std::out_of_range( "DMatrixSlice::validate_col_subscript(j) :"
00137                   "column index j is out of bounds"     );
00138 }
00139 #endif
00140 
00141 #ifdef LINALGPACK_CHECK_SLICE_SETUP
00142 void DMatrixSlice::validate_setup(size_type size) const
00143 {
00144   if( !ptr_ && !rows() && !cols() && !max_rows() )
00145       return; // an unsized matrix slice is ok.
00146   if( (rows() - 1) + (cols() - 1) * max_rows() + 1 > size )
00147     throw std::out_of_range( "DMatrixSlice::validate_setup() : "
00148                   " DMatrixSlice constructed that goes past end of array" );
00149 }
00150 #endif
00151 
00152 // /////////////////////////////////////////////////////////////////////////////////
00153 // DMatrix
00154 
00155 DVectorSlice DMatrix::p_diag(difference_type k) const { 
00156   if(k > 0) {
00157     validate_col_subscript(k+1);
00158     // upper diagonal (k > 0)
00159     return DVectorSlice( const_cast<value_type*>(col_ptr(1)) + k * rows()
00160       , cols()-k > rows() ? rows() : cols()-k, rows()+1 );
00161   }
00162   // lower diagonal (k < 0) or center diagonal (k = 0)
00163   validate_row_subscript(-k+1);
00164   return DVectorSlice( const_cast<value_type*>(col_ptr(1)) - k
00165     , rows()+k > cols() ? cols() : rows()+k, rows()+1 );
00166 }
00167 
00168 EOverLap DMatrix::overlap(const DMatrixSlice& gms) const {
00169   return (*this)().overlap(gms);
00170 }
00171 
00172 #ifdef LINALGPACK_CHECK_RANGE
00173 void DMatrix::validate_row_subscript(size_type i) const {
00174   if( i > rows() || !i )
00175     throw std::out_of_range("DMatrix::validate_row_subscript(i) : row index out of bounds");
00176 }
00177 #endif
00178 
00179 #ifdef LINALGPACK_CHECK_RANGE
00180 void DMatrix::validate_col_subscript(size_type j) const {
00181   if( j > cols() || !j )
00182     throw std::out_of_range("DMatrix::validate_col_subscript(j) : column index out of bounds");
00183 }
00184 #endif
00185 
00186 } // end namespace DenseLinAlgPack
00187 
00188 // ///////////////////////////////////////////////////////////////////////////////
00189 // Non-member funcitons
00190 
00191 void DenseLinAlgPack::assert_gms_sizes(const DMatrixSlice& gms1, BLAS_Cpp::Transp trans1
00192   , const DMatrixSlice& gms2, BLAS_Cpp::Transp trans2)
00193 {
00194   if  (
00195       (trans1 == trans2) ? 
00196         gms1.rows() == gms2.rows() && gms1.cols() == gms2.cols() 
00197         : gms1.rows() == gms2.cols() && gms1.cols() == gms2.rows()
00198     )
00199     return; // compatible sizes so exit
00200   // not compatible sizes.
00201   throw std::length_error("Matrix sizes are not the compatible");
00202 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines