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