Changeset 533


Ignore:
Timestamp:
Oct 4, 2005 11:47:13 AM (14 years ago)
Author:
andreasw
Message:
  • Changed the conditions under which symbolic factorization is performed for Pardiso; based on discussion with Olaf Schenk
Location:
branches/dev/Algorithm/LinearSolvers
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/dev/Algorithm/LinearSolvers/IpPardisoSolverInterface.cpp

    r532 r533  
    102102    dim_=0;
    103103    nonzeros_=0;
    104     PHASE_=0;
     104    have_symbolic_factorization_=false;
    105105    initialized_=false;
    106106    delete[] a_;
     
    195195                   dbg_verbosity);
    196196
    197     IpData().TimingStats().LinearSystemSymbolicFactorization().Start();
    198 
    199     // Call Pardiso to do the analysis phase
    200     ipfint PHASE = 11;
    201     ipfint N = dim_;
    202     ipfint PERM;   // This should not be accessed by Pardiso
    203     ipfint NRHS = 0;
    204     double B;  // This should not be accessed by Pardiso in analysis
    205     // phase
    206     double X;  // This should not be accessed by Pardiso in analysis
    207     // phase
    208     ipfint ERROR;
    209 
    210     F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_,
    211                               &PHASE, &N, a_, ia, ja, &PERM,
    212                               &NRHS, IPARM_, &MSGLVL_, &B, &X,
    213                               &ERROR);
    214     if (ERROR!=0) {
    215       Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
    216                      "Error in Pardiso during analysis phase.  ERROR = %d.\n",
    217                      ERROR);
    218       IpData().TimingStats().LinearSystemSymbolicFactorization().End();
    219       return SYMSOLVER_FATAL_ERROR;
    220     }
    221 
    222     IpData().TimingStats().LinearSystemSymbolicFactorization().End();
     197    // Since Pardiso requires the values of the nonzeros of the matrix
     198    // for an efficient symbolic factorization, we postpone that task
     199    // until the first call of Factorize.  All we do here is to reset
     200    // the flag (in case this interface is called for a matrix with a
     201    // new structure).
     202
     203    have_symbolic_factorization_ = false;
    223204
    224205    return SYMSOLVER_SUCCESS;
     
    233214    DBG_START_METH("PardisoSolverInterface::Factorization",dbg_verbosity);
    234215
    235     IpData().TimingStats().LinearSystemFactorization().Start();
    236 
    237216    // Call Pardiso to do the factorization
    238217    ipfint N = dim_;
     218    ipfint PHASE;
    239219    ipfint PERM;   // This should not be accessed by Pardiso
    240220    ipfint NRHS = 0;
     
    245225    ipfint ERROR;
    246226
    247     // Dump matrix to file...
    248     ipfint  NNZ = ia[N]-1;
    249     ipfint   i;
    250 
    251     if (IpData().iter_count() != debug_last_iter_) {
    252       debug_cnt_ = 0;
    253     }
    254 
    255227    if (getenv ("IPOPT_WRITE_MAT")) {
     228
     229      // Dump matrix to file...
     230      ipfint  NNZ = ia[N]-1;
     231      ipfint   i;
     232
     233      if (IpData().iter_count() != debug_last_iter_) {
     234        debug_cnt_ = 0;
     235      }
     236
    256237      /* Write header */
    257238      FILE    *mat_file;
     
    287268
    288269      fclose (mat_file);
    289     }
    290 
    291     debug_last_iter_ = IpData().iter_count();
    292     debug_cnt_++;
    293 
    294     if (PHASE_ == 0)
    295       PHASE_ = 11;
    296     F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_,
    297                               &PHASE_, &N, a_, ia, ja, &PERM,
    298                               &NRHS, IPARM_, &MSGLVL_, &B, &X,
    299                               &ERROR);
    300     if (ERROR==-4) {
    301       // I think this means that the matrix is singular
     270
     271      debug_last_iter_ = IpData().iter_count();
     272      debug_cnt_++;
     273    }
     274
     275    bool done = false;
     276    bool just_performed_symbolic_factorization = false;
     277
     278    while (!done) {
     279      if (!have_symbolic_factorization_) {
     280        IpData().TimingStats().LinearSystemSymbolicFactorization().Start();
     281        PHASE = 11;
     282        F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_,
     283                                  &PHASE, &N, a_, ia, ja, &PERM,
     284                                  &NRHS, IPARM_, &MSGLVL_, &B, &X,
     285                                  &ERROR);
     286        IpData().TimingStats().LinearSystemSymbolicFactorization().End();
     287        if (ERROR!=0 ) {
     288          Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
     289                         "Error in Pardiso during symbolic factorization phase.  ERROR = %d.\n", ERROR);
     290          return SYMSOLVER_FATAL_ERROR;
     291        }
     292        have_symbolic_factorization_ = true;
     293        just_performed_symbolic_factorization = true;
     294      }
     295
     296      PHASE = 22;
     297
     298      IpData().TimingStats().LinearSystemFactorization().Start();
     299      F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_,
     300                                &PHASE, &N, a_, ia, ja, &PERM,
     301                                &NRHS, IPARM_, &MSGLVL_, &B, &X,
     302                                &ERROR);
    302303      IpData().TimingStats().LinearSystemFactorization().End();
    303       return SYMSOLVER_SINGULAR;
    304     }
    305     else if (ERROR!=0 ) {
    306       Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
    307                      "Error in Pardiso during factorization phase.  ERROR = %d.\n", ERROR);
    308       IpData().TimingStats().LinearSystemFactorization().End();
    309       return SYMSOLVER_FATAL_ERROR;
    310     }
    311     PHASE_ = 22;
    312 
    313     /* ToDo ask Olaf what this means. */
    314     if (IPARM_[13] != 0) {
    315       Jnlst().Printf(J_WARNING, J_LINEAR_ALGEBRA,
    316                      "Number of perturbed pivots in factorization phase = %d.\n", IPARM_[13]);
    317       PHASE_ = 12;
    318     }
    319 
    320     negevals_ = IPARM_[22];
     304
     305      if (ERROR==-4) {
     306        // I think this means that the matrix is singular
     307        // OLAF said that this will never happen (ToDo)
     308        return SYMSOLVER_SINGULAR;
     309      }
     310      else if (ERROR!=0 ) {
     311        Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
     312                       "Error in Pardiso during factorization phase.  ERROR = %d.\n", ERROR);
     313        return SYMSOLVER_FATAL_ERROR;
     314      }
     315
     316      negevals_ = IPARM_[22];
     317      if (IPARM_[13] != 0) {
     318        Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
     319                       "Number of perturbed pivots in factorization phase = %d.\n", IPARM_[13]);
     320        IpData().Append_info_string("Pp");
     321        // trigger a new symblic factorization for the next factorize
     322        // call, even if the current inertia estimate is correct
     323        // OLAF MICHAEL : Maybe we need to discuss this
     324        have_symbolic_factorization_ = false;
     325        // If we there were perturbed pivots, and the inertia of the
     326        // system is correct, and we haven't redone the symblic
     327        // factorization during this call yet, trigger a new
     328        // factorization after a symblic factorization
     329        done = (just_performed_symbolic_factorization || !check_NegEVals ||
     330                numberOfNegEVals==negevals_);
     331      }
     332      else {
     333        done = true;
     334      }
     335    }
    321336
    322337    DBG_ASSERT(IPARM_[21]+IPARM_[22] == dim_);
     
    328343                     "Wrong inertia: required are %d, but we got %d.\n",
    329344                     numberOfNegEVals, negevals_);
    330       IpData().TimingStats().LinearSystemFactorization().End();
    331345      return SYMSOLVER_WRONG_INERTIA;
    332346    }
    333347
    334     IpData().TimingStats().LinearSystemFactorization().End();
    335348    return SYMSOLVER_SUCCESS;
    336349  }
     
    358371
    359372    if (IPARM_[6] != 0) {
    360       Jnlst().Printf(J_WARNING, J_LINEAR_ALGEBRA,
     373      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
    361374                     "Number of iterative refinement steps = %d.\n", IPARM_[6]);
     375      IpData().Append_info_string("Pi");
    362376    }
    363377
  • branches/dev/Algorithm/LinearSolvers/IpPardisoSolverInterface.hpp

    r532 r533  
    124124     *  For initialization, this object needs to have seen a matrix */
    125125    bool initialized_;
     126    /** Flag indicating if symbolic factorization has already been
     127     *  performed. */
     128    bool have_symbolic_factorization_;
    126129    //@}
    127130
     
    141144    /** Message level. */
    142145    ipfint MSGLVL_;
    143     /** Flag indicating in what phase Pardiso is */
    144     ipfint PHASE_;
    145146    //@}
    146147
Note: See TracChangeset for help on using the changeset viewer.