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