AbstractLinAlgPack_MatrixSymDiagSparse.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 <assert.h>
00030 
00031 #include <fstream>    // For debugging only
00032 #include <limits>
00033 
00034 #include "AbstractLinAlgPack_MatrixSymDiagSparse.hpp"
00035 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00036 #include "AbstractLinAlgPack_EtaVector.hpp"
00037 #include "AbstractLinAlgPack_AssertOp.hpp"
00038 #include "AbstractLinAlgPack_SpVectorOut.hpp"
00039 #include "DenseLinAlgPack_DVectorClass.hpp"
00040 #include "DenseLinAlgPack_DMatrixAssign.hpp"
00041 #include "DenseLinAlgPack_DMatrixClass.hpp"
00042 #include "DenseLinAlgPack_DMatrixOp.hpp"
00043 #include "DenseLinAlgPack_assert_print_nan_inf.hpp"
00044 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00045 #include "Teuchos_TestForException.hpp"
00046 
00047 namespace {
00048 template< class T >
00049 inline
00050 T my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
00051 } // end namespace
00052 
00053 namespace AbstractLinAlgPack {
00054 
00055 MatrixSymDiagSparse::MatrixSymDiagSparse()
00056   : num_updates_at_once_(0) // Flag that it is to be determined internally.
00057 {}
00058 
00059 // Overridden from MatrixBase
00060 
00061 size_type MatrixSymDiagSparse::rows() const
00062 {
00063   return diag().dim();
00064 }
00065 
00066 // Overridden from MatrixOp
00067 
00068 std::ostream& MatrixSymDiagSparse::output(std::ostream& out) const
00069 {
00070   out << "*** Sparse diagonal matrix ***\n"
00071     << "diag =\n" << diag();
00072   return out;
00073 }
00074 
00075 // Overridden from MatrixOpSerial
00076 
00077 void MatrixSymDiagSparse::Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha
00078   , BLAS_Cpp::Transp trans_rhs1, const DVectorSlice& vs_rhs2, value_type beta) const
00079 {
00080   const SpVectorSlice &diag = this->diag();
00081 
00082   size_type n = diag.dim();
00083 
00084   // Assert that the dimensions of the aruments match up and if not
00085   // then throw an excption.
00086   DenseLinAlgPack::Vp_MtV_assert_sizes( vs_lhs->dim(), n, n, trans_rhs1, vs_rhs2.dim() );
00087 
00088   // y = b*y + a * op(A) * x
00089   //
00090   // A is symmetric and diagonal A = diag(diag) so:
00091   //
00092   // y(j) = b*y(j) + a * diag(j) * x(j), for j = 1...n
00093 
00094   for( SpVectorSlice::const_iterator d_itr = diag.begin(); d_itr != diag.end(); ++d_itr )
00095   {
00096     const size_type i = d_itr->index(); 
00097     (*vs_lhs)(i) = beta * (*vs_lhs)(i) + alpha * d_itr->value() * vs_rhs2(i);
00098   }
00099 }
00100 
00101 // Overridden from MatrixSymOpSerial
00102 
00103 void MatrixSymDiagSparse::Mp_StMtMtM(
00104   DMatrixSliceSym* B, value_type alpha
00105   ,EMatRhsPlaceHolder dummy_place_holder
00106   ,const MatrixOpSerial& A, BLAS_Cpp::Transp A_trans
00107   ,value_type b
00108   ) const
00109 {
00110   using BLAS_Cpp::rows;
00111   using BLAS_Cpp::cols;
00112   using BLAS_Cpp::trans_not;
00113 
00114   using DenseLinAlgPack::nonconst_tri_ele;
00115   using DenseLinAlgPack::assign;
00116   using DenseLinAlgPack::syrk;
00117   using DenseLinAlgPack::assert_print_nan_inf;
00118 
00119   using LinAlgOpPack::V_MtV;
00120 
00121   typedef EtaVector eta;
00122 
00123   // Assert the size matching of M * op(A)
00124   DenseLinAlgPack::MtV_assert_sizes(
00125       this->rows(), this->cols(), BLAS_Cpp::no_trans
00126     , rows( A.rows(), A.cols(), A_trans ) );
00127 
00128   // Assert size matchin of B = (op(A') * M * op(A))
00129   DenseLinAlgPack::Vp_V_assert_sizes(
00130       B->cols(), cols( A.rows(), A.cols(), A_trans ) );
00131 
00132   //
00133   // B = a * op(A') * M * op(A)
00134   //
00135   //   = a * op(A') * M^(1/2) * M^(1/2) * op(A)
00136   //
00137   //   = a * E * E'
00138   //
00139   // E = M^(1/2) * op(A)
00140   //
00141   //     [ .                                                 ] [ .              ]
00142   //     [   sqrt(M(j(1)))                                   ] [ op(A)(j(1),:)  ]
00143   //     [                .                                  ] [ .              ]
00144   //   = [                  sqrt(M(j(i))                     ] [ op(A)(j(i),:)  ]
00145   //     [                              .                    ] [ .              ]
00146   //     [                                sqrt(M(j(nz))      ] [ op(A)(j(nz),:) ]
00147   //     [                                               .   ] [ .              ]
00148   //
00149   //
00150   //     [ .        ]
00151   //     [ d(j(1))' ]
00152   //     [ .        ]
00153   //   = [ d(j(i))' ]
00154   //     [ .        ]
00155   //     [ d(j(1))' ]
00156   //     [ .        ]
00157   //
00158   //     where: d(j(i)) = sqrt(M(j(i)) * op(A')(:,j(i))    <: R^m 
00159   //                    = sqrt(M(j(i)) * op(A') * e(j(i))  <: R^m
00160   //
00161   //  Above M^(1/2) only has nz nonzero elements sqrt(M(j(i)), i = 1..nz and only
00162   //  the corresponding rows of op(A)(j(i),:), i = 1..nz are shown.  A may in fact
00163   //  dense matrix but the columns are extracted through op(A)*eta(j(i)), i=1..nz.
00164   //
00165   //  The above product B = a * E * E' is a set of nz rank-1 updates and can be written
00166   //  in the form:
00167   //
00168   //  B = sum( a * d(j(i)) * d(j(i))', i = 1..nz )
00169   //
00170   //  Since it is more efficient to perform several rank-1 updates at a time we will
00171   //  perform them in blocks.
00172   //
00173   //  B = B + D(k) * D(k)', k = 1 .. num_blocks
00174   //
00175   //      where:
00176   //         num_blocks = nz / num_updates_at_once + 1 (integer division)
00177   //         D(k) = [ d(j(i1)) ... d(j(i2)) ]
00178   //         i1 = (k-1) * num_updates_at_once + 1
00179   //         i2 = i1 + num_updates_at_once - 1
00180   
00181   const SpVectorSlice
00182     &diag = this->diag();
00183 
00184   const size_type
00185     n = this->rows(),
00186     m = cols(A.rows(),A.cols(),A_trans);
00187 
00188   // Get the actual number of updates to use per rank-(num_updates) update
00189   const size_type
00190     num_updates
00191       = my_min( num_updates_at_once()
00192               ? num_updates_at_once()
00193               : 20  // There may be a better default value for this?
00194             , diag.nz()
00195             );
00196 
00197   // Get the number of blocks of rank-(num_updates) updates
00198   size_type
00199     num_blocks = diag.nz() / num_updates;
00200   if( diag.nz() % num_updates > 0 )
00201     num_blocks++;
00202 
00203   // Initialize B = b*B
00204   if( b != 1.0 )
00205     assign( &nonconst_tri_ele( B->gms(), B->uplo() ), 0.0 );
00206 
00207   // Perform the rank-(num_updates) updates
00208   DMatrix D(m,num_updates);
00209   for( size_type k = 1; k <= num_blocks; ++k ) {
00210     const size_type
00211       i1 = (k-1) * num_updates + 1,
00212       i2 = my_min( diag.nz(), i1 + num_updates - 1 );
00213     // Generate the colunns of D(k)
00214     SpVectorSlice::const_iterator
00215       m_itr = diag.begin() + (i1-1);
00216     for( size_type l = 1; l <= i2-i1+1; ++l, ++m_itr ) {
00217       assert( m_itr < diag.end() );
00218       assert( m_itr->value() >= 0.0 );
00219       V_MtV( &D.col(l), A, trans_not(A_trans)
00220         , eta( m_itr->index(), n, ::sqrt(m_itr->value()) )() );
00221     }
00222     const DMatrixSlice
00223       D_update = D(1,m,1,i2-i1+1);
00224 
00225 
00226 //    // For debugging only
00227 //    std::ofstream ofile("MatrixSymDiagonalSparse_Error.out");
00228 //    assert_print_nan_inf( D_update, "D", true, &ofile );
00229     // Perform the rank-(num_updates) update
00230     syrk( BLAS_Cpp::no_trans, alpha, D_update, 1.0, B );
00231   }
00232 }
00233 
00234 // Overridden from MatrixConvertToSparse
00235 
00236 index_type
00237 MatrixSymDiagSparse::num_nonzeros(
00238   EExtractRegion        extract_region
00239   ,EElementUniqueness   element_uniqueness
00240   ) const
00241 {
00242   return diag().nz();
00243 }
00244 
00245 void MatrixSymDiagSparse::coor_extract_nonzeros(
00246   EExtractRegion                extract_region
00247   ,EElementUniqueness           element_uniqueness
00248   ,const index_type             len_Aval
00249   ,value_type                   Aval[]
00250   ,const index_type             len_Aij
00251   ,index_type                   Arow[]
00252   ,index_type                   Acol[]
00253   ,const index_type             row_offset
00254   ,const index_type             col_offset
00255   ) const
00256 {
00257   const SpVectorSlice
00258     &diag = this->diag();
00259 
00260   TEST_FOR_EXCEPTION(
00261     (len_Aval != 0 ? len_Aval != diag.nz() : Aval != NULL)
00262     ,std::invalid_argument
00263     ,"MatrixSymDiagSparse::coor_extract_nonzeros(...): Error!" );
00264   TEST_FOR_EXCEPTION(
00265     (len_Aij != 0 ? len_Aij != diag.nz() : (Acol != NULL || Acol != NULL) )
00266     ,std::invalid_argument
00267     ,"MatrixSymDiagSparse::coor_extract_nonzeros(...): Error!" );
00268 
00269   if( len_Aval > 0 ) {
00270     SpVectorSlice::const_iterator
00271       itr;
00272     FortranTypes::f_dbl_prec
00273       *l_Aval;
00274     for( itr = diag.begin(), l_Aval = Aval; itr != diag.end(); ++itr ) {
00275       *l_Aval++ = itr->value();
00276     }     
00277   }
00278   if( len_Aij > 0 ) {
00279     SpVectorSlice::const_iterator
00280       itr;
00281     FortranTypes::f_int
00282       *l_Arow, *l_Acol;
00283     for( itr = diag.begin(), l_Arow = Arow, l_Acol = Acol; itr != diag.end(); ++itr ) {
00284       const FortranTypes::f_int
00285         ij = itr->index() + diag.offset();
00286       *l_Arow++ = ij + row_offset;
00287       *l_Acol++ = ij + col_offset;
00288     }     
00289   }
00290 }
00291 
00292 } // end namespace AbstractLinAlgPack

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