1 // Copyright (C) 2008, 2011 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: Andreas Waechter IBM 2008-09-19 8 // 9 // based on IpPardisoSolverInterface.cpp rev 1219 10 11 #include "IpoptConfig.h" 12 #include "IpIterativePardisoSolverInterface.hpp" 13 #include "IpBlas.hpp" 14 # include <math.h> 15 16 #ifdef HAVE_CSTDIO 17 # include <cstdio> 18 #else 19 # ifdef HAVE_STDIO_H 20 # include <stdio.h> 21 # else 22 # error "don't have header file for stdio" 23 # endif 24 #endif 25 26 #ifdef HAVE_CSTDLIB 27 # include <cstdlib> 28 #else 29 # ifdef HAVE_STDLIB_H 30 # include <stdlib.h> 31 # else 32 # error "don't have header file for stdlib" 33 # endif 34 #endif 35 36 #ifdef HAVE_CSTRING 37 # include <cstring> 38 #else 39 # ifdef HAVE_STRING_H 40 # include <string.h> 41 # else 42 # error "don't have header file for string" 43 # endif 44 #endif 45 46 47 Ipopt::IterativeSolverTerminationTester* global_tester_ptr_; 48 Ipopt::IterativeSolverTerminationTester::ETerminationTest test_result_; 49 extern "C" 50 { IpoptTerminationTest(int n,double * sol,double * resid,int iter,double norm2_rhs)51 int IpoptTerminationTest(int n, double* sol, double* resid, int iter, double norm2_rhs) { 52 fflush(stdout); 53 fflush(stderr); 54 // only do the termination test for PD system 55 test_result_ = global_tester_ptr_->TestTermination(n, sol, resid, iter, norm2_rhs); 56 global_tester_ptr_->GetJnlst().Printf(Ipopt::J_DETAILED, Ipopt::J_LINEAR_ALGEBRA, 57 "Termination Tester Result = %d.\n", 58 test_result_); 59 switch (test_result_) { 60 case Ipopt::IterativeSolverTerminationTester::CONTINUE: 61 return false; 62 break; 63 default: 64 return true; 65 break; 66 } 67 } 68 69 // The following global function pointer is defined in the Pardiso library 70 void SetIpoptCallbackFunction(int (*IpoptFunction)(int n, double* x, double* r, int k, double b)); 71 } 72 73 /** Prototypes for Pardiso's subroutines */ 74 extern "C" 75 { 76 void F77_FUNC(pardisoinit,PARDISOINIT)(void* PT, const ipfint* MTYPE, 77 const ipfint* SOLVER, 78 ipfint* IPARM, 79 double* DPARM, 80 ipfint* ERROR); 81 void F77_FUNC(pardiso,PARDISO)(void** PT, const ipfint* MAXFCT, 82 const ipfint* MNUM, const ipfint* MTYPE, 83 const ipfint* PHASE, const ipfint* N, 84 const double* A, const ipfint* IA, 85 const ipfint* JA, const ipfint* PERM, 86 const ipfint* NRHS, ipfint* IPARM, 87 const ipfint* MSGLVL, double* B, double* X, 88 ipfint* ERROR, double* DPARM); 89 } 90 91 namespace Ipopt 92 { 93 #if COIN_IPOPT_VERBOSITY > 0 94 static const Index dbg_verbosity = 0; 95 #endif 96 97 IterativePardisoSolverInterface:: IterativePardisoSolverInterface(IterativeSolverTerminationTester & normal_tester,IterativeSolverTerminationTester & pd_tester)98 IterativePardisoSolverInterface(IterativeSolverTerminationTester& normal_tester, 99 IterativeSolverTerminationTester& pd_tester) 100 : 101 a_(NULL), 102 negevals_(-1), 103 initialized_(false), 104 105 MAXFCT_(1), 106 MNUM_(1), 107 MTYPE_(-2), 108 MSGLVL_(0), 109 debug_last_iter_(-1), 110 normal_tester_(&normal_tester), 111 pd_tester_(&pd_tester) 112 { 113 DBG_START_METH("IterativePardisoSolverInterface::IterativePardisoSolverInterface()",dbg_verbosity); 114 115 PT_ = new void*[64]; 116 IPARM_ = new ipfint[64]; 117 DPARM_ = new double[64]; 118 } 119 ~IterativePardisoSolverInterface()120 IterativePardisoSolverInterface::~IterativePardisoSolverInterface() 121 { 122 DBG_START_METH("IterativePardisoSolverInterface::~IterativePardisoSolverInterface()", 123 dbg_verbosity); 124 125 // Tell Pardiso to release all memory 126 if (initialized_) { 127 ipfint PHASE = -1; 128 ipfint N = dim_; 129 ipfint NRHS = 0; 130 ipfint ERROR; 131 ipfint idmy; 132 double ddmy; 133 F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_, &PHASE, &N, 134 &ddmy, &idmy, &idmy, &idmy, &NRHS, IPARM_, 135 &MSGLVL_, &ddmy, &ddmy, &ERROR, DPARM_); 136 DBG_ASSERT(ERROR==0); 137 } 138 139 delete[] PT_; 140 delete[] IPARM_; 141 delete[] DPARM_; 142 delete[] a_; 143 } 144 RegisterOptions(SmartPtr<RegisteredOptions> roptions)145 void IterativePardisoSolverInterface::RegisterOptions(SmartPtr<RegisteredOptions> roptions) 146 {} 147 InitializeImpl(const OptionsList & options,const std::string & prefix)148 bool IterativePardisoSolverInterface::InitializeImpl(const OptionsList& options, 149 const std::string& prefix) 150 { 151 #ifdef HAVE_PARDISO_OLDINTERFACE 152 THROW_EXCEPTION(OPTION_INVALID, "The inexact version works only with a new version of Pardiso (at least 4.0)"); 153 #endif 154 155 Index enum_int; 156 options.GetEnumValue("pardiso_matching_strategy", enum_int, prefix); 157 match_strat_ = PardisoMatchingStrategy(enum_int); 158 options.GetBoolValue("pardiso_redo_symbolic_fact_only_if_inertia_wrong", 159 pardiso_redo_symbolic_fact_only_if_inertia_wrong_, 160 prefix); 161 options.GetBoolValue("pardiso_repeated_perturbation_means_singular", 162 pardiso_repeated_perturbation_means_singular_, 163 prefix); 164 //Index pardiso_out_of_core_power; 165 //options.GetIntegerValue("pardiso_out_of_core_power", 166 // pardiso_out_of_core_power, prefix); 167 options.GetBoolValue("pardiso_skip_inertia_check", 168 skip_inertia_check_, prefix); 169 int max_iterref_steps; 170 options.GetIntegerValue("pardiso_max_iterative_refinement_steps", max_iterref_steps, prefix); 171 172 // PD system 173 options.GetIntegerValue("pardiso_max_iter", pardiso_max_iter_, prefix); 174 options.GetNumericValue("pardiso_iter_relative_tol", 175 pardiso_iter_relative_tol_, prefix); 176 options.GetIntegerValue("pardiso_iter_coarse_size", 177 pardiso_iter_coarse_size_, prefix); 178 options.GetIntegerValue("pardiso_iter_max_levels", 179 pardiso_iter_max_levels_, prefix); 180 options.GetNumericValue("pardiso_iter_dropping_factor", 181 pardiso_iter_dropping_factor_, prefix); 182 options.GetNumericValue("pardiso_iter_dropping_schur", 183 pardiso_iter_dropping_schur_, prefix); 184 options.GetIntegerValue("pardiso_iter_max_row_fill", 185 pardiso_iter_max_row_fill_, prefix); 186 options.GetNumericValue("pardiso_iter_inverse_norm_factor", 187 pardiso_iter_inverse_norm_factor_, prefix); 188 // Normal system 189 options.GetIntegerValue("pardiso_max_iter", normal_pardiso_max_iter_, 190 prefix+"normal."); 191 options.GetNumericValue("pardiso_iter_relative_tol", 192 normal_pardiso_iter_relative_tol_, 193 prefix+"normal."); 194 options.GetIntegerValue("pardiso_iter_coarse_size", 195 normal_pardiso_iter_coarse_size_, 196 prefix+"normal."); 197 options.GetIntegerValue("pardiso_iter_max_levels", 198 normal_pardiso_iter_max_levels_, 199 prefix+"normal."); 200 options.GetNumericValue("pardiso_iter_dropping_factor", 201 normal_pardiso_iter_dropping_factor_, 202 prefix+"normal."); 203 options.GetNumericValue("pardiso_iter_dropping_schur", 204 normal_pardiso_iter_dropping_schur_, 205 prefix+"normal."); 206 options.GetIntegerValue("pardiso_iter_max_row_fill", 207 normal_pardiso_iter_max_row_fill_, 208 prefix+"normal."); 209 options.GetNumericValue("pardiso_iter_inverse_norm_factor", 210 normal_pardiso_iter_inverse_norm_factor_, 211 prefix+"normal."); 212 213 int pardiso_msglvl; 214 options.GetIntegerValue("pardiso_msglvl", pardiso_msglvl, prefix); 215 int order; 216 options.GetEnumValue("pardiso_order", order, prefix); 217 options.GetIntegerValue("pardiso_max_droptol_corrections", 218 pardiso_max_droptol_corrections_, prefix); 219 220 // Number value = 0.0; 221 222 // Tell Pardiso to release all memory if it had been used before 223 if (initialized_) { 224 ipfint PHASE = -1; 225 ipfint N = dim_; 226 ipfint NRHS = 0; 227 ipfint ERROR; 228 ipfint idmy; 229 double ddmy; 230 F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_, &PHASE, &N, 231 &ddmy, &idmy, &idmy, &idmy, &NRHS, IPARM_, 232 &MSGLVL_, &ddmy, &ddmy, &ERROR, DPARM_) ; 233 DBG_ASSERT(ERROR==0); 234 } 235 236 // Reset all private data 237 dim_=0; 238 nonzeros_=0; 239 have_symbolic_factorization_=false; 240 initialized_=false; 241 delete[] a_; 242 a_ = NULL; 243 244 // Call Pardiso's initialization routine 245 IPARM_[0] = 0; // Tell it to fill IPARM with default values(?) 246 ipfint ERROR = 0; 247 ipfint SOLVER = 1; // initialze only direct solver 248 F77_FUNC(pardisoinit,PARDISOINIT)(PT_, &MTYPE_, &SOLVER, 249 IPARM_, DPARM_, &ERROR); 250 251 // Set some parameters for Pardiso 252 IPARM_[0] = 1; // Don't use the default values 253 254 #if defined(HAVE_PARDISO_PARALLEL) || ! defined(HAVE_PARDISO) 255 // Obtain the numbers of processors from the value of OMP_NUM_THREADS 256 char *var = getenv("OMP_NUM_THREADS"); 257 int num_procs = 1; 258 if (var != NULL) { 259 sscanf( var, "%d", &num_procs ); 260 if (num_procs < 1) { 261 Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA, 262 "Invalid value for OMP_NUM_THREADS (\"%s\").\n", var); 263 return false; 264 } 265 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 266 "Using environment OMP_NUM_THREADS = %d as the number of processors.\n", num_procs); 267 } 268 #ifdef HAVE_PARDISO 269 // If we run Pardiso through the linear solver loader, 270 // we do not know whether it is the parallel version, so we do not report an error if OMP_NUM_THREADS is not set. 271 else { 272 Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA, 273 "You need to set environment variable OMP_NUM_THREADS to the number of processors used in Pardiso (e.g., 1).\n\n"); 274 return false; 275 } 276 #endif 277 IPARM_[2] = num_procs; // Set the number of processors 278 #else 279 280 IPARM_[2] = 1; 281 #endif 282 283 IPARM_[1] = order; 284 IPARM_[5] = 1; // Overwrite right-hand side 285 IPARM_[7] = max_iterref_steps; 286 287 // Options suggested by Olaf Schenk 288 IPARM_[9] = 12; 289 IPARM_[10] = 2; // Results in better scaling 290 // Matching information: IPARM_[12] = 1 seems ok, but results in a 291 // large number of pivot perturbation 292 // Matching information: IPARM_[12] = 2 robust, but more expensive method 293 IPARM_[12] = (int)match_strat_; 294 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 295 "Pardiso matching strategy (IPARM(13)): %d\n", IPARM_[12]); 296 297 IPARM_[20] = 3; // Results in better accuracy 298 IPARM_[23] = 1; // parallel fac 299 IPARM_[24] = 1; // parallel solve 300 IPARM_[28] = 0; // 32-bit factorization 301 IPARM_[29] = 1; // we need this for IPOPT interface 302 303 IPARM_[31] = 1 ; // iterative solver 304 MSGLVL_ = pardiso_msglvl; 305 306 pardiso_iter_dropping_factor_used_ = pardiso_iter_dropping_factor_; 307 pardiso_iter_dropping_schur_used_ = pardiso_iter_dropping_schur_; 308 normal_pardiso_iter_dropping_factor_used_ = normal_pardiso_iter_dropping_factor_; 309 normal_pardiso_iter_dropping_schur_used_ = normal_pardiso_iter_dropping_schur_; 310 311 // TODO Make option 312 decr_factor_ = 1./3.; 313 314 // Option for the out of core variant 315 // IPARM_[49] = pardiso_out_of_core_power; 316 317 SetIpoptCallbackFunction(&IpoptTerminationTest); 318 319 bool retval = normal_tester_->Initialize(Jnlst(), IpNLP(), IpData(), 320 IpCq(), options, prefix); 321 if (retval) { 322 retval = pd_tester_->Initialize(Jnlst(), IpNLP(), IpData(), 323 IpCq(), options, prefix); 324 } 325 326 return retval; 327 } 328 MultiSolve(bool new_matrix,const Index * ia,const Index * ja,Index nrhs,double * rhs_vals,bool check_NegEVals,Index numberOfNegEVals)329 ESymSolverStatus IterativePardisoSolverInterface::MultiSolve(bool new_matrix, 330 const Index* ia, 331 const Index* ja, 332 Index nrhs, 333 double* rhs_vals, 334 bool check_NegEVals, 335 Index numberOfNegEVals) 336 { 337 DBG_START_METH("IterativePardisoSolverInterface::MultiSolve",dbg_verbosity); 338 DBG_ASSERT(!check_NegEVals || ProvidesInertia()); 339 DBG_ASSERT(initialized_); 340 341 // check if a factorization has to be done 342 if (new_matrix) { 343 // perform the factorization 344 ESymSolverStatus retval; 345 retval = Factorization(ia, ja, check_NegEVals, numberOfNegEVals); 346 if (retval!=SYMSOLVER_SUCCESS) { 347 DBG_PRINT((1, "FACTORIZATION FAILED!\n")); 348 return retval; // Matrix singular or error occurred 349 } 350 } 351 352 // do the solve 353 return Solve(ia, ja, nrhs, rhs_vals); 354 } 355 GetValuesArrayPtr()356 double* IterativePardisoSolverInterface::GetValuesArrayPtr() 357 { 358 DBG_ASSERT(initialized_); 359 DBG_ASSERT(a_); 360 return a_; 361 } 362 363 /** Initialize the local copy of the positions of the nonzero 364 elements */ InitializeStructure(Index dim,Index nonzeros,const Index * ia,const Index * ja)365 ESymSolverStatus IterativePardisoSolverInterface::InitializeStructure 366 (Index dim, Index nonzeros, 367 const Index* ia, 368 const Index* ja) 369 { 370 DBG_START_METH("IterativePardisoSolverInterface::InitializeStructure",dbg_verbosity); 371 dim_ = dim; 372 nonzeros_ = nonzeros; 373 374 // Make space for storing the matrix elements 375 delete[] a_; 376 a_ = NULL; 377 a_ = new double[nonzeros_]; 378 379 // Do the symbolic facotrization 380 ESymSolverStatus retval = SymbolicFactorization(ia, ja); 381 if (retval != SYMSOLVER_SUCCESS) { 382 return retval; 383 } 384 385 initialized_ = true; 386 387 return retval; 388 } 389 390 ESymSolverStatus SymbolicFactorization(const Index * ia,const Index * ja)391 IterativePardisoSolverInterface::SymbolicFactorization(const Index* ia, 392 const Index* ja) 393 { 394 DBG_START_METH("IterativePardisoSolverInterface::SymbolicFactorization", 395 dbg_verbosity); 396 397 // Since Pardiso requires the values of the nonzeros of the matrix 398 // for an efficient symbolic factorization, we postpone that task 399 // until the first call of Factorize. All we do here is to reset 400 // the flag (in case this interface is called for a matrix with a 401 // new structure). 402 403 have_symbolic_factorization_ = false; 404 405 return SYMSOLVER_SUCCESS; 406 } 407 408 static void write_iajaa_matrix(int N,const Index * ia,const Index * ja,double * a_,double * rhs_vals,int iter_cnt,int sol_cnt)409 write_iajaa_matrix (int N, 410 const Index* ia, 411 const Index* ja, 412 double* a_, 413 double* rhs_vals, 414 int iter_cnt, 415 int sol_cnt) 416 { 417 if (getenv ("IPOPT_WRITE_MAT")) { 418 /* Write header */ 419 FILE *mat_file; 420 char mat_name[128]; 421 char mat_pref[32]; 422 423 ipfint NNZ = ia[N]-1; 424 ipfint i; 425 426 if (getenv ("IPOPT_WRITE_PREFIX")) 427 strcpy (mat_pref, getenv ("IPOPT_WRITE_PREFIX")); 428 else 429 strcpy (mat_pref, "mat-ipopt"); 430 431 Snprintf (mat_name, 127, "%s_%03d-%02d.iajaa", 432 mat_pref, iter_cnt, sol_cnt); 433 434 // Open and write matrix file. 435 mat_file = fopen (mat_name, "w"); 436 437 fprintf (mat_file, "%d\n", N); 438 fprintf (mat_file, "%d\n", NNZ); 439 440 for (i = 0; i < N+1; i++) 441 fprintf (mat_file, "%d\n", ia[i]); 442 for (i = 0; i < NNZ; i++) 443 fprintf (mat_file, "%d\n", ja[i]); 444 for (i = 0; i < NNZ; i++) 445 fprintf (mat_file, "%32.24e\n", a_[i]); 446 447 /* Right hand side. */ 448 if (rhs_vals) 449 for (i = 0; i < N; i++) 450 //FIXME: PUT BACK ORIGINAL: fprintf (mat_file, "%32.24e\n", rhs_vals[i]); 451 fprintf (mat_file, "%32.24e\n", -rhs_vals[i]); 452 453 fclose (mat_file); 454 } 455 } 456 457 ESymSolverStatus Factorization(const Index * ia,const Index * ja,bool check_NegEVals,Index numberOfNegEVals)458 IterativePardisoSolverInterface::Factorization(const Index* ia, 459 const Index* ja, 460 bool check_NegEVals, 461 Index numberOfNegEVals) 462 { 463 DBG_START_METH("IterativePardisoSolverInterface::Factorization",dbg_verbosity); 464 465 // Call Pardiso to do the factorization 466 ipfint PHASE ; 467 ipfint N = dim_; 468 ipfint PERM; // This should not be accessed by Pardiso 469 ipfint NRHS = 0; 470 double B; // This should not be accessed by Pardiso in factorization 471 // phase 472 double X; // This should not be accessed by Pardiso in factorization 473 // phase 474 ipfint ERROR; 475 476 bool done = false; 477 bool just_performed_symbolic_factorization = false; 478 479 while (!done) { 480 bool is_normal = false; 481 if (IsNull(InexData().normal_x()) && InexData().compute_normal()) { 482 is_normal = true; 483 } 484 if (!have_symbolic_factorization_) { 485 if (HaveIpData()) { 486 IpData().TimingStats().LinearSystemSymbolicFactorization().Start(); 487 } 488 PHASE = 11; 489 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 490 "Calling Pardiso for symbolic factorization.\n"); 491 492 if (is_normal) { 493 DPARM_[ 0] = normal_pardiso_max_iter_; 494 DPARM_[ 1] = normal_pardiso_iter_relative_tol_; 495 DPARM_[ 2] = normal_pardiso_iter_coarse_size_; 496 DPARM_[ 3] = normal_pardiso_iter_max_levels_; 497 DPARM_[ 4] = normal_pardiso_iter_dropping_factor_used_; 498 DPARM_[ 5] = normal_pardiso_iter_dropping_schur_used_; 499 DPARM_[ 6] = normal_pardiso_iter_max_row_fill_; 500 DPARM_[ 7] = normal_pardiso_iter_inverse_norm_factor_; 501 } 502 else { 503 DPARM_[ 0] = pardiso_max_iter_; 504 DPARM_[ 1] = pardiso_iter_relative_tol_; 505 DPARM_[ 2] = pardiso_iter_coarse_size_; 506 DPARM_[ 3] = pardiso_iter_max_levels_; 507 DPARM_[ 4] = pardiso_iter_dropping_factor_used_; 508 DPARM_[ 5] = pardiso_iter_dropping_schur_used_; 509 DPARM_[ 6] = pardiso_iter_max_row_fill_; 510 DPARM_[ 7] = pardiso_iter_inverse_norm_factor_; 511 } 512 DPARM_[ 8] = 25; // maximum number of non-improvement steps 513 514 F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_, 515 &PHASE, &N, a_, ia, ja, &PERM, 516 &NRHS, IPARM_, &MSGLVL_, &B, &X, 517 &ERROR, DPARM_); 518 if (HaveIpData()) { 519 IpData().TimingStats().LinearSystemSymbolicFactorization().End(); 520 } 521 if (ERROR==-7) { 522 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 523 "Pardiso symbolic factorization returns ERROR = %d. Matrix is singular.\n", ERROR); 524 return SYMSOLVER_SINGULAR; 525 } 526 else if (ERROR!=0) { 527 Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA, 528 "Error in Pardiso during symbolic factorization phase. ERROR = %d.\n", ERROR); 529 return SYMSOLVER_FATAL_ERROR; 530 } 531 have_symbolic_factorization_ = true; 532 just_performed_symbolic_factorization = true; 533 534 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 535 "Memory in KB required for the symbolic factorization = %d.\n", IPARM_[14]); 536 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 537 "Integer memory in KB required for the numerical factorization = %d.\n", IPARM_[15]); 538 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 539 "Double memory in KB required for the numerical factorization = %d.\n", IPARM_[16]); 540 } 541 542 PHASE = 22; 543 544 if (HaveIpData()) { 545 IpData().TimingStats().LinearSystemFactorization().Start(); 546 } 547 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 548 "Calling Pardiso for factorization.\n"); 549 // Dump matrix to file, and count number of solution steps. 550 if (HaveIpData()) { 551 if (IpData().iter_count() != debug_last_iter_) 552 debug_cnt_ = 0; 553 debug_last_iter_ = IpData().iter_count(); 554 debug_cnt_ ++; 555 } 556 else { 557 debug_cnt_ = 0; 558 debug_last_iter_ = 0; 559 } 560 561 if (is_normal) { 562 DPARM_[ 0] = normal_pardiso_max_iter_; 563 DPARM_[ 1] = normal_pardiso_iter_relative_tol_; 564 DPARM_[ 2] = normal_pardiso_iter_coarse_size_; 565 DPARM_[ 3] = normal_pardiso_iter_max_levels_; 566 DPARM_[ 4] = normal_pardiso_iter_dropping_factor_used_; 567 DPARM_[ 5] = normal_pardiso_iter_dropping_schur_used_; 568 DPARM_[ 6] = normal_pardiso_iter_max_row_fill_; 569 DPARM_[ 7] = normal_pardiso_iter_inverse_norm_factor_; 570 } 571 else { 572 DPARM_[ 0] = pardiso_max_iter_; 573 DPARM_[ 1] = pardiso_iter_relative_tol_; 574 DPARM_[ 2] = pardiso_iter_coarse_size_; 575 DPARM_[ 3] = pardiso_iter_max_levels_; 576 DPARM_[ 4] = pardiso_iter_dropping_factor_used_; 577 DPARM_[ 5] = pardiso_iter_dropping_schur_used_; 578 DPARM_[ 6] = pardiso_iter_max_row_fill_; 579 DPARM_[ 7] = pardiso_iter_inverse_norm_factor_; 580 } 581 DPARM_[ 8] = 25; // maximum number of non-improvement steps 582 583 F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_, 584 &PHASE, &N, a_, ia, ja, &PERM, 585 &NRHS, IPARM_, &MSGLVL_, &B, &X, 586 &ERROR, DPARM_); 587 if (HaveIpData()) { 588 IpData().TimingStats().LinearSystemFactorization().End(); 589 } 590 591 if (ERROR==-7) { 592 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 593 "Pardiso factorization returns ERROR = %d. Matrix is singular.\n", ERROR); 594 return SYMSOLVER_SINGULAR; 595 } 596 else if (ERROR==-4) { 597 // I think this means that the matrix is singular 598 // OLAF said that this will never happen (ToDo) 599 return SYMSOLVER_SINGULAR; 600 } 601 else if (ERROR!=0 ) { 602 Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA, 603 "Error in Pardiso during factorization phase. ERROR = %d.\n", ERROR); 604 return SYMSOLVER_FATAL_ERROR; 605 } 606 607 negevals_ = Max(IPARM_[22], numberOfNegEVals); 608 if (IPARM_[13] != 0) { 609 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 610 "Number of perturbed pivots in factorization phase = %d.\n", IPARM_[13]); 611 if (HaveIpData()) { 612 IpData().Append_info_string("Pp"); 613 } 614 done = true; 615 } 616 else { 617 done = true; 618 } 619 } 620 621 // DBG_ASSERT(IPARM_[21]+IPARM_[22] == dim_); 622 623 // Check whether the number of negative eigenvalues matches the requested 624 // count 625 if (skip_inertia_check_) numberOfNegEVals=negevals_; 626 627 if (check_NegEVals && (numberOfNegEVals!=negevals_)) { 628 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 629 "Wrong inertia: required are %d, but we got %d.\n", 630 numberOfNegEVals, negevals_); 631 return SYMSOLVER_WRONG_INERTIA; 632 } 633 634 return SYMSOLVER_SUCCESS; 635 } 636 Solve(const Index * ia,const Index * ja,Index nrhs,double * rhs_vals)637 ESymSolverStatus IterativePardisoSolverInterface::Solve(const Index* ia, 638 const Index* ja, 639 Index nrhs, 640 double *rhs_vals) 641 { 642 DBG_START_METH("IterativePardisoSolverInterface::Solve",dbg_verbosity); 643 644 DBG_ASSERT(nrhs==1); 645 646 if (HaveIpData()) { 647 IpData().TimingStats().LinearSystemBackSolve().Start(); 648 } 649 // Call Pardiso to do the solve for the given right-hand sides 650 ipfint PHASE = 33; 651 ipfint N = dim_; 652 ipfint PERM; // This should not be accessed by Pardiso 653 ipfint NRHS = nrhs; 654 double* X = new double[nrhs*dim_]; 655 double* ORIG_RHS = new double[nrhs*dim_]; 656 ipfint ERROR; 657 658 // Initialize solution with zero and save right hand side 659 for (int i = 0; i < N; i++) { 660 X[i] = 0; 661 ORIG_RHS[i] = rhs_vals[i]; 662 } 663 664 // Dump matrix to file if requested 665 Index iter_count = 0; 666 if (HaveIpData()) { 667 iter_count = IpData().iter_count(); 668 } 669 write_iajaa_matrix (N, ia, ja, a_, rhs_vals, iter_count, debug_cnt_); 670 671 IterativeSolverTerminationTester* tester; 672 673 int attempts = 0; 674 const int max_attempts = pardiso_max_droptol_corrections_+1; 675 676 bool is_normal = false; 677 if (IsNull(InexData().normal_x()) && InexData().compute_normal()) { 678 tester = GetRawPtr(normal_tester_); 679 is_normal = true; 680 } 681 else { 682 tester = GetRawPtr(pd_tester_); 683 } 684 global_tester_ptr_ = tester; 685 686 while (attempts<max_attempts) { 687 bool retval = tester->InitializeSolve(); 688 ASSERT_EXCEPTION(retval, INTERNAL_ABORT, "tester->InitializeSolve(); returned false"); 689 690 for (int i = 0; i < N; i++) { 691 rhs_vals[i] = ORIG_RHS[i]; 692 } 693 694 DPARM_[ 8] = 25; // non_improvement in SQMR iteration 695 F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_, 696 &PHASE, &N, a_, ia, ja, &PERM, 697 &NRHS, IPARM_, &MSGLVL_, rhs_vals, X, 698 &ERROR, DPARM_); 699 700 if (ERROR <= -100 && ERROR >= -110) { 701 Jnlst().Printf(J_WARNING, J_LINEAR_ALGEBRA, 702 "Iterative solver in Pardiso did not converge (ERROR = %d)\n", ERROR); 703 Jnlst().Printf(J_WARNING, J_LINEAR_ALGEBRA, 704 " Decreasing drop tolerances from DPARM_[ 4] = %e and DPARM_[ 5] = %e ", DPARM_[ 4], DPARM_[ 5]); 705 if (is_normal) { 706 Jnlst().Printf(J_WARNING, J_LINEAR_ALGEBRA, 707 "(normal step)\n"); 708 } 709 else { 710 Jnlst().Printf(J_WARNING, J_LINEAR_ALGEBRA, 711 "(PD step)\n"); 712 } 713 PHASE = 23; 714 DPARM_[ 4] *= decr_factor_; 715 DPARM_[ 5] *= decr_factor_; 716 Jnlst().Printf(J_WARNING, J_LINEAR_ALGEBRA, 717 " to DPARM_[ 4] = %e and DPARM_[ 5] = %e\n", DPARM_[ 4], DPARM_[ 5]); 718 // Copy solution back to y to get intial values for the next iteration 719 attempts++; 720 ERROR = 0; 721 } 722 else { 723 attempts = max_attempts; 724 Index iterations_used = tester->GetSolverIterations(); 725 if (is_normal) { 726 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 727 "Number of iterations in Pardiso iterative solver for normal step = %d.\n", iterations_used); 728 } 729 else { 730 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 731 "Number of iterations in Pardiso iterative solver for PD step = %d.\n", iterations_used); 732 } 733 } 734 tester->Clear(); 735 } 736 737 if (is_normal) { 738 if (DPARM_[4] < normal_pardiso_iter_dropping_factor_) { 739 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 740 "Increasing drop tolerances from DPARM_[ 4] = %e and DPARM_[ 5] = %e (normal step\n", DPARM_[ 4], DPARM_[ 5]); 741 } 742 normal_pardiso_iter_dropping_factor_used_ = 743 Min(DPARM_[4]/decr_factor_, normal_pardiso_iter_dropping_factor_); 744 normal_pardiso_iter_dropping_schur_used_ = 745 Min(DPARM_[5]/decr_factor_, normal_pardiso_iter_dropping_schur_); 746 if (DPARM_[4] < normal_pardiso_iter_dropping_factor_) { 747 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 748 " to DPARM_[ 4] = %e and DPARM_[ 5] = %e for next iteration.\n", normal_pardiso_iter_dropping_factor_used_, normal_pardiso_iter_dropping_schur_used_); 749 } 750 } 751 else { 752 if (DPARM_[4] < pardiso_iter_dropping_factor_) { 753 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 754 "Increasing drop tolerances from DPARM_[ 4] = %e and DPARM_[ 5] = %e (PD step\n", DPARM_[ 4], DPARM_[ 5]); 755 } 756 pardiso_iter_dropping_factor_used_ = 757 Min(DPARM_[4]/decr_factor_, pardiso_iter_dropping_factor_); 758 pardiso_iter_dropping_schur_used_ = 759 Min(DPARM_[5]/decr_factor_, pardiso_iter_dropping_schur_); 760 if (DPARM_[4] < pardiso_iter_dropping_factor_) { 761 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 762 " to DPARM_[ 4] = %e and DPARM_[ 5] = %e for next iteration.\n", pardiso_iter_dropping_factor_used_, pardiso_iter_dropping_schur_used_); 763 } 764 } 765 766 delete [] X; 767 delete [] ORIG_RHS; 768 769 if (IPARM_[6] != 0) { 770 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 771 "Number of iterative refinement steps = %d.\n", IPARM_[6]); 772 if (HaveIpData()) { 773 IpData().Append_info_string("Pi"); 774 } 775 } 776 777 if (HaveIpData()) { 778 IpData().TimingStats().LinearSystemBackSolve().End(); 779 } 780 if (ERROR!=0 ) { 781 Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA, 782 "Error in Pardiso during solve phase. ERROR = %d.\n", ERROR); 783 return SYMSOLVER_FATAL_ERROR; 784 } 785 if (test_result_ == IterativeSolverTerminationTester::MODIFY_HESSIAN) { 786 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 787 "Termination tester requests modification of Hessian\n"); 788 return SYMSOLVER_WRONG_INERTIA; 789 } 790 #if 0 791 // FRANK: look at this: 792 if (test_result_ == IterativeSolverTerminationTester::CONTINUE) { 793 if (InexData().compute_normal()) { 794 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 795 "Termination tester not satisfied!!! Pretend singular\n"); 796 return SYMSOLVER_SINGULAR; 797 } 798 } 799 #endif 800 if (test_result_ == IterativeSolverTerminationTester::TEST_2_SATISFIED) { 801 // Termination Test 2 is satisfied, set the step for the primal 802 // iterates to zero 803 Index nvars = IpData().curr()->x()->Dim() + IpData().curr()->s()->Dim(); 804 const Number zero = 0.; 805 IpBlasDcopy(nvars, &zero, 0, rhs_vals, 1); 806 } 807 return SYMSOLVER_SUCCESS; 808 } 809 NumberOfNegEVals() const810 Index IterativePardisoSolverInterface::NumberOfNegEVals() const 811 { 812 DBG_START_METH("IterativePardisoSolverInterface::NumberOfNegEVals",dbg_verbosity); 813 DBG_ASSERT(negevals_>=0); 814 return negevals_; 815 } 816 IncreaseQuality()817 bool IterativePardisoSolverInterface::IncreaseQuality() 818 { 819 // At the moment, I don't see how we could tell Pardiso to do better 820 // (maybe switch from IPARM[20]=1 to IPARM[20]=2?) 821 return false; 822 } 823 824 } // namespace Ipopt 825