1 // Copyright (C) 2004, 2009 International Business Machines and others. 2 // All Rights Reserved. 3 // This code is published under the Eclipse Public License. 4 // 5 // $Id$ 6 // 7 // Authors: Carl Laird, Andreas Waechter IBM 2004-08-13 8 9 #include "IpoptConfig.h" 10 #include "IpCompoundMatrix.hpp" 11 #include "IpCompoundVector.hpp" 12 13 #ifdef HAVE_CSTDIO 14 # include <cstdio> 15 #else 16 # ifdef HAVE_STDIO_H 17 # include <stdio.h> 18 # else 19 # error "don't have header file for stdio" 20 # endif 21 #endif 22 23 #define ALLOW_NESTED 24 25 namespace Ipopt 26 { 27 CompoundMatrix(const CompoundMatrixSpace * owner_space)28 CompoundMatrix::CompoundMatrix(const CompoundMatrixSpace* owner_space) 29 : 30 Matrix(owner_space), 31 owner_space_(owner_space), 32 matrices_valid_(false) 33 { 34 std::vector<SmartPtr<Matrix> > row(NComps_Cols()); 35 std::vector<SmartPtr<const Matrix> > const_row(NComps_Cols()); 36 for (Index irow=0; irow<NComps_Rows(); irow++) { 37 const_comps_.push_back(const_row); 38 comps_.push_back(row); 39 } 40 } 41 42 ~CompoundMatrix()43 CompoundMatrix::~CompoundMatrix() 44 {} 45 SetComp(Index irow,Index jcol,const Matrix & matrix)46 void CompoundMatrix::SetComp(Index irow, Index jcol, const Matrix& matrix) 47 { 48 DBG_ASSERT(irow<NComps_Rows()); 49 DBG_ASSERT(jcol<NComps_Cols()); 50 DBG_ASSERT(owner_space_->GetCompSpace(irow, jcol)->IsMatrixFromSpace(matrix)); 51 52 comps_[irow][jcol] = NULL; 53 const_comps_[irow][jcol] = &matrix; 54 ObjectChanged(); 55 } 56 SetCompNonConst(Index irow,Index jcol,Matrix & matrix)57 void CompoundMatrix::SetCompNonConst(Index irow, Index jcol, Matrix& matrix) 58 { 59 DBG_ASSERT(irow < NComps_Rows()); 60 DBG_ASSERT(jcol < NComps_Cols()); 61 DBG_ASSERT(owner_space_->GetCompSpace(irow, jcol)->IsMatrixFromSpace(matrix)); 62 63 const_comps_[irow][jcol] = NULL; 64 comps_[irow][jcol] = &matrix; 65 ObjectChanged(); 66 } 67 CreateBlockFromSpace(Index irow,Index jcol)68 void CompoundMatrix::CreateBlockFromSpace(Index irow, Index jcol) 69 { 70 DBG_ASSERT(irow < NComps_Rows()); 71 DBG_ASSERT(jcol < NComps_Cols()); 72 DBG_ASSERT(IsValid(owner_space_->GetCompSpace(irow, jcol))); 73 SetCompNonConst(irow, jcol, *owner_space_->GetCompSpace(irow,jcol)->MakeNew()); 74 } 75 MultVectorImpl(Number alpha,const Vector & x,Number beta,Vector & y) const76 void CompoundMatrix::MultVectorImpl(Number alpha, const Vector &x, 77 Number beta, Vector &y) const 78 { 79 if (!matrices_valid_) { 80 matrices_valid_ = MatricesValid(); 81 } 82 DBG_ASSERT(matrices_valid_); 83 84 // The vectors are assumed to be compound Vectors as well 85 const CompoundVector* comp_x = dynamic_cast<const CompoundVector*>(&x); 86 CompoundVector* comp_y = dynamic_cast<CompoundVector*>(&y); 87 88 #ifndef ALLOW_NESTED 89 // A few sanity checks 90 if (comp_x) { 91 DBG_ASSERT(NComps_Cols()==comp_x->NComps()); 92 } 93 else { 94 DBG_ASSERT(NComps_Cols() == 1); 95 } 96 97 if (comp_y) { 98 DBG_ASSERT(NComps_Rows()==comp_y->NComps()); 99 } 100 else { 101 DBG_ASSERT(NComps_Rows() == 1); 102 } 103 #endif 104 if (comp_x) { 105 if (NComps_Cols()!=comp_x->NComps()) { 106 comp_x = NULL; 107 } 108 } 109 if (comp_y) { 110 if (NComps_Rows()!=comp_y->NComps()) { 111 comp_y = NULL; 112 } 113 } 114 115 // Take care of the y part of the addition 116 if ( beta!=0.0 ) { 117 y.Scal(beta); 118 } 119 else { 120 y.Set(0.0); // In case y hasn't been initialized yet 121 } 122 123 for ( Index irow = 0; irow < NComps_Rows(); irow++ ) { 124 SmartPtr<Vector> y_i; 125 if (comp_y) { 126 y_i = comp_y->GetCompNonConst(irow); 127 } 128 else { 129 y_i = &y; 130 } 131 DBG_ASSERT(IsValid(y_i)); 132 133 for ( Index jcol = 0; jcol < NComps_Cols(); jcol++ ) { 134 if ( (owner_space_->Diagonal() && irow == jcol) 135 || (!owner_space_->Diagonal() && ConstComp(irow,jcol)) ) { 136 SmartPtr<const Vector> x_j; 137 if (comp_x) { 138 x_j = comp_x->GetComp(jcol); 139 } 140 else if (NComps_Cols() == 1) { 141 x_j = &x; 142 } 143 DBG_ASSERT(IsValid(x_j)); 144 145 ConstComp(irow, jcol)->MultVector(alpha, *x_j, 146 1., *y_i); 147 } 148 } 149 } 150 } 151 TransMultVectorImpl(Number alpha,const Vector & x,Number beta,Vector & y) const152 void CompoundMatrix::TransMultVectorImpl(Number alpha, const Vector &x, 153 Number beta, Vector &y) const 154 { 155 if (!matrices_valid_) { 156 matrices_valid_ = MatricesValid(); 157 } 158 DBG_ASSERT(matrices_valid_); 159 160 // The vectors are assumed to be compound Vectors as well 161 const CompoundVector* comp_x = dynamic_cast<const CompoundVector*>(&x); 162 CompoundVector* comp_y = dynamic_cast<CompoundVector*>(&y); 163 164 #ifndef ALLOW_NESTED 165 // A few sanity checks 166 if (comp_y) { 167 DBG_ASSERT(NComps_Cols()==comp_y->NComps()); 168 } 169 else { 170 DBG_ASSERT(NComps_Cols() == 1); 171 } 172 173 if (comp_x) { 174 DBG_ASSERT(NComps_Rows()==comp_x->NComps()); 175 } 176 else { 177 DBG_ASSERT(NComps_Rows() == 1); 178 } 179 #endif 180 if (comp_y) { 181 if (NComps_Cols()!=comp_y->NComps()) { 182 comp_y = NULL; 183 } 184 } 185 if (comp_x) { 186 if (NComps_Rows()!=comp_x->NComps()) { 187 comp_x = NULL; 188 } 189 } 190 191 // Take care of the y part of the addition 192 if ( beta!=0.0 ) { 193 y.Scal(beta); 194 } 195 else { 196 y.Set(0.0); // In case y hasn't been initialized yet 197 } 198 199 for ( Index irow = 0; irow < NComps_Cols(); irow++ ) { 200 SmartPtr<Vector> y_i; 201 if (comp_y) { 202 y_i = comp_y->GetCompNonConst(irow); 203 } 204 else { 205 y_i = &y; 206 } 207 DBG_ASSERT(IsValid(y_i)); 208 209 for ( Index jcol = 0; jcol < NComps_Rows(); jcol++ ) { 210 if ( (owner_space_->Diagonal() && irow == jcol) 211 || (!owner_space_->Diagonal() && ConstComp(jcol, irow)) ) { 212 SmartPtr<const Vector> x_j; 213 if (comp_x) { 214 x_j = comp_x->GetComp(jcol); 215 } 216 else { 217 x_j = &x; 218 } 219 DBG_ASSERT(IsValid(x_j)); 220 221 ConstComp(jcol, irow)->TransMultVector(alpha, *x_j, 222 1., *y_i); 223 } 224 } 225 } 226 } 227 228 // Specialized method (overloaded from IpMatrix) AddMSinvZImpl(Number alpha,const Vector & S,const Vector & Z,Vector & X) const229 void CompoundMatrix::AddMSinvZImpl(Number alpha, const Vector& S, 230 const Vector& Z, Vector& X) const 231 { 232 // The vectors are assumed to be compound Vectors as well (unless they 233 // are assumed to consist of only one component 234 const CompoundVector* comp_S = dynamic_cast<const CompoundVector*>(&S); 235 const CompoundVector* comp_Z = dynamic_cast<const CompoundVector*>(&Z); 236 CompoundVector* comp_X = dynamic_cast<CompoundVector*>(&X); 237 238 #ifndef ALLOW_NESTED 239 // A few sanity checks for sizes 240 if (comp_S) { 241 DBG_ASSERT(NComps_Cols()==comp_S->NComps()); 242 } 243 else { 244 DBG_ASSERT(NComps_Cols() == 1); 245 } 246 if (comp_Z) { 247 DBG_ASSERT(NComps_Cols()==comp_Z->NComps()); 248 } 249 else { 250 DBG_ASSERT(NComps_Cols() == 1); 251 } 252 if (comp_X) { 253 DBG_ASSERT(NComps_Rows()==comp_X->NComps()); 254 } 255 else { 256 DBG_ASSERT(NComps_Rows() == 1); 257 } 258 #endif 259 if (comp_S) { 260 if (NComps_Cols()!=comp_S->NComps()) { 261 comp_S = NULL; 262 } 263 } 264 if (comp_Z) { 265 if (NComps_Cols()!=comp_Z->NComps()) { 266 comp_Z = NULL; 267 } 268 } 269 if (comp_X) { 270 if (NComps_Rows()!=comp_X->NComps()) { 271 comp_X = NULL; 272 } 273 } 274 275 for ( Index irow = 0; irow < NComps_Rows(); irow++ ) { 276 SmartPtr<Vector> X_i; 277 if (comp_X) { 278 X_i = comp_X->GetCompNonConst(irow); 279 } 280 else { 281 X_i = &X; 282 } 283 DBG_ASSERT(IsValid(X_i)); 284 285 for ( Index jcol = 0; jcol < NComps_Cols(); jcol++ ) { 286 if ( (owner_space_->Diagonal() && irow == jcol) 287 || (!owner_space_->Diagonal() && ConstComp(irow, jcol)) ) { 288 SmartPtr<const Vector> S_j; 289 if (comp_S) { 290 S_j = comp_S->GetComp(jcol); 291 } 292 else { 293 S_j = &S; 294 } 295 DBG_ASSERT(IsValid(S_j)); 296 SmartPtr<const Vector> Z_j; 297 if (comp_Z) { 298 Z_j = comp_Z->GetComp(jcol); 299 } 300 else { 301 Z_j = &Z; 302 } 303 DBG_ASSERT(IsValid(Z_j)); 304 305 ConstComp(irow, jcol)->AddMSinvZ(alpha, *S_j, *Z_j, *X_i); 306 } 307 } 308 } 309 } 310 311 // Specialized method (overloaded from IpMatrix) SinvBlrmZMTdBrImpl(Number alpha,const Vector & S,const Vector & R,const Vector & Z,const Vector & D,Vector & X) const312 void CompoundMatrix::SinvBlrmZMTdBrImpl(Number alpha, const Vector& S, 313 const Vector& R, const Vector& Z, 314 const Vector& D, Vector& X) const 315 { 316 // First check if the matrix is indeed such that we can use the 317 // special methods from the component spaces (this only works if 318 // we have exactly one submatrix per column) 319 320 bool fast_SinvBlrmZMTdBr = false; 321 322 if (!owner_space_->Diagonal()) { 323 fast_SinvBlrmZMTdBr = true; 324 for (Index jcol=0; jcol < NComps_Cols(); jcol++ ) { 325 Index nblocks = 0; 326 for (Index irow=0; irow < NComps_Rows(); irow++ ) { 327 if (ConstComp(irow, jcol)) { 328 nblocks++; 329 if (nblocks>1) { 330 break; 331 } 332 } 333 } 334 if (nblocks!=1) { 335 fast_SinvBlrmZMTdBr = false; 336 break; 337 } 338 } 339 } 340 341 if (!owner_space_->Diagonal() && !fast_SinvBlrmZMTdBr) { 342 // Use the standard replacement implementation 343 Matrix::SinvBlrmZMTdBrImpl(alpha, S, R, Z, D, X); 344 DBG_ASSERT(false && "Found a matrix where we can't use the fast SinvBlrmZMTdBr implementation in CompoundMatrix"); 345 } 346 else { 347 // The vectors are assumed to be compound Vectors as well (unless they 348 // are assumed to consist of only one component 349 const CompoundVector* comp_S = dynamic_cast<const CompoundVector*>(&S); 350 const CompoundVector* comp_R = dynamic_cast<const CompoundVector*>(&R); 351 const CompoundVector* comp_Z = dynamic_cast<const CompoundVector*>(&Z); 352 const CompoundVector* comp_D = dynamic_cast<const CompoundVector*>(&D); 353 CompoundVector* comp_X = dynamic_cast<CompoundVector*>(&X); 354 355 #ifndef ALLOW_NESTED 356 // A few sanity checks for sizes 357 if (comp_S) { 358 DBG_ASSERT(NComps_Cols()==comp_S->NComps()); 359 } 360 else { 361 DBG_ASSERT(NComps_Cols() == 1); 362 } 363 if (comp_Z) { 364 DBG_ASSERT(NComps_Cols()==comp_Z->NComps()); 365 } 366 else { 367 DBG_ASSERT(NComps_Cols() == 1); 368 } 369 if (comp_R) { 370 DBG_ASSERT(NComps_Cols()==comp_R->NComps()); 371 } 372 else { 373 DBG_ASSERT(NComps_Cols() == 1); 374 } 375 if (comp_D) { 376 DBG_ASSERT(NComps_Rows()==comp_D->NComps()); 377 } 378 else { 379 DBG_ASSERT(NComps_Rows() == 1); 380 } 381 if (comp_X) { 382 DBG_ASSERT(NComps_Cols()==comp_X->NComps()); 383 } 384 else { 385 DBG_ASSERT(NComps_Cols() == 1); 386 } 387 #endif 388 if (comp_S) { 389 if (NComps_Cols()!=comp_S->NComps()) { 390 comp_S = NULL; 391 } 392 } 393 if (comp_Z) { 394 if (NComps_Cols()!=comp_Z->NComps()) { 395 comp_Z = NULL; 396 } 397 } 398 if (comp_R) { 399 if (NComps_Cols()!=comp_R->NComps()) { 400 comp_R = NULL; 401 } 402 } 403 if (comp_D) { 404 if (NComps_Rows()!=comp_D->NComps()) { 405 comp_D = NULL; 406 } 407 } 408 if (comp_X) { 409 if (NComps_Cols()!=comp_X->NComps()) { 410 comp_X = NULL; 411 } 412 } 413 414 for (Index irow=0; irow<NComps_Cols(); irow++) { 415 Index jcol = irow; 416 if (!owner_space_->Diagonal()) { 417 for (Index j=0; j<NComps_Rows(); j++) { 418 if (ConstComp(j, irow)) { 419 jcol = j; 420 break; 421 } 422 } 423 } 424 SmartPtr<const Vector> S_i; 425 if (comp_S) { 426 S_i = comp_S->GetComp(irow); 427 } 428 else { 429 S_i = &S; 430 } 431 DBG_ASSERT(IsValid(S_i)); 432 SmartPtr<const Vector> Z_i; 433 if (comp_Z) { 434 Z_i = comp_Z->GetComp(irow); 435 } 436 else { 437 Z_i = &Z; 438 } 439 DBG_ASSERT(IsValid(Z_i)); 440 SmartPtr<const Vector> R_i; 441 if (comp_R) { 442 R_i = comp_R->GetComp(irow); 443 } 444 else { 445 R_i = &R; 446 } 447 DBG_ASSERT(IsValid(R_i)); 448 SmartPtr<const Vector> D_i; 449 if (comp_D) { 450 D_i = comp_D->GetComp(jcol); 451 } 452 else { 453 D_i = &D; 454 } 455 DBG_ASSERT(IsValid(D_i)); 456 SmartPtr<Vector> X_i; 457 if (comp_X) { 458 X_i = comp_X->GetCompNonConst(irow); 459 } 460 else { 461 X_i = &X; 462 } 463 DBG_ASSERT(IsValid(X_i)); 464 465 ConstComp(jcol, irow)->SinvBlrmZMTdBr(alpha, *S_i, *R_i, *Z_i, 466 *D_i, *X_i); 467 } 468 } 469 } 470 HasValidNumbersImpl() const471 bool CompoundMatrix::HasValidNumbersImpl() const 472 { 473 if (!matrices_valid_) { 474 matrices_valid_ = MatricesValid(); 475 } 476 DBG_ASSERT(matrices_valid_); 477 478 for ( Index irow = 0; irow < NComps_Rows(); irow++ ) { 479 for ( Index jcol = 0; jcol < NComps_Cols(); jcol++ ) { 480 if ( (owner_space_->Diagonal() && irow == jcol) 481 || (!owner_space_->Diagonal() && ConstComp(irow,jcol)) ) { 482 if (!ConstComp(irow, jcol)->HasValidNumbers()) { 483 return false; 484 } 485 } 486 } 487 } 488 return true; 489 } 490 ComputeRowAMaxImpl(Vector & rows_norms,bool init) const491 void CompoundMatrix::ComputeRowAMaxImpl(Vector& rows_norms, bool init) const 492 { 493 if (!matrices_valid_) { 494 matrices_valid_ = MatricesValid(); 495 } 496 DBG_ASSERT(matrices_valid_); 497 498 // The vector is assumed to be compound Vectors as well except if 499 // there is only one component 500 CompoundVector* comp_vec = dynamic_cast<CompoundVector*>(&rows_norms); 501 502 #ifndef ALLOW_NESTED 503 // A few sanity checks 504 if (comp_vec) { 505 DBG_ASSERT(NComps_Rows()==comp_vec->NComps()); 506 } 507 else { 508 DBG_ASSERT(NComps_Rows() == 1); 509 } 510 #endif 511 if (comp_vec) { 512 if (NComps_Rows()!=comp_vec->NComps()) { 513 comp_vec = NULL; 514 } 515 } 516 517 for (Index jcol = 0; jcol < NComps_Cols(); jcol++) { 518 for (Index irow = 0; irow < NComps_Rows(); irow++) { 519 if (ConstComp(irow, jcol)) { 520 SmartPtr<Vector> vec_i; 521 if (comp_vec) { 522 vec_i = comp_vec->GetCompNonConst(irow); 523 } 524 else { 525 vec_i = &rows_norms; 526 } 527 DBG_ASSERT(IsValid(vec_i)); 528 ConstComp(irow, jcol)->ComputeRowAMax(*vec_i, false); 529 } 530 } 531 } 532 } 533 ComputeColAMaxImpl(Vector & cols_norms,bool init) const534 void CompoundMatrix::ComputeColAMaxImpl(Vector& cols_norms, bool init) const 535 { 536 if (!matrices_valid_) { 537 matrices_valid_ = MatricesValid(); 538 } 539 DBG_ASSERT(matrices_valid_); 540 541 // The vector is assumed to be compound Vectors as well except if 542 // there is only one component 543 CompoundVector* comp_vec = dynamic_cast<CompoundVector*>(&cols_norms); 544 545 #ifndef ALLOW_NESTED 546 // A few sanity checks 547 if (comp_vec) { 548 DBG_ASSERT(NComps_Cols()==comp_vec->NComps()); 549 } 550 else { 551 DBG_ASSERT(NComps_Cols() == 1); 552 } 553 #endif 554 if (comp_vec) { 555 if (NComps_Cols()!=comp_vec->NComps()) { 556 comp_vec = NULL; 557 } 558 } 559 560 for (Index irow = 0; irow < NComps_Rows(); irow++) { 561 for (Index jcol = 0; jcol < NComps_Cols(); jcol++) { 562 if (ConstComp(irow, jcol)) { 563 SmartPtr<Vector> vec_i; 564 if (comp_vec) { 565 vec_i = comp_vec->GetCompNonConst(irow); 566 } 567 else { 568 vec_i = &cols_norms; 569 } 570 DBG_ASSERT(IsValid(vec_i)); 571 ConstComp(irow, jcol)->ComputeColAMax(*vec_i, false); 572 } 573 } 574 } 575 } 576 PrintImpl(const Journalist & jnlst,EJournalLevel level,EJournalCategory category,const std::string & name,Index indent,const std::string & prefix) const577 void CompoundMatrix::PrintImpl(const Journalist& jnlst, 578 EJournalLevel level, 579 EJournalCategory category, 580 const std::string& name, 581 Index indent, 582 const std::string& prefix) const 583 { 584 jnlst.Printf(level, category, "\n"); 585 jnlst.PrintfIndented(level, category, indent, 586 "%sCompoundMatrix \"%s\" with %d row and %d columns components:\n", 587 prefix.c_str(), name.c_str(), 588 NComps_Rows(), NComps_Cols()); 589 for (Index irow = 0; irow < NComps_Rows(); irow++ ) { 590 for (Index jcol = 0; jcol < NComps_Cols(); jcol++ ) { 591 jnlst.PrintfIndented(level, category, indent, 592 "%sComponent for row %d and column %d:\n", 593 prefix.c_str(), irow, jcol); 594 if (ConstComp(irow, jcol)) { 595 DBG_ASSERT(name.size()<200); 596 char buffer[256]; 597 Snprintf(buffer, 255, "%s[%2d][%2d]", name.c_str(), irow, jcol); 598 std::string term_name = buffer; 599 ConstComp(irow, jcol)->Print(&jnlst, level, category, term_name, 600 indent+1, prefix); 601 } 602 else { 603 jnlst.PrintfIndented(level, category, indent, 604 "%sComponent has not been set.\n", 605 prefix.c_str()); 606 } 607 } 608 } 609 } 610 MatricesValid() const611 bool CompoundMatrix::MatricesValid() const 612 { 613 // Check to make sure we have matrices everywhere the space has matrices 614 // We already check that the matrix agrees with the block space 615 // in the SetComp methods 616 bool retValue = true; 617 for (Index i=0; i<NComps_Rows(); i++) { 618 for (Index j=0; j<NComps_Cols(); j++) { 619 if ( (!ConstComp(i,j) && IsValid(owner_space_->GetCompSpace(i,j)) && 620 owner_space_->GetCompSpace(i,j)->NRows()>0 && 621 owner_space_->GetCompSpace(i,j)->NCols()>0) 622 || (ConstComp(i,j) && IsNull(owner_space_->GetCompSpace(i,j))) ) { 623 retValue = false; 624 break; 625 } 626 } 627 } 628 return retValue; 629 } 630 631 CompoundMatrixSpace(Index ncomps_rows,Index ncomps_cols,Index total_nRows,Index total_nCols)632 CompoundMatrixSpace::CompoundMatrixSpace(Index ncomps_rows, Index ncomps_cols, Index total_nRows, Index total_nCols) 633 : 634 MatrixSpace(total_nRows, total_nCols), 635 ncomps_rows_(ncomps_rows), 636 ncomps_cols_(ncomps_cols), 637 dimensions_set_(false), 638 block_rows_(ncomps_rows, -1), 639 block_cols_(ncomps_cols, -1), 640 diagonal_(false) 641 { 642 DBG_START_METH("CompoundMatrixSpace::CompoundMatrixSpace", 0); 643 std::vector<SmartPtr<const MatrixSpace> > row(ncomps_cols_); 644 std::vector< bool > allocate_row(ncomps_cols_, false); 645 for (Index i=0; i<ncomps_rows_; i++) { 646 DBG_PRINT((1, "block_rows_[%d] = %d\n", i, block_rows_[i])); 647 comp_spaces_.push_back(row); 648 allocate_block_.push_back(allocate_row); 649 } 650 } 651 SetBlockCols(Index jcol,Index ncols)652 void CompoundMatrixSpace::SetBlockCols(Index jcol, Index ncols) 653 { 654 DBG_ASSERT(!dimensions_set_ && "for now, if the dimensions have all been set, they cannot be changed"); 655 DBG_ASSERT(jcol < ncomps_cols_); 656 DBG_ASSERT(block_cols_[jcol] == -1 && "This dimension has already been set - sanity check"); 657 block_cols_[jcol] = ncols; 658 } 659 SetBlockRows(Index irow,Index nrows)660 void CompoundMatrixSpace::SetBlockRows(Index irow, Index nrows) 661 { 662 DBG_ASSERT(!dimensions_set_ && "for now, if the dimensions have all been set, they cannot be changed"); 663 DBG_ASSERT(irow < ncomps_rows_); 664 DBG_ASSERT(block_rows_[irow] == -1 && "This dimension has already been set - sanity check"); 665 block_rows_[irow] = nrows; 666 } 667 GetBlockRows(Index irow) const668 Index CompoundMatrixSpace::GetBlockRows(Index irow) const 669 { 670 DBG_ASSERT(dimensions_set_); 671 DBG_ASSERT(irow < ncomps_rows_); 672 return block_rows_[irow]; 673 } 674 GetBlockCols(Index jcol) const675 Index CompoundMatrixSpace::GetBlockCols(Index jcol) const 676 { 677 DBG_ASSERT(dimensions_set_); 678 DBG_ASSERT(jcol < ncomps_cols_); 679 return block_cols_[jcol]; 680 } 681 SetCompSpace(Index irow,Index jcol,const MatrixSpace & mat_space,bool auto_allocate)682 void CompoundMatrixSpace::SetCompSpace(Index irow, Index jcol, 683 const MatrixSpace& mat_space, 684 bool auto_allocate /*=false*/) 685 { 686 if (!dimensions_set_) { 687 dimensions_set_ = DimensionsSet(); 688 } 689 DBG_ASSERT(dimensions_set_); 690 DBG_ASSERT(irow<ncomps_rows_); 691 DBG_ASSERT(jcol<ncomps_cols_); 692 DBG_ASSERT(IsNull(comp_spaces_[irow][jcol])); 693 DBG_ASSERT(block_cols_[jcol] != -1 && block_cols_[jcol] == mat_space.NCols()); 694 DBG_ASSERT(block_rows_[irow] != -1 && block_rows_[irow] == mat_space.NRows()); 695 696 comp_spaces_[irow][jcol] = &mat_space; 697 allocate_block_[irow][jcol] = auto_allocate; 698 699 diagonal_ = true; 700 for (Index i=0; i < NComps_Rows(); i++) { 701 for (Index j=0; j < NComps_Cols(); j++) { 702 if ( (i == j && IsNull(GetCompSpace(i,j))) 703 || (i != j && IsValid(GetCompSpace(i,j)))) { 704 diagonal_ = false; 705 break; 706 } 707 } 708 } 709 } 710 MakeNewCompoundMatrix() const711 CompoundMatrix* CompoundMatrixSpace::MakeNewCompoundMatrix() const 712 { 713 if (!dimensions_set_) { 714 dimensions_set_ = DimensionsSet(); 715 } 716 DBG_ASSERT(dimensions_set_); 717 718 CompoundMatrix* mat = new CompoundMatrix(this); 719 for (Index i=0; i<ncomps_rows_; i++) { 720 for (Index j=0; j<ncomps_cols_; j++) { 721 if (allocate_block_[i][j]) { 722 mat->SetCompNonConst(i, j, *GetCompSpace(i, j)->MakeNew()); 723 } 724 } 725 } 726 727 return mat; 728 } 729 DimensionsSet() const730 bool CompoundMatrixSpace::DimensionsSet() const 731 { 732 DBG_START_METH("CompoundMatrixSpace::DimensionsSet", 0); 733 Index total_nrows = 0; 734 Index total_ncols = 0; 735 bool valid = true; 736 for (Index i=0; i<ncomps_rows_; i++) { 737 if (block_rows_[i] == -1) { 738 valid = false; 739 break; 740 } 741 total_nrows += block_rows_[i]; 742 } 743 if (valid) { 744 for (Index j=0; j<ncomps_cols_; j++) { 745 if (block_cols_[j] == -1) { 746 valid = false; 747 break; 748 } 749 total_ncols += block_cols_[j]; 750 } 751 } 752 753 if (valid) { 754 DBG_ASSERT(total_nrows == NRows() && total_ncols == NCols()); 755 } 756 757 return valid; 758 } 759 760 } // namespace Ipopt 761