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