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 // This library is free software; you can redistribute it and/or modify
00013 // it under the terms of the GNU Lesser General Public License as
00014 // published by the Free Software Foundation; either version 2.1 of the
00015 // License, or (at your option) any later version.
00016 //  
00017 // This library is distributed in the hope that it will be useful, but
00018 // WITHOUT ANY WARRANTY; without even the implied warranty of
00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 // Lesser General Public License for more details.
00021 //  
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License along with this library; if not, write to the Free Software
00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 // USA
00026 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00027 // 
00028 // ***********************************************************************
00029 // @HEADER
00030 
00031 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixSymDiagStd.hpp"
00032 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_SpVectorClass.hpp"
00033 #include "DenseLinAlgPack_DVectorOp.hpp"
00034 #include "DenseLinAlgPack_DMatrixOp.hpp"
00035 #include "DenseLinAlgPack_AssertOp.hpp"
00036 
00037 namespace AbstractLinAlgPack {
00038 
00039 MatrixSymDiagStd::MatrixSymDiagStd()
00040 {}
00041 
00042 DVectorSlice MatrixSymDiagStd::diag()
00043 {
00044   return diag_();
00045 }
00046 
00047 const DVectorSlice MatrixSymDiagStd::diag() const
00048 {
00049   return diag_();
00050 }
00051 
00052 // Overridden from MatrixSymInitDiag
00053 
00054 void MatrixSymDiagStd::init_identity( size_type n, value_type alpha = 1.0 )
00055 {
00056   diag_.resize(n);
00057   if(n)
00058     diag_ = alpha;
00059 }
00060 
00061 void MatrixSymDiagStd::init_diagonal( const DVectorSlice& diag )
00062 {
00063   diag_ = diag;
00064 }
00065 
00066 // Overridden from Matrix
00067 
00068 size_type MatrixSymDiagStd::rows() const
00069 {
00070   return diag_.size();
00071 }
00072 
00073 // Overridden from MatrixOp
00074 
00075 void MatrixSymDiagStd::Mp_StM(
00076   DMatrixSlice* gms_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs) const
00077 {
00078   using DenseLinAlgPack::Vp_StV;
00079   // Assert that the dimensions match up.
00080   DenseLinAlgPack::Mp_M_assert_sizes( gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
00081     , rows(), cols(), trans_rhs );
00082   // Add to the diagonal
00083   Vp_StV( &gms_lhs->diag(), alpha, diag_ );
00084 }
00085 
00086 void MatrixSymDiagStd::Vp_StMtV(
00087   DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00088   , const DVectorSlice& vs_rhs2, value_type beta) const
00089 {
00090   const DVectorSlice diag = this->diag();
00091   size_type n = diag.size();
00092 
00093   //
00094   // y = b*y + a * op(A) * x
00095   //
00096   DenseLinAlgPack::Vp_MtV_assert_sizes(
00097     vs_lhs->size(), n, n, trans_rhs1, vs_rhs2.size() );
00098   //
00099   // A is symmetric and diagonal A = diag(diag) so:
00100   //
00101   // y(j) += a * diag(j) * x(j), for j = 1...n
00102   //
00103   if( vs_rhs2.stride() == 1 && vs_lhs->stride() == 1 ) {
00104     // Optimized implementation
00105     const value_type
00106       *d_itr      = diag.raw_ptr(),
00107       *x_itr      = vs_rhs2.raw_ptr();
00108     value_type
00109       *y_itr      = vs_lhs->raw_ptr(),
00110       *y_end      = y_itr + vs_lhs->size();
00111 
00112     if( beta == 0.0 ) {
00113       while( y_itr != y_end )
00114         *y_itr++ = alpha * (*d_itr++) * (*x_itr++);
00115     }
00116     else if( beta == 1.0 ) {
00117       while( y_itr != y_end )
00118         *y_itr++ += alpha * (*d_itr++) * (*x_itr++);
00119     }
00120     else {
00121       for( ; y_itr != y_end; ++y_itr )
00122         *y_itr = beta * (*y_itr) + alpha * (*d_itr++) * (*x_itr++);
00123     }
00124   }
00125   else {
00126     // Generic implementation
00127     DVectorSlice::const_iterator
00128       d_itr = diag.begin(),
00129       x_itr = vs_rhs2.begin();
00130     DVectorSlice::iterator
00131       y_itr = vs_lhs->begin(),
00132       y_end = vs_lhs->end();
00133     for( ; y_itr != y_end; ++y_itr, ++d_itr, ++x_itr ) {
00134 #ifdef LINALGPACK_CHECK_RANGE
00135       TEST_FOR_EXCEPT( !(  d_itr < diag.end()  ) );
00136       TEST_FOR_EXCEPT( !(  x_itr < vs_rhs2.end()  ) );
00137       TEST_FOR_EXCEPT( !(  y_itr < vs_lhs->end()  ) );
00138 #endif
00139       *y_itr = beta * (*y_itr) + alpha * (*d_itr) * (*x_itr);
00140     }
00141   }
00142 }
00143 
00144 void MatrixSymDiagStd::Vp_StMtV(
00145   DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00146   , const SpVectorSlice& sv_rhs2, value_type beta) const
00147 {
00148   const DVectorSlice diag = this->diag();
00149   size_type n = diag.size();
00150 
00151   DenseLinAlgPack::Vp_MtV_assert_sizes( vs_lhs->size()
00152     , n, n, trans_rhs1, sv_rhs2.size() );
00153   //
00154   // y = b*y + a * op(A) * x
00155   //
00156   DenseLinAlgPack::Vt_S(vs_lhs,beta); // y = b * y
00157   //
00158   // A is symmetric and diagonal A = diag(diag) so:
00159   //
00160   // y(j) += a * diag(j) * x(j), for j = 1...n
00161   //
00162   // x is sparse so take account of this.
00163 
00164   for(   SpVectorSlice::const_iterator x_itr = sv_rhs2.begin()
00165      ; x_itr != sv_rhs2.end()
00166      ; ++x_itr )
00167   {
00168     (*vs_lhs)(x_itr->indice() + sv_rhs2.offset())
00169       += alpha * diag(x_itr->indice() + sv_rhs2.offset()) * x_itr->value();
00170       // Note: The indice x(i) invocations are ranged check
00171       // if this is compiled into the code.
00172   }
00173 }
00174 
00175 // Overridden from MatrixWithOpFactorized
00176 
00177 void MatrixSymDiagStd::V_InvMtV(
00178   DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
00179   , const DVectorSlice& vs_rhs2) const
00180 {
00181   const DVectorSlice diag = this->diag();
00182   size_type n = diag.size();
00183 
00184   // y = inv(op(A)) * x
00185   //
00186   // A is symmetric and diagonal (A = diag(diag)) so:
00187   //
00188   // y(j) = x(j) / diag(j), for j = 1...n
00189 
00190   DenseLinAlgPack::Vp_MtV_assert_sizes( vs_lhs->size()
00191     , n, n, trans_rhs1, vs_rhs2.size() );
00192   
00193   if( vs_rhs2.stride() == 1 && vs_lhs->stride() == 1 ) {
00194     // Optimized implementation
00195     const value_type
00196       *d_itr      = diag.raw_ptr(),
00197       *x_itr      = vs_rhs2.raw_ptr();
00198     value_type
00199       *y_itr      = vs_lhs->raw_ptr(),
00200       *y_end      = y_itr + vs_lhs->size();
00201     while( y_itr != y_end )
00202       *y_itr++ = (*x_itr++) / (*d_itr++);
00203   }
00204   else {
00205     // Generic implementation
00206     DVectorSlice::const_iterator
00207       d_itr = diag.begin(),
00208       x_itr = vs_rhs2.begin();
00209     DVectorSlice::iterator
00210       y_itr = vs_lhs->begin(),
00211       y_end = vs_lhs->end();
00212     for( ; y_itr != y_end; ++y_itr, ++d_itr, ++x_itr ) {
00213       TEST_FOR_EXCEPT( !(  d_itr < diag.end()  ) );
00214       TEST_FOR_EXCEPT( !(  x_itr < vs_rhs2.end()  ) );
00215       TEST_FOR_EXCEPT( !(  y_itr < vs_lhs->end()  ) );
00216       *y_itr = (*x_itr)/(*d_itr);
00217     }
00218   }
00219 }
00220 
00221 void MatrixSymDiagStd::V_InvMtV(
00222   DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
00223   , const SpVectorSlice& sv_rhs2) const
00224 {
00225   const DVectorSlice diag = this->diag();
00226   size_type n = diag.size();
00227 
00228   // y = inv(op(A)) * x
00229   //
00230   // A is symmetric and diagonal A = diag(diag) so:
00231   //
00232   // y(j) = x(j) / diag(j), for j = 1...n
00233   //
00234   // x is sparse so take account of this.
00235   
00236   DenseLinAlgPack::Vp_MtV_assert_sizes( vs_lhs->size()
00237     , n, n, trans_rhs1, sv_rhs2.size() );
00238 
00239   for(   SpVectorSlice::const_iterator x_itr = sv_rhs2.begin()
00240      ; x_itr != sv_rhs2.end()
00241      ; ++x_itr )
00242   {
00243     (*vs_lhs)(x_itr->indice() + sv_rhs2.offset())
00244       = x_itr->value() / diag(x_itr->indice() + sv_rhs2.offset());
00245       // Note: The indice x(i) invocations are ranged check
00246       // if this is compiled into the code.
00247   }
00248 }
00249 
00250 } // end namespace AbstractLinAlgPack
00251 
00252 #endif // 0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:10:56 2011 for MOOCHO (Single Doxygen Collection) by  doxygen 1.6.3