Changeset 531


Ignore:
Timestamp:
Oct 3, 2005 4:06:26 PM (14 years ago)
Author:
andreasw
Message:

Included Olaf Schenk's modifications (choice of Pardiso options, performing pivoting analysis with
values)

Location:
branches/dev/Algorithm/LinearSolvers
Files:
2 edited

Legend:

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

    r501 r531  
    66//
    77// Authors:  Carl Laird, Andreas Waechter     IBM    2005-03-17
     8//
     9//           Olaf Schenk                      Univ of Basel 2005-09-20
     10//                  - changed options, added PHASE_ flag
    811
    912#include "IpPardisoSolverInterface.hpp"
     13
    1014
    1115#ifdef HAVE_CMATH
     
    4246  PardisoSolverInterface::PardisoSolverInterface()
    4347      :
    44       dim_(0),
    45       nonzeros_(0),
     48      a_(NULL),
    4649      initialized_(false),
    47       negevals_(-1),
    4850
    4951      MAXFCT_(1),
    5052      MNUM_(1),
    5153      MTYPE_(-2),
    52       MSGLVL_(0),
    53 
    54       a_(NULL)
     54      MSGLVL_(0)
    5555  {
    5656    DBG_START_METH("PardisoSolverInterface::PardisoSolverInterface()",dbg_verbosity);
     
    8585      const std::string& prefix)
    8686  {
    87     Number value = 0.0;
     87    // Number value = 0.0;
    8888
    8989    // Tell Pardiso to release all memory if it had been used before
     
    102102    dim_=0;
    103103    nonzeros_=0;
     104    PHASE_=0;
    104105    initialized_=false;
    105106    delete[] a_;
     
    113114    IPARM_[2] = 1;  // Only one CPU for now
    114115    IPARM_[5] = 1;  // Overwrite right-hand side
    115 
    116116    // ToDo: decide if we need iterative refinement in Pardiso.  For
    117117    // now, switch it off ?
    118118    IPARM_[7] = 0;
    119119
    120     // IPARM_[20] = 2;
     120    // Options suggested by Olaf Schenk
     121    IPARM_[9] = 12;
     122    IPARM_[10] = 1;
     123    IPARM_[12] = 1;
     124
     125    IPARM_[20] = 1;
    121126
    122127    return true;
     
    224229
    225230    // Call Pardiso to do the factorization
    226     ipfint PHASE = 22;
    227231    ipfint N = dim_;
    228232    ipfint PERM;   // This should not be accessed by Pardiso
     
    234238    ipfint ERROR;
    235239
     240    // Dump matrix to file...
     241    ipfint  NNZ = ia[N]-1;
     242    ipfint   i;
     243
     244    if (IpData().iter_count() != debug_last_iter) {
     245      debug_cnt = 0;
     246    }
     247
     248    if (getenv ("IPOPT_WRITE_MAT")) {
     249      /* Write header */
     250      FILE    *mat_file;
     251      char     mat_name[128];
     252      char     mat_pref[32];
     253
     254      if (getenv ("IPOPT_WRITE_PREFIX"))
     255        strcpy (mat_pref, getenv ("IPOPT_WRITE_PREFIX"));
     256      else
     257        strcpy (mat_pref, "mat-ipopt");
     258
     259      sprintf (mat_name, "%s_%03d-%02d.iajaa", mat_pref, IpData().iter_count(), debug_cnt);
     260      mat_file = fopen (mat_name, "w");
     261
     262      fprintf (mat_file, "%d\n", N);
     263      fprintf (mat_file, "%d\n", NNZ);
     264
     265      /* Write ia's */
     266      for (i = 0; i < N+1; i++)
     267        fprintf (mat_file, "%d\n", ia[i]);
     268
     269      /* Write ja's */
     270      for (i = 0; i < NNZ; i++)
     271        fprintf (mat_file, "%d\n", ja[i]);
     272
     273      /* Write values */
     274      for (i = 0; i < NNZ; i++)
     275        fprintf (mat_file, "%32.24e\n", a_[i]);
     276
     277      // /* Write RHS */
     278      // for (i = 0; i < N; i++)
     279      //        fprintf (mat_file, "%32.24e\n", rhs_vals[i]);
     280
     281      fclose (mat_file);
     282    }
     283
     284    debug_last_iter = IpData().iter_count();
     285    debug_cnt ++;
     286
     287
     288
     289    if (PHASE_ == 0)
     290      PHASE_ = 11;
    236291    F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_,
    237                               &PHASE, &N, a_, ia, ja, &PERM,
     292                              &PHASE_, &N, a_, ia, ja, &PERM,
    238293                              &NRHS, IPARM_, &MSGLVL_, &B, &X,
    239294                              &ERROR);
     
    247302      return SYMSOLVER_FATAL_ERROR;
    248303    }
    249 
    250     /* ToDo ask Olaf what this means
     304    PHASE_ = 22;
     305
     306    /* ToDo ask Olaf what this means. */
    251307    if (IPARM_[13] != 0) {
    252308      Jnlst().Printf(J_WARNING, J_LINEAR_ALGEBRA,
    253309                     "Number of perturbed pivots in factorization phase = %d.\n", IPARM_[13]);
    254     }
    255     */
     310      PHASE_ = 12;
     311    }
    256312
    257313    negevals_ = IPARM_[22];
     
    286342    ipfint ERROR;
    287343
     344
     345    // // Dump matrix to file...
     346    // ipfint  NNZ = ia[N]-1;
     347    // ipfint   i;
     348
     349    // if (IpData().iter_count() != debug_last_iter) {
     350    //   debug_cnt = 0;
     351    // }
     352
     353    // if (getenv ("IPOPT_WRITE_MAT"))
     354    // {
     355    //   /* Write header */
     356    //   FILE    *mat_file;
     357    //   char     mat_name[128];
     358    //   char     mat_pref[32];
     359
     360    //   if (getenv ("IPOPT_WRITE_PREFIX"))
     361    //  strcpy (mat_pref, getenv ("IPOPT_WRITE_PREFIX"));
     362    //   else
     363    //  strcpy (mat_pref, "mat-ipopt");
     364
     365    //   sprintf (mat_name, "%s_%03d-%02d.iajaa", mat_pref, IpData().iter_count(), debug_cnt);
     366    //   mat_file = fopen (mat_name, "w");
     367
     368    //   fprintf (mat_file, "%d\n", N);
     369    //   fprintf (mat_file, "%d\n", NNZ);
     370
     371    //   /* Write ia's */
     372    //   for (i = 0; i < N+1; i++)
     373    //  fprintf (mat_file, "%d\n", ia[i]);
     374
     375    //   /* Write ja's */
     376    //   for (i = 0; i < NNZ; i++)
     377    //  fprintf (mat_file, "%d\n", ja[i]);
     378
     379    //   /* Write values */
     380    //   for (i = 0; i < NNZ; i++)
     381    //  fprintf (mat_file, "%32.24e\n", a_[i]);
     382
     383    //   /* Write RHS */
     384    //   for (i = 0; i < N; i++)
     385    //  fprintf (mat_file, "%32.24e\n", rhs_vals[i]);
     386
     387    //   fclose (mat_file);
     388    // }
     389
     390    // debug_last_iter = IpData().iter_count();
     391    // debug_cnt ++;
     392
     393
    288394    F77_FUNC(pardiso,PARDISO)(PT_, &MAXFCT_, &MNUM_, &MTYPE_,
    289395                              &PHASE, &N, a_, ia, ja, &PERM,
    290396                              &NRHS, IPARM_, &MSGLVL_, rhs_vals, X,
    291397                              &ERROR);
     398
     399    if (IPARM_[6] != 0) {
     400      Jnlst().Printf(J_WARNING, J_LINEAR_ALGEBRA,
     401                     "Number of iterative refinement steps = %d.\n", IPARM_[6]);
     402    }
     403
    292404    if (ERROR!=0 ) {
    293405      Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
  • branches/dev/Algorithm/LinearSolvers/IpPardisoSolverInterface.hpp

    r501 r531  
    141141    /** Message level. */
    142142    ipfint MSGLVL_;
     143    /** Flag indicating in what phase Pardiso is */
     144    ipfint PHASE_;
    143145    //@}
    144146
Note: See TracChangeset for help on using the changeset viewer.