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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 //
00029 
00030 #include <iomanip>
00031 
00032 #include "DenseLinAlgPack_DMatrixClass.hpp"
00033 
00034 namespace DenseLinAlgPack {
00035 
00036 // ////////////////////////////////////////////////////////////////////////////////
00037 // DMatrixSlice
00038 
00039 DVectorSlice DMatrixSlice::p_diag(difference_type k) const {
00040   if(k > 0) {
00041     validate_col_subscript(k+1);
00042     // upper diagonal (k > 0)
00043     return DVectorSlice( const_cast<value_type*>(col_ptr(1)) + k*max_rows()
00044       , cols()-k > rows() ? rows() : cols()-k, max_rows()+1 );
00045   }
00046   // lower diagonal (k < 0) or center diagonal (k = 0)
00047   validate_row_subscript(-k+1);
00048   return DVectorSlice( const_cast<value_type*>(col_ptr(1)) - k
00049     , rows()+k > cols() ? cols() : rows()+k, max_rows()+1 );
00050 }
00051 
00052 EOverLap DMatrixSlice::overlap(const DMatrixSlice& gms) const
00053 {
00054   typedef DMatrixSlice::size_type size_type;
00055   
00056   const DVectorSlice::value_type
00057     *raw_ptr1 = col_ptr(1),
00058     *raw_ptr2 = gms.col_ptr(1);
00059 
00060   if( !raw_ptr1 || !raw_ptr2 )
00061     return NO_OVERLAP;    // If one of the views is unbound there can't be any overlap
00062 
00063   DVectorSlice::size_type
00064     max_rows1 = max_rows(),
00065     max_rows2 = gms.max_rows(),
00066     rows1 = rows(),
00067     rows2 = gms.rows(),
00068     cols1 = cols(),
00069     cols2 = gms.cols();
00070 
00071   // Establish a frame of reference where raw_ptr1 < raw_ptr2
00072   if(raw_ptr1 > raw_ptr2) {
00073     std::swap(raw_ptr1,raw_ptr2);
00074     std::swap(max_rows1,max_rows2);
00075     std::swap(rows1,rows2);
00076     std::swap(cols1,cols2);
00077   }
00078 
00079   if( raw_ptr2 > (raw_ptr1 + (cols1 - 1) * max_rows1 + (rows1 - 1)) ) {
00080     return NO_OVERLAP; // can't be any overlap
00081   }
00082 
00083   DVectorSlice::size_type
00084     start1 = 0,
00085     start2 = raw_ptr2 - raw_ptr1;
00086 
00087   if(start1 == start2 && max_rows1 == max_rows2 && rows1 == rows2 && cols1 == cols2)
00088     return SAME_MEM;
00089   if(start1 + (rows1 - 1) + (cols1 - 1) * max_rows1 < start2)
00090     return NO_OVERLAP;  // start2 is past the last element in m1 so no overlap.
00091   // There may be some overlap so determine if start2 lays in the region of m1.
00092   // Determine the number of elements in v that start2 is past the start of a
00093   // column of m1.  If start2 was the first element in one of m1's cols
00094   // row_i would be 1, and if it was just before for first element of the next
00095   // column of m1 then row_i would equal to max_rows1.
00096   size_type row_i = (start2 - start1 + 1) % max_rows1;
00097   if(row_i <= rows1)
00098     return SOME_OVERLAP; // start2 is in a column of m1 so there is some overlap
00099   // Determine how many rows in the original matrix are below the last row in m1
00100   size_type lower_rows = max_rows1 - (start1 % max_rows1 + rows1);
00101   if(row_i < rows1 + lower_rows)
00102     return NO_OVERLAP; // m2 lays below m1 in the original matrix
00103   // If you get here start2 lays above m1 in original matix so if m2 has enough
00104   // rows then the lower rows of m2 will overlap the upper rows of m1.
00105   if(row_i + rows2 - 1 <= max_rows1)
00106     return NO_OVERLAP; // m2 completely lays above m1
00107   return SOME_OVERLAP; // Some lower rows of m2 extend into m1 
00108 }
00109 
00110 #ifdef LINALGPACK_CHECK_RANGE
00111 void DMatrixSlice::validate_row_subscript(size_type i) const
00112 {
00113   if( i > rows() || !i )
00114     throw std::out_of_range( "DMatrixSlice::validate_row_subscript(i) :"
00115                   "row index i is out of bounds"        );
00116 }
00117 #endif
00118 
00119 #ifdef LINALGPACK_CHECK_RANGE
00120 void DMatrixSlice::validate_col_subscript(size_type j) const
00121 {
00122   if( j > cols() || !j )
00123     throw std::out_of_range( "DMatrixSlice::validate_col_subscript(j) :"
00124                   "column index j is out of bounds"     );
00125 }
00126 #endif
00127 
00128 #ifdef LINALGPACK_CHECK_SLICE_SETUP
00129 void DMatrixSlice::validate_setup(size_type size) const
00130 {
00131   if( !ptr_ && !rows() && !cols() && !max_rows() )
00132       return; // an unsized matrix slice is ok.
00133   if( (rows() - 1) + (cols() - 1) * max_rows() + 1 > size )
00134     throw std::out_of_range( "DMatrixSlice::validate_setup() : "
00135                   " DMatrixSlice constructed that goes past end of array" );
00136 }
00137 #endif
00138 
00139 // /////////////////////////////////////////////////////////////////////////////////
00140 // DMatrix
00141 
00142 DVectorSlice DMatrix::p_diag(difference_type k) const { 
00143   if(k > 0) {
00144     validate_col_subscript(k+1);
00145     // upper diagonal (k > 0)
00146     return DVectorSlice( const_cast<value_type*>(col_ptr(1)) + k * rows()
00147       , cols()-k > rows() ? rows() : cols()-k, rows()+1 );
00148   }
00149   // lower diagonal (k < 0) or center diagonal (k = 0)
00150   validate_row_subscript(-k+1);
00151   return DVectorSlice( const_cast<value_type*>(col_ptr(1)) - k
00152     , rows()+k > cols() ? cols() : rows()+k, rows()+1 );
00153 }
00154 
00155 EOverLap DMatrix::overlap(const DMatrixSlice& gms) const {
00156   return (*this)().overlap(gms);
00157 }
00158 
00159 #ifdef LINALGPACK_CHECK_RANGE
00160 void DMatrix::validate_row_subscript(size_type i) const {
00161   if( i > rows() || !i )
00162     throw std::out_of_range("DMatrix::validate_row_subscript(i) : row index out of bounds");
00163 }
00164 #endif
00165 
00166 #ifdef LINALGPACK_CHECK_RANGE
00167 void DMatrix::validate_col_subscript(size_type j) const {
00168   if( j > cols() || !j )
00169     throw std::out_of_range("DMatrix::validate_col_subscript(j) : column index out of bounds");
00170 }
00171 #endif
00172 
00173 } // end namespace DenseLinAlgPack
00174 
00175 // ///////////////////////////////////////////////////////////////////////////////
00176 // Non-member funcitons
00177 
00178 void DenseLinAlgPack::assert_gms_sizes(const DMatrixSlice& gms1, BLAS_Cpp::Transp trans1
00179   , const DMatrixSlice& gms2, BLAS_Cpp::Transp trans2)
00180 {
00181   if  (
00182       (trans1 == trans2) ? 
00183         gms1.rows() == gms2.rows() && gms1.cols() == gms2.cols() 
00184         : gms1.rows() == gms2.cols() && gms1.cols() == gms2.rows()
00185     )
00186     return; // compatible sizes so exit
00187   // not compatible sizes.
00188   throw std::length_error("Matrix sizes are not the compatible");
00189 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:10:58 2011 for MOOCHO (Single Doxygen Collection) by  doxygen 1.6.3