1 // Copyright (C) 2004, 2006 International Business Machines and others. 2 // All Rights Reserved. 3 // This code is published under the Common Public License. 4 // 5 // $Id: IpScaledMatrix.cpp 759 2006-07-07 03:07:08Z andreasw $ 6 // 7 // Authors: Carl Laird, Andreas Waechter IBM 2004-08-13 8 9 #include "IpScaledMatrix.hpp" 10 11 namespace SimTKIpopt 12 { 13 ScaledMatrix(const ScaledMatrixSpace * owner_space)14 ScaledMatrix::ScaledMatrix(const ScaledMatrixSpace* owner_space) 15 : 16 Matrix(owner_space), 17 owner_space_(owner_space) 18 {} 19 20 ~ScaledMatrix()21 ScaledMatrix::~ScaledMatrix() 22 {} 23 MultVectorImpl(Number alpha,const Vector & x,Number beta,Vector & y) const24 void ScaledMatrix::MultVectorImpl(Number alpha, const Vector &x, 25 Number beta, Vector &y) const 26 { 27 DBG_ASSERT(IsValid(matrix_)); 28 29 // Take care of the y part of the addition 30 if( beta!=0.0 ) { 31 y.Scal(beta); 32 } 33 else { 34 y.Set(0.0); // In case y hasn't been initialized yet 35 } 36 37 // need some temporary vectors 38 SmartPtr<Vector> tmp_x = x.MakeNewCopy(); 39 SmartPtr<Vector> tmp_y = y.MakeNew(); 40 41 if (IsValid(owner_space_->ColumnScaling())) { 42 tmp_x->ElementWiseMultiply(*owner_space_->ColumnScaling()); 43 } 44 45 matrix_->MultVector(1.0, *tmp_x, 0.0, *tmp_y); 46 47 if (IsValid(owner_space_->RowScaling())) { 48 tmp_y->ElementWiseMultiply(*owner_space_->RowScaling()); 49 } 50 51 y.Axpy(1.0, *tmp_y); 52 } 53 TransMultVectorImpl(Number alpha,const Vector & x,Number beta,Vector & y) const54 void ScaledMatrix::TransMultVectorImpl(Number alpha, const Vector &x, 55 Number beta, Vector &y) const 56 { 57 DBG_ASSERT(IsValid(matrix_)); 58 59 // Take care of the y part of the addition 60 if( beta!=0.0 ) { 61 y.Scal(beta); 62 } 63 else { 64 y.Set(0.0); // In case y hasn't been initialized yet 65 } 66 67 // need some temporary vectors 68 SmartPtr<Vector> tmp_x = x.MakeNewCopy(); 69 SmartPtr<Vector> tmp_y = y.MakeNew(); 70 71 if (IsValid(owner_space_->RowScaling())) { 72 tmp_x->ElementWiseMultiply(*owner_space_->RowScaling()); 73 } 74 75 matrix_->TransMultVector(1.0, *tmp_x, 0.0, *tmp_y); 76 77 if (IsValid(owner_space_->ColumnScaling())) { 78 tmp_y->ElementWiseMultiply(*owner_space_->ColumnScaling()); 79 } 80 81 y.Axpy(1.0, *tmp_y); 82 } 83 HasValidNumbersImpl() const84 bool ScaledMatrix::HasValidNumbersImpl() const 85 { 86 DBG_ASSERT(IsValid(matrix_)); 87 return matrix_->HasValidNumbers(); 88 } 89 PrintImpl(const Journalist & jnlst,EJournalLevel level,EJournalCategory category,const std::string & name,Index indent,const std::string & prefix) const90 void ScaledMatrix::PrintImpl(const Journalist& jnlst, 91 EJournalLevel level, 92 EJournalCategory category, 93 const std::string& name, 94 Index indent, 95 const std::string& prefix) const 96 { 97 jnlst.Printf(level, category, "\n"); 98 jnlst.PrintfIndented(level, category, indent, 99 "%sScaledMatrix \"%s\" of dimension %d x %d:\n", 100 prefix.c_str(), name.c_str(), NRows(), NCols()); 101 if (IsValid(owner_space_->RowScaling())) { 102 owner_space_->RowScaling()->Print(&jnlst, level, category, 103 name+"_row_scaling", indent+1, prefix); 104 } 105 else { 106 jnlst.PrintfIndented(level, category, indent+1, "RowScaling is NULL\n"); 107 } 108 if (IsValid(matrix_)) { 109 matrix_->Print(&jnlst, level, category, name+"_unscaled_matrix", 110 indent+1, prefix); 111 } 112 else { 113 jnlst.PrintfIndented(level, category, indent+1, 114 "%sunscaled matrix is NULL\n", prefix.c_str()); 115 } 116 if (IsValid(owner_space_->ColumnScaling())) { 117 owner_space_->ColumnScaling()->Print(&jnlst, level, category, 118 name+"_column_scaling", 119 indent+1, prefix); 120 } 121 else { 122 jnlst.PrintfIndented(level, category, indent+1, 123 "%sColumnScaling is NULL\n", prefix.c_str()); 124 } 125 } 126 AddMSinvZImpl(Number alpha,const Vector & S,const Vector & Z,Vector & X) const127 void ScaledMatrix::AddMSinvZImpl(Number alpha, const Vector& S, 128 const Vector& Z, Vector& X) const 129 { 130 DBG_ASSERT(false && "Got the ScaledMatrix::AddMSinvZImpl. Should implement specialized method!"); 131 132 SmartPtr<Vector> tmp = S.MakeNew(); 133 tmp->AddVectorQuotient(1., Z, S, 0.); 134 MultVector(alpha, *tmp, 1., X); 135 } 136 SinvBlrmZMTdBrImpl(Number alpha,const Vector & S,const Vector & R,const Vector & Z,const Vector & D,Vector & X) const137 void ScaledMatrix::SinvBlrmZMTdBrImpl(Number alpha, const Vector& S, 138 const Vector& R, const Vector& Z, 139 const Vector& D, Vector& X) const 140 { 141 DBG_ASSERT(false && "Got the ScaledMatrix::SinvBlrmZMTdBrImpl. Should implement specialized method!"); 142 143 TransMultVector(alpha, D, 0., X); 144 X.ElementWiseMultiply(Z); 145 X.Axpy(1., R); 146 X.ElementWiseDivide(S); 147 } 148 149 ScaledMatrixSpace(const SmartPtr<const Vector> & row_scaling,bool row_scaling_reciprocal,const SmartPtr<const MatrixSpace> & unscaled_matrix_space,const SmartPtr<const Vector> & column_scaling,bool column_scaling_reciprocal)150 ScaledMatrixSpace::ScaledMatrixSpace( 151 const SmartPtr<const Vector>& row_scaling, 152 bool row_scaling_reciprocal, 153 const SmartPtr<const MatrixSpace>& unscaled_matrix_space, 154 const SmartPtr<const Vector>& column_scaling, 155 bool column_scaling_reciprocal) 156 : 157 MatrixSpace(unscaled_matrix_space->NRows(), 158 unscaled_matrix_space->NCols()), 159 unscaled_matrix_space_(unscaled_matrix_space) 160 { 161 if (IsValid(row_scaling)) { 162 row_scaling_ = row_scaling->MakeNewCopy(); 163 if (row_scaling_reciprocal) { 164 row_scaling_->ElementWiseReciprocal(); 165 } 166 } 167 else { 168 row_scaling_ = NULL; 169 } 170 171 if (IsValid(column_scaling)) { 172 column_scaling_ = column_scaling->MakeNewCopy(); 173 if (column_scaling_reciprocal) { 174 column_scaling_->ElementWiseReciprocal(); 175 } 176 } 177 else { 178 column_scaling_ = NULL; 179 } 180 } 181 } // namespace Ipopt 182