AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects Version of the Day
AbstractLinAlgPack_GenPermMatrixSliceOp.cpp
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 #include <assert.h>
00043 
00044 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
00045 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00046 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00047 #include "DenseLinAlgPack_DVectorClass.hpp"
00048 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
00049 #include "DenseLinAlgPack_AssertOp.hpp"
00050 
00051 void AbstractLinAlgPack::V_StMtV(
00052     SpVector* y, value_type a, const GenPermMatrixSlice& P
00053   , BLAS_Cpp::Transp P_trans, const DVectorSlice& x )
00054 {
00055   using BLAS_Cpp::no_trans;
00056   using BLAS_Cpp::trans;
00057   namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
00058   using DenseLinAlgPack::MtV_assert_sizes;
00059 
00060   MtV_assert_sizes( P.rows(), P.cols(), P_trans, x.dim() );
00061 
00062   y->resize( BLAS_Cpp::rows( P.rows(), P.cols(), P_trans ), P.nz() );
00063 
00064   typedef SpVector::element_type ele_t;
00065 
00066   if( P.is_identity() ) {
00067     for( size_type i = 1; i <= P.nz(); ++i ) {
00068       const value_type x_j = x(i);
00069       if( x_j != 0.0 )
00070         y->add_element( ele_t( i, a * x_j ) );
00071     }
00072   }   
00073   else if( P_trans == no_trans ) {
00074     for( GenPermMatrixSlice::const_iterator itr = P.begin(); itr != P.end(); ++itr ) {
00075       const size_type
00076         i = itr->row_i(),
00077         j = itr->col_j();
00078       const value_type x_j = x(j);
00079       if( x_j != 0.0 )
00080         y->add_element( ele_t( i, a * x_j ) );
00081     }
00082   }
00083   else {
00084     for( GenPermMatrixSlice::const_iterator itr = P.begin(); itr != P.end(); ++itr ) {
00085       const size_type
00086         j = itr->row_i(),
00087         i = itr->col_j();
00088       const value_type x_j = x(j);
00089       if( x_j != 0.0 )
00090         y->add_element( ele_t( i, a * x_j ) );
00091     }
00092   }
00093   if( P.ordered_by() == GPMSIP::BY_ROW_AND_COL
00094     || ( P_trans == no_trans  && P.ordered_by() == GPMSIP::BY_ROW )
00095     ||  ( P_trans == trans    && P.ordered_by() == GPMSIP::BY_COL ) )
00096   {
00097     y->assume_sorted(true);
00098   }
00099 }
00100 
00101 void AbstractLinAlgPack::V_StMtV(
00102     SpVector* y, value_type a, const GenPermMatrixSlice& P
00103   , BLAS_Cpp::Transp P_trans, const SpVectorSlice& x )
00104 {
00105   using BLAS_Cpp::no_trans;
00106   using BLAS_Cpp::trans;
00107   namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
00108   using DenseLinAlgPack::MtV_assert_sizes;
00109   MtV_assert_sizes( P.rows(), P.cols(), P_trans, x.dim() );
00110 
00111   y->resize( BLAS_Cpp::rows( P.rows(), P.cols(), P_trans ), P.nz() );
00112 
00113   typedef SpVector::element_type ele_t;
00114   const SpVectorSlice::element_type *ele_ptr;
00115 
00116   if( P.is_identity() ) {
00117     AbstractLinAlgPack::add_elements( y, 1.0, x(1,P.nz()) );
00118     AbstractLinAlgPack::Vt_S( &(*y)(), a );
00119   }   
00120   else if( x.is_sorted() ) {
00121     if( P_trans == no_trans ) {
00122       for( GenPermMatrixSlice::const_iterator itr = P.begin(); itr != P.end(); ++itr ) {
00123         const size_type
00124           i = itr->row_i(),
00125           j = itr->col_j();
00126         if( ele_ptr = x.lookup_element(j) )
00127           y->add_element( ele_t( i, a * ele_ptr->value() ) );
00128       }
00129     }
00130     else {
00131       for( GenPermMatrixSlice::const_iterator itr = P.begin(); itr != P.end(); ++itr ) {
00132         const size_type
00133           j = itr->row_i(),
00134           i = itr->col_j();
00135         if( ele_ptr = x.lookup_element(j) )
00136           y->add_element( ele_t( i, a * ele_ptr->value() ) );
00137       }
00138     }
00139   }
00140   else {
00141     TEUCHOS_TEST_FOR_EXCEPT(true);  // ToDo: Implement the other cases!
00142   }
00143   if(  P.ordered_by() == GPMSIP::BY_ROW_AND_COL
00144     ||  ( P_trans == no_trans && P.ordered_by() == GPMSIP::BY_ROW )
00145     ||  ( P_trans == trans    && P.ordered_by() == GPMSIP::BY_COL ) )
00146   {
00147     y->assume_sorted(true);
00148   }
00149 }
00150 
00151 void AbstractLinAlgPack::Vp_StMtV(
00152     SpVector* y, value_type a, const GenPermMatrixSlice& P
00153   , BLAS_Cpp::Transp P_trans, const DVectorSlice& x )
00154 {
00155   using BLAS_Cpp::no_trans;
00156   using BLAS_Cpp::trans;
00157   namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
00158   using DenseLinAlgPack::Vp_MtV_assert_sizes;
00159 
00160   Vp_MtV_assert_sizes( y->dim(), P.rows(), P.cols(), P_trans, x.dim() );
00161 
00162   typedef SpVector::element_type ele_t;
00163 
00164   const bool was_sorted = y->is_sorted();
00165 
00166   if( P.is_identity() ) {
00167     for( size_type i = 1; i <= P.nz(); ++i ) {
00168       const value_type x_j = x(i);
00169       if( x_j != 0.0 )
00170         y->add_element( ele_t( i, a * x_j ) );
00171     }
00172   }   
00173   else if( P_trans == no_trans ) {
00174     for( GenPermMatrixSlice::const_iterator itr = P.begin(); itr != P.end(); ++itr ) {
00175       const size_type
00176         i = itr->row_i(),
00177         j = itr->col_j();
00178       const value_type x_j = x(j);
00179       if( x_j != 0.0 )
00180         y->add_element( ele_t( i, a * x_j ) );
00181     }
00182   }
00183   else {
00184     for( GenPermMatrixSlice::const_iterator itr = P.begin(); itr != P.end(); ++itr ) {
00185       const size_type
00186         j = itr->row_i(),
00187         i = itr->col_j();
00188       const value_type x_j = x(j);
00189       if( x_j != 0.0 )
00190         y->add_element( ele_t( i, a * x_j ) );
00191     }
00192   }
00193   if( was_sorted &&
00194     ( P.ordered_by() == GPMSIP::BY_ROW_AND_COL
00195       || ( P_trans == no_trans  && P.ordered_by() == GPMSIP::BY_ROW )
00196       ||  ( P_trans == trans    && P.ordered_by() == GPMSIP::BY_COL ) ) )
00197   {
00198     y->assume_sorted(true);
00199   }
00200 }
00201 
00202 void AbstractLinAlgPack::Vp_StMtV(
00203     DVectorSlice* y, value_type a, const GenPermMatrixSlice& P
00204   , BLAS_Cpp::Transp P_trans, const DVectorSlice& x, value_type b )
00205 {
00206   using DenseLinAlgPack::Vt_S;
00207   using DenseLinAlgPack::Vp_MtV_assert_sizes;
00208   Vp_MtV_assert_sizes( y->dim(), P.rows(), P.cols(), P_trans, x.dim() );
00209   // y = b*y
00210   if( b == 0.0 )
00211     *y = 0.0;
00212   else
00213     Vt_S(y,b);  
00214   // y += a*op(P)*x
00215   if( P.is_identity() ) {
00216     if( b == 0.0 )
00217       *y = 0.0;
00218     else
00219       DenseLinAlgPack::Vt_S( y, b );
00220     DenseLinAlgPack::Vp_StV( &(*y)(1,P.nz()), a, x(1,P.nz()) );
00221   }   
00222   else if( P_trans == BLAS_Cpp::no_trans ) {
00223     for( GenPermMatrixSlice::const_iterator itr = P.begin(); itr != P.end(); ++itr ) {
00224       const size_type
00225         i = itr->row_i(),
00226         j = itr->col_j();
00227       (*y)(i) += a * x(j);
00228     }
00229   }
00230   else {
00231     for( GenPermMatrixSlice::const_iterator itr = P.begin(); itr != P.end(); ++itr ) {
00232       const size_type
00233         j = itr->row_i(),
00234         i = itr->col_j();
00235       (*y)(i) += a * x(j);
00236     }
00237   }
00238 }
00239 
00240 void AbstractLinAlgPack::Vp_StMtV(
00241     DVectorSlice* y, value_type a, const GenPermMatrixSlice& P
00242   , BLAS_Cpp::Transp P_trans, const SpVectorSlice& x, value_type b )
00243 {
00244   using BLAS_Cpp::no_trans;
00245   using BLAS_Cpp::trans;
00246   using BLAS_Cpp::rows;
00247   using BLAS_Cpp::cols;
00248   namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
00249   using DenseLinAlgPack::Vt_S;
00250   using DenseLinAlgPack::Vp_MtV_assert_sizes;
00251   
00252   Vp_MtV_assert_sizes( y->dim(), P.rows(), P.cols(), P_trans, x.dim() );
00253   // y = b*y
00254   if( b == 0.0 )
00255     *y = 0.0;
00256   else
00257     Vt_S(y,b);
00258   // y += a*op(P)*x
00259   if( P.is_identity() ) {
00260     DenseLinAlgPack::Vt_S( y, b ); // takes care of b == 0.0 and y == NaN
00261     AbstractLinAlgPack::Vp_StV( &(*y)(1,P.nz()), a, x(1,P.nz()) );
00262   }   
00263   else if( x.is_sorted() ) {
00264     const SpVectorSlice::difference_type x_off = x.offset();
00265     if( P_trans == no_trans && P.ordered_by() == GPMSIP::BY_COL ) {
00266       TEUCHOS_TEST_FOR_EXCEPT(true);  // ToDo: implement this!
00267     }
00268     else if( ( P_trans == trans && P.ordered_by() == GPMSIP::BY_ROW )
00269       || P.ordered_by() == GPMSIP::BY_ROW_AND_COL )
00270     {
00271       GenPermMatrixSlice::const_iterator
00272         P_itr = P.begin(),
00273         P_end = P.end();
00274       SpVectorSlice::const_iterator
00275         x_itr = x.begin(),
00276         x_end = x.end();
00277       while( P_itr != P_end && x_itr != x_end ) {
00278         const size_type
00279           i = rows(P_itr->row_i(),P_itr->col_j(),P_trans),
00280           j = cols(P_itr->row_i(),P_itr->col_j(),P_trans);
00281         if( j < x_itr->index() + x_off ) {
00282           ++P_itr;
00283           continue;
00284         }
00285         else if( j > x_itr->index() + x_off ) {
00286           ++x_itr;
00287           continue;
00288         }
00289         else {  // they are equal
00290           (*y)(i) += a * x_itr->value();
00291           ++P_itr;
00292           ++x_itr;
00293         }
00294       }
00295     }
00296     else {
00297       TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement what ever this case is? 
00298     }
00299   }
00300   else {
00301     // Since things do not match up we will have to create a
00302     // temporary dense copy of x to operate on.
00303     TEUCHOS_TEST_FOR_EXCEPT(true);  // ToDo: Implement this!
00304   }
00305 }
00306 
00307 namespace {
00308 
00309 AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::EOrderedBy
00310 ordered_by(
00311   AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::EOrderedBy P_ordered_by
00312   , BLAS_Cpp::Transp P_trans
00313   )
00314 {
00315   using BLAS_Cpp::no_trans;
00316   namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
00317   GPMSIP::EOrderedBy
00318     opP_ordered_by;
00319   switch( P_ordered_by ) {
00320       case GPMSIP::BY_ROW_AND_COL:
00321       opP_ordered_by = GPMSIP::BY_ROW_AND_COL;
00322       break;
00323       case GPMSIP::BY_ROW:
00324       opP_ordered_by = P_trans == no_trans ? GPMSIP::BY_ROW : GPMSIP::BY_COL;
00325       break;
00326       case GPMSIP::BY_COL:
00327       opP_ordered_by = P_trans == no_trans ? GPMSIP::BY_COL : GPMSIP::BY_COL;
00328       break;
00329       case GPMSIP::UNORDERED:
00330       opP_ordered_by = GPMSIP::UNORDERED;
00331       break;
00332         default:
00333       TEUCHOS_TEST_FOR_EXCEPT(true); // Should never happen
00334   }
00335   return opP_ordered_by;
00336 }
00337 
00338 } // end namespace
00339 
00340 void AbstractLinAlgPack::intersection(
00341   const GenPermMatrixSlice     &P1
00342   ,BLAS_Cpp::Transp            P1_trans
00343   ,const GenPermMatrixSlice    &P2
00344   ,BLAS_Cpp::Transp            P2_trans
00345   ,size_type                   *Q_nz
00346   ,const size_type             Q_max_nz
00347   ,size_type                   Q_row_i[]
00348   ,size_type                   Q_col_j[]
00349   ,GenPermMatrixSlice          *Q
00350   )
00351 {
00352   using BLAS_Cpp::no_trans;
00353   using BLAS_Cpp::trans;
00354   using BLAS_Cpp::trans_not;
00355   using BLAS_Cpp::rows;
00356   using BLAS_Cpp::cols;
00357   namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
00358   //
00359   // Q = op(P1)*op(P2)
00360   //
00361   // There are several different possibilities for how to compute this
00362   // intersection.
00363   //
00364   DenseLinAlgPack::MtM_assert_sizes(
00365     P1.rows(), P1.cols() , P1_trans, P2.rows(), P2.cols() , P2_trans );
00366   //
00367   const size_type
00368     opP1_rows = rows(P1.rows(),P1.cols(),P1_trans),
00369     opP1_cols = cols(P1.rows(),P1.cols(),P1_trans),
00370     opP2_rows = rows(P2.rows(),P2.cols(),P2_trans),
00371     opP2_cols = cols(P2.rows(),P2.cols(),P2_trans);
00372   GPMSIP::EOrderedBy
00373     opP1_ordered_by = ordered_by(P1.ordered_by(),P1_trans),
00374     opP2_ordered_by = ordered_by(P2.ordered_by(),P2_trans);
00375   //
00376   *Q_nz = 0;
00377   // Either is zero?
00378   if( !P1.nz() || !P2.nz() ) {
00379     *Q_nz = 0;
00380     if(Q)
00381       Q->initialize(opP1_rows,opP2_cols,GenPermMatrixSlice::ZERO_MATRIX);
00382     return;
00383   }
00384   // Both are identity?
00385   if( P1.is_identity() && P2.is_identity() ) {
00386     *Q_nz = P1.nz(); // Should be the same as P2.nz();
00387     TEUCHOS_TEST_FOR_EXCEPT( !(  P1.nz() == P2.nz()  ) );
00388     if(Q)
00389       Q->initialize(opP1_rows,opP2_cols,GenPermMatrixSlice::IDENTITY_MATRIX);
00390     return;
00391   }
00392   // One is identity?
00393   if( P1.is_identity() || P2.is_identity() ) {
00394     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement this but it is a little tricking?
00395     return;
00396   }
00397   //
00398   // Both are not identity or zero
00399   //
00400   if( ( opP1_ordered_by == GPMSIP::BY_COL || opP1_ordered_by == GPMSIP::BY_ROW_AND_COL )
00401     && ( opP2_ordered_by == GPMSIP::BY_ROW || opP2_ordered_by == GPMSIP::BY_ROW_AND_COL ) )
00402   {
00403     // This is great!  Both of the matrices are sorted and we don't need any temparary storage
00404     if( Q_max_nz ) {
00405       TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement initializing Q_row_i, Q_col_j
00406     }
00407     else {
00408       GenPermMatrixSlice::const_iterator
00409         P1_itr = P1.begin(), // Should not throw exception!
00410         P1_end = P1.end(),
00411         P2_itr = P2.begin(), // Should not throw exception!
00412         P2_end = P2.end();
00413       while( P1_itr != P1_end && P2_itr != P2_end ) {
00414         const size_type
00415           opP1_col_j = cols(P1_itr->row_i(),P1_itr->col_j(),P1_trans),
00416           opP2_row_i = rows(P2_itr->row_i(),P2_itr->col_j(),P2_trans);
00417         if( opP1_col_j < opP2_row_i ) {
00418           ++P1_itr;
00419           continue;
00420         }
00421         if( opP1_col_j > opP2_row_i ) {
00422           ++P2_itr;
00423           continue;
00424         }
00425         ++(*Q_nz);
00426         ++P1_itr;
00427         ++P2_itr;
00428       }
00429     }
00430   }
00431   else {
00432     //
00433     // We will load the row indices of op(P1) or the column indices op(P1)
00434     // indexed by column or row indices (whichever is smaller)
00435     // into a temporary sorted buffer and then loop through the nonzeros of the other.
00436     //
00437     // First let's get reorder P1 and P2 so that we put the rows of P2 into a buffer
00438     //
00439     const GenPermMatrixSlice
00440       &oP1      = opP1_cols > opP2_rows ? P1 : P2,
00441       &oP2      = opP1_cols > opP2_rows ? P2 : P1;
00442     const BLAS_Cpp::Transp
00443       oP1_trans = opP1_cols > opP2_rows ? P1_trans : trans_not(P1_trans),
00444       oP2_trans = opP1_cols > opP2_rows ? P2_trans : trans_not(P2_trans);
00445     // Load the column indices of op(op(P2)) into a lookup array
00446     typedef std::vector<size_type> oP2_col_j_lookup_t;      // Todo: use tempoarary workspace
00447     oP2_col_j_lookup_t oP2_col_j_lookup(rows(oP2.rows(),oP2.rows(),oP2_trans));
00448     std::fill( oP2_col_j_lookup.begin(), oP2_col_j_lookup.end(), 0 );
00449     {
00450       GenPermMatrixSlice::const_iterator
00451         itr     = oP2.begin(), // Should not throw exception!
00452         itr_end = oP2.end();
00453       while( itr != itr_end ) {
00454         oP2_col_j_lookup[rows(itr->row_i(),itr->col_j(),oP2_trans)]
00455           = cols(itr->row_i(),itr->col_j(),oP2_trans);
00456       }
00457     }
00458     // Loop through the columns of op(op(P1)) and look for matches
00459     if( Q_max_nz ) {
00460       TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement initializing Q_row_i, Q_col_j
00461     }
00462     else {
00463       GenPermMatrixSlice::const_iterator
00464         itr     = oP1.begin(), // Should not throw exception!
00465         itr_end = oP1.end();
00466       while( itr != itr_end ) {
00467         const size_type
00468           oP2_col_j = oP2_col_j_lookup[cols(itr->row_i(),itr->col_j(),oP1_trans)];
00469         if(oP2_col_j)
00470           ++(*Q_nz);
00471       }
00472     }
00473 
00474   }
00475   // Setup Q
00476   TEUCHOS_TEST_FOR_EXCEPT( !( Q == NULL ) ); // ToDo: Implement initializing Q
00477 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends