MOOCHO (Single Doxygen Collection) Version of the Day
AbstractLinAlgPack_MatrixSymDiagonalStd.cpp
Go to the documentation of this file.
00001 #if 0
00002 
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00007 //                  Copyright (2003) Sandia Corporation
00008 // 
00009 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00010 // license for use of this work by or on behalf of the U.S. Government.
00011 // 
00012 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission.
00026 //
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //
00039 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00040 // 
00041 // ***********************************************************************
00042 // @HEADER
00043 
00044 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixSymDiagStd.hpp"
00045 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_SpVectorClass.hpp"
00046 #include "DenseLinAlgPack_DVectorOp.hpp"
00047 #include "DenseLinAlgPack_DMatrixOp.hpp"
00048 #include "DenseLinAlgPack_AssertOp.hpp"
00049 
00050 namespace AbstractLinAlgPack {
00051 
00052 MatrixSymDiagStd::MatrixSymDiagStd()
00053 {}
00054 
00055 DVectorSlice MatrixSymDiagStd::diag()
00056 {
00057   return diag_();
00058 }
00059 
00060 const DVectorSlice MatrixSymDiagStd::diag() const
00061 {
00062   return diag_();
00063 }
00064 
00065 // Overridden from MatrixSymInitDiag
00066 
00067 void MatrixSymDiagStd::init_identity( size_type n, value_type alpha = 1.0 )
00068 {
00069   diag_.resize(n);
00070   if(n)
00071     diag_ = alpha;
00072 }
00073 
00074 void MatrixSymDiagStd::init_diagonal( const DVectorSlice& diag )
00075 {
00076   diag_ = diag;
00077 }
00078 
00079 // Overridden from Matrix
00080 
00081 size_type MatrixSymDiagStd::rows() const
00082 {
00083   return diag_.size();
00084 }
00085 
00086 // Overridden from MatrixOp
00087 
00088 void MatrixSymDiagStd::Mp_StM(
00089   DMatrixSlice* gms_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs) const
00090 {
00091   using DenseLinAlgPack::Vp_StV;
00092   // Assert that the dimensions match up.
00093   DenseLinAlgPack::Mp_M_assert_sizes( gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
00094     , rows(), cols(), trans_rhs );
00095   // Add to the diagonal
00096   Vp_StV( &gms_lhs->diag(), alpha, diag_ );
00097 }
00098 
00099 void MatrixSymDiagStd::Vp_StMtV(
00100   DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00101   , const DVectorSlice& vs_rhs2, value_type beta) const
00102 {
00103   const DVectorSlice diag = this->diag();
00104   size_type n = diag.size();
00105 
00106   //
00107   // y = b*y + a * op(A) * x
00108   //
00109   DenseLinAlgPack::Vp_MtV_assert_sizes(
00110     vs_lhs->size(), n, n, trans_rhs1, vs_rhs2.size() );
00111   //
00112   // A is symmetric and diagonal A = diag(diag) so:
00113   //
00114   // y(j) += a * diag(j) * x(j), for j = 1...n
00115   //
00116   if( vs_rhs2.stride() == 1 && vs_lhs->stride() == 1 ) {
00117     // Optimized implementation
00118     const value_type
00119       *d_itr      = diag.raw_ptr(),
00120       *x_itr      = vs_rhs2.raw_ptr();
00121     value_type
00122       *y_itr      = vs_lhs->raw_ptr(),
00123       *y_end      = y_itr + vs_lhs->size();
00124 
00125     if( beta == 0.0 ) {
00126       while( y_itr != y_end )
00127         *y_itr++ = alpha * (*d_itr++) * (*x_itr++);
00128     }
00129     else if( beta == 1.0 ) {
00130       while( y_itr != y_end )
00131         *y_itr++ += alpha * (*d_itr++) * (*x_itr++);
00132     }
00133     else {
00134       for( ; y_itr != y_end; ++y_itr )
00135         *y_itr = beta * (*y_itr) + alpha * (*d_itr++) * (*x_itr++);
00136     }
00137   }
00138   else {
00139     // Generic implementation
00140     DVectorSlice::const_iterator
00141       d_itr = diag.begin(),
00142       x_itr = vs_rhs2.begin();
00143     DVectorSlice::iterator
00144       y_itr = vs_lhs->begin(),
00145       y_end = vs_lhs->end();
00146     for( ; y_itr != y_end; ++y_itr, ++d_itr, ++x_itr ) {
00147 #ifdef LINALGPACK_CHECK_RANGE
00148       TEUCHOS_TEST_FOR_EXCEPT( !(  d_itr < diag.end()  ) );
00149       TEUCHOS_TEST_FOR_EXCEPT( !(  x_itr < vs_rhs2.end()  ) );
00150       TEUCHOS_TEST_FOR_EXCEPT( !(  y_itr < vs_lhs->end()  ) );
00151 #endif
00152       *y_itr = beta * (*y_itr) + alpha * (*d_itr) * (*x_itr);
00153     }
00154   }
00155 }
00156 
00157 void MatrixSymDiagStd::Vp_StMtV(
00158   DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00159   , const SpVectorSlice& sv_rhs2, value_type beta) const
00160 {
00161   const DVectorSlice diag = this->diag();
00162   size_type n = diag.size();
00163 
00164   DenseLinAlgPack::Vp_MtV_assert_sizes( vs_lhs->size()
00165     , n, n, trans_rhs1, sv_rhs2.size() );
00166   //
00167   // y = b*y + a * op(A) * x
00168   //
00169   DenseLinAlgPack::Vt_S(vs_lhs,beta); // y = b * y
00170   //
00171   // A is symmetric and diagonal A = diag(diag) so:
00172   //
00173   // y(j) += a * diag(j) * x(j), for j = 1...n
00174   //
00175   // x is sparse so take account of this.
00176 
00177   for(   SpVectorSlice::const_iterator x_itr = sv_rhs2.begin()
00178      ; x_itr != sv_rhs2.end()
00179      ; ++x_itr )
00180   {
00181     (*vs_lhs)(x_itr->indice() + sv_rhs2.offset())
00182       += alpha * diag(x_itr->indice() + sv_rhs2.offset()) * x_itr->value();
00183       // Note: The indice x(i) invocations are ranged check
00184       // if this is compiled into the code.
00185   }
00186 }
00187 
00188 // Overridden from MatrixWithOpFactorized
00189 
00190 void MatrixSymDiagStd::V_InvMtV(
00191   DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
00192   , const DVectorSlice& vs_rhs2) const
00193 {
00194   const DVectorSlice diag = this->diag();
00195   size_type n = diag.size();
00196 
00197   // y = inv(op(A)) * x
00198   //
00199   // A is symmetric and diagonal (A = diag(diag)) so:
00200   //
00201   // y(j) = x(j) / diag(j), for j = 1...n
00202 
00203   DenseLinAlgPack::Vp_MtV_assert_sizes( vs_lhs->size()
00204     , n, n, trans_rhs1, vs_rhs2.size() );
00205   
00206   if( vs_rhs2.stride() == 1 && vs_lhs->stride() == 1 ) {
00207     // Optimized implementation
00208     const value_type
00209       *d_itr      = diag.raw_ptr(),
00210       *x_itr      = vs_rhs2.raw_ptr();
00211     value_type
00212       *y_itr      = vs_lhs->raw_ptr(),
00213       *y_end      = y_itr + vs_lhs->size();
00214     while( y_itr != y_end )
00215       *y_itr++ = (*x_itr++) / (*d_itr++);
00216   }
00217   else {
00218     // Generic implementation
00219     DVectorSlice::const_iterator
00220       d_itr = diag.begin(),
00221       x_itr = vs_rhs2.begin();
00222     DVectorSlice::iterator
00223       y_itr = vs_lhs->begin(),
00224       y_end = vs_lhs->end();
00225     for( ; y_itr != y_end; ++y_itr, ++d_itr, ++x_itr ) {
00226       TEUCHOS_TEST_FOR_EXCEPT( !(  d_itr < diag.end()  ) );
00227       TEUCHOS_TEST_FOR_EXCEPT( !(  x_itr < vs_rhs2.end()  ) );
00228       TEUCHOS_TEST_FOR_EXCEPT( !(  y_itr < vs_lhs->end()  ) );
00229       *y_itr = (*x_itr)/(*d_itr);
00230     }
00231   }
00232 }
00233 
00234 void MatrixSymDiagStd::V_InvMtV(
00235   DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
00236   , const SpVectorSlice& sv_rhs2) const
00237 {
00238   const DVectorSlice diag = this->diag();
00239   size_type n = diag.size();
00240 
00241   // y = inv(op(A)) * x
00242   //
00243   // A is symmetric and diagonal A = diag(diag) so:
00244   //
00245   // y(j) = x(j) / diag(j), for j = 1...n
00246   //
00247   // x is sparse so take account of this.
00248   
00249   DenseLinAlgPack::Vp_MtV_assert_sizes( vs_lhs->size()
00250     , n, n, trans_rhs1, sv_rhs2.size() );
00251 
00252   for(   SpVectorSlice::const_iterator x_itr = sv_rhs2.begin()
00253      ; x_itr != sv_rhs2.end()
00254      ; ++x_itr )
00255   {
00256     (*vs_lhs)(x_itr->indice() + sv_rhs2.offset())
00257       = x_itr->value() / diag(x_itr->indice() + sv_rhs2.offset());
00258       // Note: The indice x(i) invocations are ranged check
00259       // if this is compiled into the code.
00260   }
00261 }
00262 
00263 } // end namespace AbstractLinAlgPack
00264 
00265 #endif // 0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines