AbstractLinAlgPack_MatrixSymDiagonalStd.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 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixSymDiagStd.hpp"
00030 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_SpVectorClass.hpp"
00031 #include "DenseLinAlgPack_DVectorOp.hpp"
00032 #include "DenseLinAlgPack_DMatrixOp.hpp"
00033 #include "DenseLinAlgPack_AssertOp.hpp"
00034 
00035 namespace AbstractLinAlgPack {
00036 
00037 MatrixSymDiagStd::MatrixSymDiagStd()
00038 {}
00039 
00040 DVectorSlice MatrixSymDiagStd::diag()
00041 {
00042   return diag_();
00043 }
00044 
00045 const DVectorSlice MatrixSymDiagStd::diag() const
00046 {
00047   return diag_();
00048 }
00049 
00050 // Overridden from MatrixSymInitDiag
00051 
00052 void MatrixSymDiagStd::init_identity( size_type n, value_type alpha = 1.0 )
00053 {
00054   diag_.resize(n);
00055   if(n)
00056     diag_ = alpha;
00057 }
00058 
00059 void MatrixSymDiagStd::init_diagonal( const DVectorSlice& diag )
00060 {
00061   diag_ = diag;
00062 }
00063 
00064 // Overridden from Matrix
00065 
00066 size_type MatrixSymDiagStd::rows() const
00067 {
00068   return diag_.size();
00069 }
00070 
00071 // Overridden from MatrixOp
00072 
00073 void MatrixSymDiagStd::Mp_StM(
00074   DMatrixSlice* gms_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs) const
00075 {
00076   using DenseLinAlgPack::Vp_StV;
00077   // Assert that the dimensions match up.
00078   DenseLinAlgPack::Mp_M_assert_sizes( gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
00079     , rows(), cols(), trans_rhs );
00080   // Add to the diagonal
00081   Vp_StV( &gms_lhs->diag(), alpha, diag_ );
00082 }
00083 
00084 void MatrixSymDiagStd::Vp_StMtV(
00085   DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00086   , const DVectorSlice& vs_rhs2, value_type beta) const
00087 {
00088   const DVectorSlice diag = this->diag();
00089   size_type n = diag.size();
00090 
00091   //
00092   // y = b*y + a * op(A) * x
00093   //
00094   DenseLinAlgPack::Vp_MtV_assert_sizes(
00095     vs_lhs->size(), n, n, trans_rhs1, vs_rhs2.size() );
00096   //
00097   // A is symmetric and diagonal A = diag(diag) so:
00098   //
00099   // y(j) += a * diag(j) * x(j), for j = 1...n
00100   //
00101   if( vs_rhs2.stride() == 1 && vs_lhs->stride() == 1 ) {
00102     // Optimized implementation
00103     const value_type
00104       *d_itr      = diag.raw_ptr(),
00105       *x_itr      = vs_rhs2.raw_ptr();
00106     value_type
00107       *y_itr      = vs_lhs->raw_ptr(),
00108       *y_end      = y_itr + vs_lhs->size();
00109 
00110     if( beta == 0.0 ) {
00111       while( y_itr != y_end )
00112         *y_itr++ = alpha * (*d_itr++) * (*x_itr++);
00113     }
00114     else if( beta == 1.0 ) {
00115       while( y_itr != y_end )
00116         *y_itr++ += alpha * (*d_itr++) * (*x_itr++);
00117     }
00118     else {
00119       for( ; y_itr != y_end; ++y_itr )
00120         *y_itr = beta * (*y_itr) + alpha * (*d_itr++) * (*x_itr++);
00121     }
00122   }
00123   else {
00124     // Generic implementation
00125     DVectorSlice::const_iterator
00126       d_itr = diag.begin(),
00127       x_itr = vs_rhs2.begin();
00128     DVectorSlice::iterator
00129       y_itr = vs_lhs->begin(),
00130       y_end = vs_lhs->end();
00131     for( ; y_itr != y_end; ++y_itr, ++d_itr, ++x_itr ) {
00132 #ifdef LINALGPACK_CHECK_RANGE
00133       assert( d_itr < diag.end() );
00134       assert( x_itr < vs_rhs2.end() );
00135       assert( y_itr < vs_lhs->end() );
00136 #endif
00137       *y_itr = beta * (*y_itr) + alpha * (*d_itr) * (*x_itr);
00138     }
00139   }
00140 }
00141 
00142 void MatrixSymDiagStd::Vp_StMtV(
00143   DVectorSlice* vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00144   , const SpVectorSlice& sv_rhs2, value_type beta) const
00145 {
00146   const DVectorSlice diag = this->diag();
00147   size_type n = diag.size();
00148 
00149   DenseLinAlgPack::Vp_MtV_assert_sizes( vs_lhs->size()
00150     , n, n, trans_rhs1, sv_rhs2.size() );
00151   //
00152   // y = b*y + a * op(A) * x
00153   //
00154   DenseLinAlgPack::Vt_S(vs_lhs,beta); // y = b * y
00155   //
00156   // A is symmetric and diagonal A = diag(diag) so:
00157   //
00158   // y(j) += a * diag(j) * x(j), for j = 1...n
00159   //
00160   // x is sparse so take account of this.
00161 
00162   for(   SpVectorSlice::const_iterator x_itr = sv_rhs2.begin()
00163      ; x_itr != sv_rhs2.end()
00164      ; ++x_itr )
00165   {
00166     (*vs_lhs)(x_itr->indice() + sv_rhs2.offset())
00167       += alpha * diag(x_itr->indice() + sv_rhs2.offset()) * x_itr->value();
00168       // Note: The indice x(i) invocations are ranged check
00169       // if this is compiled into the code.
00170   }
00171 }
00172 
00173 // Overridden from MatrixWithOpFactorized
00174 
00175 void MatrixSymDiagStd::V_InvMtV(
00176   DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
00177   , const DVectorSlice& vs_rhs2) const
00178 {
00179   const DVectorSlice diag = this->diag();
00180   size_type n = diag.size();
00181 
00182   // y = inv(op(A)) * x
00183   //
00184   // A is symmetric and diagonal (A = diag(diag)) so:
00185   //
00186   // y(j) = x(j) / diag(j), for j = 1...n
00187 
00188   DenseLinAlgPack::Vp_MtV_assert_sizes( vs_lhs->size()
00189     , n, n, trans_rhs1, vs_rhs2.size() );
00190   
00191   if( vs_rhs2.stride() == 1 && vs_lhs->stride() == 1 ) {
00192     // Optimized implementation
00193     const value_type
00194       *d_itr      = diag.raw_ptr(),
00195       *x_itr      = vs_rhs2.raw_ptr();
00196     value_type
00197       *y_itr      = vs_lhs->raw_ptr(),
00198       *y_end      = y_itr + vs_lhs->size();
00199     while( y_itr != y_end )
00200       *y_itr++ = (*x_itr++) / (*d_itr++);
00201   }
00202   else {
00203     // Generic implementation
00204     DVectorSlice::const_iterator
00205       d_itr = diag.begin(),
00206       x_itr = vs_rhs2.begin();
00207     DVectorSlice::iterator
00208       y_itr = vs_lhs->begin(),
00209       y_end = vs_lhs->end();
00210     for( ; y_itr != y_end; ++y_itr, ++d_itr, ++x_itr ) {
00211       assert( d_itr < diag.end() );
00212       assert( x_itr < vs_rhs2.end() );
00213       assert( y_itr < vs_lhs->end() );
00214       *y_itr = (*x_itr)/(*d_itr);
00215     }
00216   }
00217 }
00218 
00219 void MatrixSymDiagStd::V_InvMtV(
00220   DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1
00221   , const SpVectorSlice& sv_rhs2) const
00222 {
00223   const DVectorSlice diag = this->diag();
00224   size_type n = diag.size();
00225 
00226   // y = inv(op(A)) * x
00227   //
00228   // A is symmetric and diagonal A = diag(diag) so:
00229   //
00230   // y(j) = x(j) / diag(j), for j = 1...n
00231   //
00232   // x is sparse so take account of this.
00233   
00234   DenseLinAlgPack::Vp_MtV_assert_sizes( vs_lhs->size()
00235     , n, n, trans_rhs1, sv_rhs2.size() );
00236 
00237   for(   SpVectorSlice::const_iterator x_itr = sv_rhs2.begin()
00238      ; x_itr != sv_rhs2.end()
00239      ; ++x_itr )
00240   {
00241     (*vs_lhs)(x_itr->indice() + sv_rhs2.offset())
00242       = x_itr->value() / diag(x_itr->indice() + sv_rhs2.offset());
00243       // Note: The indice x(i) invocations are ranged check
00244       // if this is compiled into the code.
00245   }
00246 }
00247 
00248 } // end namespace AbstractLinAlgPack

Generated on Thu Sep 18 12:35:12 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1