1 // Copyright (C) 2005, 2006 International Business Machines and others.
2 // All Rights Reserved.
3 // This code is published under the Common Public License.
4 //
5 // $Id: IpDenseSymMatrix.cpp 759 2006-07-07 03:07:08Z andreasw $
6 //
7 // Authors:  Andreas Waechter             IBM    2005-12-25
8 
9 #include "IpDenseSymMatrix.hpp"
10 #include "IpDenseVector.hpp"
11 #include "IpDenseGenMatrix.hpp"
12 #include "IpBlas.hpp"
13 
14 namespace SimTKIpopt
15 {
16 
17 #ifdef IP_DEBUG
18   static const Index dbg_verbosity = 0;
19 #endif
20 
DenseSymMatrix(const DenseSymMatrixSpace * owner_space)21   DenseSymMatrix::DenseSymMatrix(const DenseSymMatrixSpace* owner_space)
22       :
23       SymMatrix(owner_space),
24       owner_space_(owner_space),
25       values_(new Number[NCols()*NRows()]),
26       initialized_(false)
27   {}
28 
~DenseSymMatrix()29   DenseSymMatrix::~DenseSymMatrix()
30   {
31     delete [] values_;
32   }
33 
MultVectorImpl(Number alpha,const Vector & x,Number beta,Vector & y) const34   void DenseSymMatrix::MultVectorImpl(Number alpha, const Vector &x,
35                                       Number beta, Vector &y) const
36   {
37     //  A few sanity checks
38     DBG_ASSERT(NCols()==x.Dim());
39     DBG_ASSERT(NRows()==y.Dim());
40     DBG_ASSERT(initialized_);
41 
42     // See if we can understand the data
43     const DenseVector* dense_x = dynamic_cast<const DenseVector*>(&x);
44     DBG_ASSERT(dense_x); /* ToDo: Implement others */
45     DenseVector* dense_y = dynamic_cast<DenseVector*>(&y);
46     DBG_ASSERT(dense_y); /* ToDo: Implement others */
47 
48     IpBlasDsymv(Dim(), alpha, values_, NRows(),
49                 dense_x->Values(), 1, beta, dense_y->Values(), 1);
50   }
51 
FillIdentity(Number factor)52   void DenseSymMatrix::FillIdentity(Number factor /*=1.*/)
53   {
54     const Index dim = Dim();
55     for (Index j=0; j<dim; j++) {
56       values_[j + j*dim] = factor;
57       for (Index i=j+1; i<dim; i++) {
58         values_[i + j*dim] = 0.;
59       }
60     }
61     ObjectChanged();
62     initialized_ = true;
63   }
64 
AddMatrix(Number alpha,const DenseSymMatrix & A,Number beta)65   void DenseSymMatrix::AddMatrix(Number alpha, const DenseSymMatrix& A,
66                                  Number beta)
67   {
68     DBG_ASSERT(beta==0. || initialized_);
69     DBG_ASSERT(Dim()==A.Dim());
70 
71     if (alpha==0.)
72       return;
73 
74     const Number* Avalues = A.Values();
75     const Index dim = Dim();
76     if (beta==0.) {
77       for (Index j=0; j<dim; j++) {
78         for (Index i=j; i<dim; i++) {
79           values_[i+j*dim] = alpha*Avalues[i+j*dim];
80         }
81       }
82     }
83     else if (beta==1.) {
84       for (Index j=0; j<dim; j++) {
85         for (Index i=j; i<dim; i++) {
86           values_[i+j*dim] += alpha*Avalues[i+j*dim];
87         }
88       }
89     }
90     else {
91       for (Index j=0; j<dim; j++) {
92         for (Index i=j; i<dim; i++) {
93           values_[i+j*dim] = alpha*Avalues[i+j*dim] + beta*values_[i+j*dim];
94         }
95       }
96     }
97     ObjectChanged();
98     initialized_ = true;
99   }
100 
HighRankUpdateTranspose(Number alpha,const MultiVectorMatrix & V1,const MultiVectorMatrix & V2,Number beta)101   void DenseSymMatrix::HighRankUpdateTranspose(Number alpha,
102       const MultiVectorMatrix& V1,
103       const MultiVectorMatrix& V2,
104       Number beta)
105   {
106     DBG_ASSERT(Dim()==V1.NCols());
107     DBG_ASSERT(Dim()==V2.NCols());
108     DBG_ASSERT(beta==0. || initialized_);
109 
110     const Index dim = Dim();
111     if (beta==0.) {
112       for (Index j=0; j<dim; j++) {
113         for (Index i=j; i<dim; i++) {
114           values_[i+j*dim] = alpha*V1.GetVector(i)->Dot(*V2.GetVector(j));
115         }
116       }
117     }
118     else {
119       for (Index j=0; j<dim; j++) {
120         for (Index i=j; i<dim; i++) {
121           values_[i+j*dim] = alpha*V1.GetVector(i)->Dot(*V2.GetVector(j))
122                              + beta*values_[i+j*dim];
123         }
124       }
125     }
126     initialized_ = true;
127     ObjectChanged();
128   }
129 
HighRankUpdate(bool trans,Number alpha,const DenseGenMatrix & V,Number beta)130   void DenseSymMatrix::HighRankUpdate(bool trans, Number alpha,
131                                       const DenseGenMatrix& V,
132                                       Number beta)
133   {
134     DBG_ASSERT((!trans && Dim()==V.NRows()) || (trans && Dim()==V.NCols()));
135     DBG_ASSERT(beta==0. || initialized_);
136 
137     Index nrank;
138     if (trans) {
139       nrank = V.NRows();
140     }
141     else {
142       nrank = V.NCols();
143     }
144 
145     IpBlasDsyrk(trans, Dim(), nrank, alpha, V.Values(), V.NRows(),
146                 beta, values_, NRows());
147 
148     initialized_ = true;
149     ObjectChanged();
150   }
151 
SpecialAddForLMSR1(const DenseVector & D,const DenseGenMatrix & L)152   void DenseSymMatrix::SpecialAddForLMSR1(const DenseVector& D,
153                                           const DenseGenMatrix& L)
154   {
155     const Index dim = Dim();
156     DBG_ASSERT(initialized_);
157     DBG_ASSERT(dim==D.Dim());
158     DBG_ASSERT(dim==L.NRows());
159     DBG_ASSERT(dim==L.NCols());
160 
161     // First add the diagonal matrix
162     const Number* Dvalues = D.Values();
163     for (Index i=0; i<dim; i++) {
164       values_[i+i*dim] += Dvalues[i];
165     }
166 
167     // Now add the strictly-lower triagular matrix L and its transpose
168     const Number* Lvalues = L.Values();
169     for (Index j=0; j<dim; j++) {
170       for (Index i=j+1; i<dim; i++) {
171         values_[i+j*dim] += Lvalues[i+j*dim];
172       }
173     }
174     ObjectChanged();
175   }
176 
HasValidNumbersImpl() const177   bool DenseSymMatrix::HasValidNumbersImpl() const
178   {
179     DBG_ASSERT(initialized_);
180     Number sum = 0.;
181     const Index dim = Dim();
182     for (Index j=0; j<dim; j++) {
183       sum += values_[j + j*dim];
184       for (Index i=j+1; i<dim; i++) {
185         sum += values_[i + j*dim];
186       }
187     }
188     return IsFiniteNumber(sum);
189   }
190 
PrintImpl(const Journalist & jnlst,EJournalLevel level,EJournalCategory category,const std::string & name,Index indent,const std::string & prefix) const191   void DenseSymMatrix::PrintImpl(const Journalist& jnlst,
192                                  EJournalLevel level,
193                                  EJournalCategory category,
194                                  const std::string& name,
195                                  Index indent,
196                                  const std::string& prefix) const
197   {
198     jnlst.Printf(level, category, "\n");
199     jnlst.PrintfIndented(level, category, indent,
200                          "%sDenseSymMatrix \"%s\" of dimension %d (only lower triangular part printed):\n",
201                          prefix.c_str(), name.c_str(), Dim());
202 
203     if (initialized_) {
204       for (Index j=0; j<NCols(); j++) {
205         for (Index i=j; i<NRows(); i++) {
206           jnlst.PrintfIndented(level, category, indent,
207                                "%s%s[%5d,%5d]=%23.16e\n",
208                                prefix.c_str(), name.c_str(), i, j, values_[i+NRows()*j]);
209         }
210       }
211     }
212     else {
213       jnlst.PrintfIndented(level, category, indent,
214                            "The matrix has not yet been initialized!\n");
215     }
216   }
217 
DenseSymMatrixSpace(Index nDim)218   DenseSymMatrixSpace::DenseSymMatrixSpace(Index nDim)
219       :
220       SymMatrixSpace(nDim)
221   {}
222 
223 } // namespace Ipopt
224