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: IpSumMatrix.cpp 759 2006-07-07 03:07:08Z andreasw $
6 //
7 // Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8 
9 #include "IpoptConfig.h"
10 #include "IpSumMatrix.hpp"
11 
12 #ifdef HAVE_CSTDIO
13 # include <cstdio>
14 #else
15 # ifdef HAVE_STDIO_H
16 #  include <stdio.h>
17 # else
18 #  error "don't have header file for stdio"
19 # endif
20 #endif
21 
22 
23 // Keeps MS VC++ 8 quiet about sprintf, strcpy, etc.
24 #ifdef _MSC_VER
25 #pragma warning(disable:4996)
26 #endif
27 
28 namespace SimTKIpopt
29 {
30 
SumMatrix(const SumMatrixSpace * owner_space)31   SumMatrix::SumMatrix(const SumMatrixSpace* owner_space)
32       :
33       Matrix(owner_space),
34       factors_(owner_space->NTerms(), 1.0),
35       matrices_(owner_space->NTerms()),
36       owner_space_(owner_space)
37   {}
38 
~SumMatrix()39   SumMatrix::~SumMatrix()
40   {}
41 
SetTerm(Index iterm,Number factor,const Matrix & matrix)42   void SumMatrix::SetTerm(Index iterm, Number factor,
43                           const Matrix& matrix)
44   {
45     DBG_ASSERT(iterm<owner_space_->NTerms());
46     factors_[iterm] = factor;
47     matrices_[iterm] = &matrix;
48   }
49 
GetTerm(Index iterm,Number & factor,SmartPtr<const Matrix> & matrix) const50   void SumMatrix::GetTerm(Index iterm, Number& factor, SmartPtr<const Matrix>& matrix) const
51   {
52     DBG_ASSERT(iterm<owner_space_->NTerms());
53     factor = factors_[iterm];
54     matrix = matrices_[iterm];
55   }
56 
NTerms() const57   Index SumMatrix::NTerms() const
58   {
59     return owner_space_->NTerms();
60   }
61 
MultVectorImpl(Number alpha,const Vector & x,Number beta,Vector & y) const62   void SumMatrix::MultVectorImpl(Number alpha, const Vector &x,
63                                  Number beta, Vector &y) const
64   {
65     //  A few sanity checks
66     DBG_ASSERT(NCols()==x.Dim());
67     DBG_ASSERT(NRows()==y.Dim());
68 
69     // Take care of the y part of the addition
70     if( beta!=0.0 ) {
71       y.Scal(beta);
72     }
73     else {
74       y.Set(0.0);  // In case y hasn't been initialized yet
75     }
76 
77     for (Index iterm=0; iterm<NTerms(); iterm++) {
78       DBG_ASSERT(IsValid(matrices_[iterm]));
79       matrices_[iterm]->MultVector(alpha*factors_[iterm], x,
80                                    1.0, y);
81     }
82   }
83 
TransMultVectorImpl(Number alpha,const Vector & x,Number beta,Vector & y) const84   void SumMatrix::TransMultVectorImpl(Number alpha, const Vector& x,
85                                       Number beta, Vector& y) const
86   {
87     //  A few sanity checks
88     DBG_ASSERT(NRows()==x.Dim());
89     DBG_ASSERT(NCols()==y.Dim());
90 
91     // Take care of the y part of the addition
92     if( beta!=0.0 ) {
93       y.Scal(beta);
94     }
95     else {
96       y.Set(0.0);  // In case y hasn't been initialized yet
97     }
98 
99     for (Index iterm=0; iterm<NTerms(); iterm++) {
100       DBG_ASSERT(IsValid(matrices_[iterm]));
101       matrices_[iterm]->TransMultVector(alpha*factors_[iterm], x,
102                                         1.0, y);
103     }
104   }
105 
HasValidNumbersImpl() const106   bool SumMatrix::HasValidNumbersImpl() const
107   {
108     for (Index iterm=0; iterm<NTerms(); iterm++) {
109       DBG_ASSERT(IsValid(matrices_[iterm]));
110       if (!matrices_[iterm]->HasValidNumbers()) {
111         return false;
112       }
113     }
114     return true;
115   }
116 
PrintImpl(const Journalist & jnlst,EJournalLevel level,EJournalCategory category,const std::string & name,Index indent,const std::string & prefix) const117   void SumMatrix::PrintImpl(const Journalist& jnlst,
118                             EJournalLevel level,
119                             EJournalCategory category,
120                             const std::string& name,
121                             Index indent,
122                             const std::string& prefix) const
123   {
124     jnlst.Printf(level, category, "\n");
125     jnlst.PrintfIndented(level, category, indent,
126                          "%sSumMatrix \"%s\" of dimension %d x %d with %d terms:\n",
127                          prefix.c_str(), name.c_str(), NRows(), NCols(), NTerms());
128     for (Index iterm=0; iterm<NTerms(); iterm++) {
129       jnlst.PrintfIndented(level, category, indent,
130                            "%sTerm %d with factor %23.16e and the following matrix:\n",
131                            prefix.c_str(), iterm, factors_[iterm]);
132       char buffer[256];
133       sprintf(buffer, "Term: %d", iterm);
134       std::string name = buffer;
135       matrices_[iterm]->Print(&jnlst, level, category, name, indent+1, prefix);
136     }
137   }
138 
SetTermSpace(Index term_idx,const MatrixSpace & mat_space)139   void SumMatrixSpace::SetTermSpace(Index term_idx, const MatrixSpace& mat_space)
140   {
141     while(term_idx >= (Index)term_spaces_.size()) {
142       term_spaces_.push_back(NULL);
143     }
144     term_spaces_[term_idx] = &mat_space;
145   }
146 
GetTermSpace(Index term_idx) const147   SmartPtr<const MatrixSpace> SumMatrixSpace::GetTermSpace(Index term_idx) const
148   {
149     if (term_idx >= 0 && term_idx < (Index)term_spaces_.size()) {
150       return term_spaces_[term_idx];
151     }
152     return NULL;
153   }
154 
MakeNewSumMatrix() const155   SumMatrix* SumMatrixSpace::MakeNewSumMatrix() const
156   {
157     DBG_ASSERT(nterms_ == (Index)term_spaces_.size());
158     return new SumMatrix(this);
159   }
160 
MakeNew() const161   Matrix* SumMatrixSpace::MakeNew() const
162   {
163     return MakeNewSumMatrix();
164   }
165 } // namespace Ipopt
166