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

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