AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects Version of the Day
AbstractLinAlgPack_MatrixSymDiagSparse.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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include <assert.h>
00043 
00044 #include <fstream>    // For debugging only
00045 #include <limits>
00046 
00047 #include "AbstractLinAlgPack_MatrixSymDiagSparse.hpp"
00048 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00049 #include "AbstractLinAlgPack_EtaVector.hpp"
00050 #include "AbstractLinAlgPack_AssertOp.hpp"
00051 #include "AbstractLinAlgPack_SpVectorOut.hpp"
00052 #include "DenseLinAlgPack_DVectorClass.hpp"
00053 #include "DenseLinAlgPack_DMatrixAssign.hpp"
00054 #include "DenseLinAlgPack_DMatrixClass.hpp"
00055 #include "DenseLinAlgPack_DMatrixOp.hpp"
00056 #include "DenseLinAlgPack_assert_print_nan_inf.hpp"
00057 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00058 #include "Teuchos_Assert.hpp"
00059 
00060 namespace {
00061 template< class T >
00062 inline
00063 T my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
00064 } // end namespace
00065 
00066 namespace AbstractLinAlgPack {
00067 
00068 MatrixSymDiagSparse::MatrixSymDiagSparse()
00069   : num_updates_at_once_(0) // Flag that it is to be determined internally.
00070 {}
00071 
00072 // Overridden from MatrixBase
00073 
00074 size_type MatrixSymDiagSparse::rows() const
00075 {
00076   return diag().dim();
00077 }
00078 
00079 // Overridden from MatrixOp
00080 
00081 std::ostream& MatrixSymDiagSparse::output(std::ostream& out) const
00082 {
00083   out << "*** Sparse diagonal matrix ***\n"
00084     << "diag =\n" << diag();
00085   return out;
00086 }
00087 
00088 // Overridden from MatrixOpSerial
00089 
00090 void MatrixSymDiagSparse::Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha
00091   , BLAS_Cpp::Transp trans_rhs1, const DVectorSlice& vs_rhs2, value_type beta) const
00092 {
00093   const SpVectorSlice &diag = this->diag();
00094 
00095   size_type n = diag.dim();
00096 
00097   // Assert that the dimensions of the aruments match up and if not
00098   // then throw an excption.
00099   DenseLinAlgPack::Vp_MtV_assert_sizes( vs_lhs->dim(), n, n, trans_rhs1, vs_rhs2.dim() );
00100 
00101   // y = b*y + a * op(A) * x
00102   //
00103   // A is symmetric and diagonal A = diag(diag) so:
00104   //
00105   // y(j) = b*y(j) + a * diag(j) * x(j), for j = 1...n
00106 
00107   for( SpVectorSlice::const_iterator d_itr = diag.begin(); d_itr != diag.end(); ++d_itr )
00108   {
00109     const size_type i = d_itr->index(); 
00110     (*vs_lhs)(i) = beta * (*vs_lhs)(i) + alpha * d_itr->value() * vs_rhs2(i);
00111   }
00112 }
00113 
00114 // Overridden from MatrixSymOpSerial
00115 
00116 void MatrixSymDiagSparse::Mp_StMtMtM(
00117   DMatrixSliceSym* B, value_type alpha
00118   ,EMatRhsPlaceHolder dummy_place_holder
00119   ,const MatrixOpSerial& A, BLAS_Cpp::Transp A_trans
00120   ,value_type b
00121   ) const
00122 {
00123   using BLAS_Cpp::rows;
00124   using BLAS_Cpp::cols;
00125   using BLAS_Cpp::trans_not;
00126 
00127   using DenseLinAlgPack::nonconst_tri_ele;
00128   using DenseLinAlgPack::assign;
00129   using DenseLinAlgPack::syrk;
00130   using DenseLinAlgPack::assert_print_nan_inf;
00131 
00132   using LinAlgOpPack::V_MtV;
00133 
00134   typedef EtaVector eta;
00135 
00136   // Assert the size matching of M * op(A)
00137   DenseLinAlgPack::MtV_assert_sizes(
00138       this->rows(), this->cols(), BLAS_Cpp::no_trans
00139     , rows( A.rows(), A.cols(), A_trans ) );
00140 
00141   // Assert size matchin of B = (op(A') * M * op(A))
00142   DenseLinAlgPack::Vp_V_assert_sizes(
00143       B->cols(), cols( A.rows(), A.cols(), A_trans ) );
00144 
00145   //
00146   // B = a * op(A') * M * op(A)
00147   //
00148   //   = a * op(A') * M^(1/2) * M^(1/2) * op(A)
00149   //
00150   //   = a * E * E'
00151   //
00152   // E = M^(1/2) * op(A)
00153   //
00154   //     [ .                                                 ] [ .              ]
00155   //     [   sqrt(M(j(1)))                                   ] [ op(A)(j(1),:)  ]
00156   //     [                .                                  ] [ .              ]
00157   //   = [                  sqrt(M(j(i))                     ] [ op(A)(j(i),:)  ]
00158   //     [                              .                    ] [ .              ]
00159   //     [                                sqrt(M(j(nz))      ] [ op(A)(j(nz),:) ]
00160   //     [                                               .   ] [ .              ]
00161   //
00162   //
00163   //     [ .        ]
00164   //     [ d(j(1))' ]
00165   //     [ .        ]
00166   //   = [ d(j(i))' ]
00167   //     [ .        ]
00168   //     [ d(j(1))' ]
00169   //     [ .        ]
00170   //
00171   //     where: d(j(i)) = sqrt(M(j(i)) * op(A')(:,j(i))    <: R^m 
00172   //                    = sqrt(M(j(i)) * op(A') * e(j(i))  <: R^m
00173   //
00174   //  Above M^(1/2) only has nz nonzero elements sqrt(M(j(i)), i = 1..nz and only
00175   //  the corresponding rows of op(A)(j(i),:), i = 1..nz are shown.  A may in fact
00176   //  dense matrix but the columns are extracted through op(A)*eta(j(i)), i=1..nz.
00177   //
00178   //  The above product B = a * E * E' is a set of nz rank-1 updates and can be written
00179   //  in the form:
00180   //
00181   //  B = sum( a * d(j(i)) * d(j(i))', i = 1..nz )
00182   //
00183   //  Since it is more efficient to perform several rank-1 updates at a time we will
00184   //  perform them in blocks.
00185   //
00186   //  B = B + D(k) * D(k)', k = 1 .. num_blocks
00187   //
00188   //      where:
00189   //         num_blocks = nz / num_updates_at_once + 1 (integer division)
00190   //         D(k) = [ d(j(i1)) ... d(j(i2)) ]
00191   //         i1 = (k-1) * num_updates_at_once + 1
00192   //         i2 = i1 + num_updates_at_once - 1
00193   
00194   const SpVectorSlice
00195     &diag = this->diag();
00196 
00197   const size_type
00198     n = this->rows(),
00199     m = cols(A.rows(),A.cols(),A_trans);
00200 
00201   // Get the actual number of updates to use per rank-(num_updates) update
00202   const size_type
00203     num_updates
00204       = my_min( num_updates_at_once()
00205               ? num_updates_at_once()
00206               : 20  // There may be a better default value for this?
00207             , diag.nz()
00208             );
00209 
00210   // Get the number of blocks of rank-(num_updates) updates
00211   size_type
00212     num_blocks = diag.nz() / num_updates;
00213   if( diag.nz() % num_updates > 0 )
00214     num_blocks++;
00215 
00216   // Initialize B = b*B
00217   if( b != 1.0 )
00218     assign( &nonconst_tri_ele( B->gms(), B->uplo() ), 0.0 );
00219 
00220   // Perform the rank-(num_updates) updates
00221   DMatrix D(m,num_updates);
00222   for( size_type k = 1; k <= num_blocks; ++k ) {
00223     const size_type
00224       i1 = (k-1) * num_updates + 1,
00225       i2 = my_min( diag.nz(), i1 + num_updates - 1 );
00226     // Generate the colunns of D(k)
00227     SpVectorSlice::const_iterator
00228       m_itr = diag.begin() + (i1-1);
00229     for( size_type l = 1; l <= i2-i1+1; ++l, ++m_itr ) {
00230       TEUCHOS_TEST_FOR_EXCEPT( !(  m_itr < diag.end()  ) );
00231       TEUCHOS_TEST_FOR_EXCEPT( !(  m_itr->value() >= 0.0  ) );
00232       V_MtV( &D.col(l), A, trans_not(A_trans)
00233         , eta( m_itr->index(), n, std::sqrt(m_itr->value()) )() );
00234     }
00235     const DMatrixSlice
00236       D_update = D(1,m,1,i2-i1+1);
00237 
00238 
00239 //    // For debugging only
00240 //    std::ofstream ofile("MatrixSymDiagonalSparse_Error.out");
00241 //    assert_print_nan_inf( D_update, "D", true, &ofile );
00242     // Perform the rank-(num_updates) update
00243     syrk( BLAS_Cpp::no_trans, alpha, D_update, 1.0, B );
00244   }
00245 }
00246 
00247 // Overridden from MatrixConvertToSparse
00248 
00249 index_type
00250 MatrixSymDiagSparse::num_nonzeros(
00251   EExtractRegion        extract_region
00252   ,EElementUniqueness   element_uniqueness
00253   ) const
00254 {
00255   return diag().nz();
00256 }
00257 
00258 void MatrixSymDiagSparse::coor_extract_nonzeros(
00259   EExtractRegion                extract_region
00260   ,EElementUniqueness           element_uniqueness
00261   ,const index_type             len_Aval
00262   ,value_type                   Aval[]
00263   ,const index_type             len_Aij
00264   ,index_type                   Arow[]
00265   ,index_type                   Acol[]
00266   ,const index_type             row_offset
00267   ,const index_type             col_offset
00268   ) const
00269 {
00270   const SpVectorSlice
00271     &diag = this->diag();
00272 
00273   TEUCHOS_TEST_FOR_EXCEPTION(
00274     (len_Aval != 0 ? len_Aval != diag.nz() : Aval != NULL)
00275     ,std::invalid_argument
00276     ,"MatrixSymDiagSparse::coor_extract_nonzeros(...): Error!" );
00277   TEUCHOS_TEST_FOR_EXCEPTION(
00278     (len_Aij != 0 ? len_Aij != diag.nz() : (Acol != NULL || Acol != NULL) )
00279     ,std::invalid_argument
00280     ,"MatrixSymDiagSparse::coor_extract_nonzeros(...): Error!" );
00281 
00282   if( len_Aval > 0 ) {
00283     SpVectorSlice::const_iterator
00284       itr;
00285     FortranTypes::f_dbl_prec
00286       *l_Aval;
00287     for( itr = diag.begin(), l_Aval = Aval; itr != diag.end(); ++itr ) {
00288       *l_Aval++ = itr->value();
00289     }     
00290   }
00291   if( len_Aij > 0 ) {
00292     SpVectorSlice::const_iterator
00293       itr;
00294     index_type
00295       *l_Arow, *l_Acol;
00296     for( itr = diag.begin(), l_Arow = Arow, l_Acol = Acol; itr != diag.end(); ++itr ) {
00297       const index_type
00298         ij = itr->index() + diag.offset();
00299       *l_Arow++ = ij + row_offset;
00300       *l_Acol++ = ij + col_offset;
00301     }     
00302   }
00303 }
00304 
00305 } // end namespace AbstractLinAlgPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends