AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects Version of the Day
AbstractLinAlgPack_MatrixOpNonsing.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 <math.h>
00043 
00044 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
00045 #include "AbstractLinAlgPack_VectorSpace.hpp"
00046 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00047 #include "Teuchos_Assert.hpp"
00048 
00049 namespace AbstractLinAlgPack {
00050 
00051 MatrixOpNonsing::mat_mwons_mut_ptr_t
00052 MatrixOpNonsing::clone_mwons()
00053 {
00054   return Teuchos::null;
00055 }
00056 
00057 MatrixOpNonsing::mat_mwons_ptr_t
00058 MatrixOpNonsing::clone_mwons() const
00059 {
00060   return Teuchos::null;
00061 }
00062 
00063 const MatrixOp::MatNorm
00064 MatrixOpNonsing::calc_cond_num(
00065   EMatNormType  requested_norm_type
00066   ,bool         allow_replacement
00067   ) const
00068 {
00069   using BLAS_Cpp::no_trans;
00070   using BLAS_Cpp::trans;
00071   using LinAlgOpPack::V_InvMtV;
00072   const VectorSpace
00073     &space_cols = this->space_cols(),
00074     &space_rows = this->space_rows();
00075   const index_type
00076     num_cols = space_rows.dim();
00077   TEUCHOS_TEST_FOR_EXCEPTION(
00078     !(requested_norm_type == MAT_NORM_1 || requested_norm_type == MAT_NORM_INF), MethodNotImplemented
00079     ,"MatrixOp::calc_norm(...): Error, This default implemenation can only "
00080     "compute the one norm or the infinity norm!"
00081     );
00082   //
00083   // Here we implement Algorithm 2.5 in "Applied Numerical Linear Algebra", Demmel (1997)
00084   // using the momenclature in the text.  This is applied to the inverse matrix.
00085   //
00086   const MatrixOpNonsing
00087     &B = *this;
00088   bool
00089     do_trans = requested_norm_type == MAT_NORM_INF;
00090   VectorSpace::vec_mut_ptr_t
00091     x    = (do_trans ? space_rows : space_cols).create_member(1.0/num_cols),
00092     w    = (do_trans ? space_cols : space_rows).create_member(),
00093     zeta = (do_trans ? space_cols : space_rows).create_member(),
00094     z    = (do_trans ? space_rows : space_cols).create_member();
00095   const index_type max_iter = 5;  // Recommended by Highman 1988, (see Demmel's reference)
00096   value_type w_nrm = 0.0;
00097   for( index_type k = 0; k <= max_iter; ++k ) {
00098     V_InvMtV( w.get(), B, !do_trans ? no_trans : trans, *x );     // w = B*x
00099     sign( *w, zeta.get() );                                       // zeta = sign(w)
00100     V_InvMtV( z.get(), B, !do_trans ? trans : no_trans, *zeta );  // z = B'*zeta
00101     value_type  z_j = 0.0;                                        // max |z(j)| = ||z||inf
00102     index_type  j   = 0;
00103     max_abs_ele( *z, &z_j, &j );
00104     const value_type zTx = dot(*z,*x);                            // z'*x
00105     w_nrm = w->norm_1();                                        // ||w||1
00106     if( ::fabs(z_j) <= zTx ) {                                    // Update
00107       break;
00108     }
00109     else {
00110       *x = 0.0;
00111       x->set_ele(j,1.0);
00112     }
00113   }
00114   const MatNorm M_nrm = this->calc_norm(requested_norm_type);
00115   return MatNorm( w_nrm * M_nrm.value ,requested_norm_type );
00116 }
00117 
00118 // Overridden from MatrixOp
00119 
00120 MatrixOpNonsing::mat_mut_ptr_t
00121 MatrixOpNonsing::clone()
00122 {
00123   return clone_mwons();
00124 }
00125 
00126 MatrixOpNonsing::mat_ptr_t
00127 MatrixOpNonsing::clone() const
00128 {
00129   return clone_mwons();
00130 }
00131 
00132 // Overridden from MatrixNonsing
00133 
00134 MatrixOpNonsing::mat_mns_mut_ptr_t
00135 MatrixOpNonsing::clone_mns()
00136 {
00137   return clone_mwons();
00138 }
00139 
00140 MatrixOpNonsing::mat_mns_ptr_t
00141 MatrixOpNonsing::clone_mns() const
00142 {
00143   return clone_mwons();
00144 }
00145 
00146 } // end namespace AbstractLinAlgPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends