Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (3 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpLsqr.cpp

    r1940 r2385  
    77#include "ClpPdco.hpp"
    88
    9 
    10 void ClpLsqr::do_lsqr( CoinDenseVector<double> &b,
    11                        double damp, double atol, double btol, double conlim, int itnlim,
    12                        bool show, Info info, CoinDenseVector<double> &x , int *istop,
    13                        int *itn, Outfo *outfo, bool precon, CoinDenseVector<double> &Pr)
    14 {
    15 
    16      /**
     9void ClpLsqr::do_lsqr(CoinDenseVector< double > &b,
     10  double damp, double atol, double btol, double conlim, int itnlim,
     11  bool show, Info info, CoinDenseVector< double > &x, int *istop,
     12  int *itn, Outfo *outfo, bool precon, CoinDenseVector< double > &Pr)
     13{
     14
     15  /**
    1716     Special version of LSQR for use with pdco.m.
    1817     It continues with a reduced atol if a pdco-specific test isn't
     
    2019     */
    2120
    22 //     Initialize.
    23      static char term_msg[8][80] = {
    24           "The exact solution is x = 0",
    25           "The residual Ax - b is small enough, given ATOL and BTOL",
    26           "The least squares error is small enough, given ATOL",
    27           "The estimated condition number has exceeded CONLIM",
    28           "The residual Ax - b is small enough, given machine precision",
    29           "The least squares error is small enough, given machine precision",
    30           "The estimated condition number has exceeded machine precision",
    31           "The iteration limit has been reached"
    32      };
    33 
    34      //  printf("***************** Entering LSQR *************\n");
    35      assert (model_);
    36 
    37      char str1[100], str2[100], str3[100], str4[100], head1[100], head2[100];
    38 
    39      int n = ncols_;     // set  m,n from lsqr object
    40 
    41      *itn     = 0;
    42      *istop   = 0;
    43      double ctol   = 0;
    44      if (conlim > 0) ctol = 1 / conlim;
    45 
    46      double anorm  = 0;
    47      double acond  = 0;
    48      double ddnorm = 0;
    49      double xnorm  = 0;
    50      double xxnorm = 0;
    51      double z      = 0;
    52      double cs2    = -1;
    53      double sn2    = 0;
    54 
    55      // Set up the first vectors u and v for the bidiagonalization.
    56      // These satisfy  beta*u = b,  alfa*v = A'u.
    57 
    58      CoinDenseVector<double> u(b);
    59      CoinDenseVector<double> v(n, 0.0);
    60      x.clear();
    61      double alfa   = 0;
    62      double beta = u.twoNorm();
    63      if (beta > 0) {
    64           u = (1 / beta) * u;
    65           matVecMult( 2, v, u );
    66           if (precon)
    67                v = v * Pr;
    68           alfa = v.twoNorm();
    69      }
    70      if (alfa > 0) {
    71           v.scale(1 / alfa);
    72      }
    73      CoinDenseVector<double> w(v);
    74 
    75      double arnorm = alfa * beta;
    76      if (arnorm == 0) {
    77           printf("  %s\n\n", term_msg[0]);
    78           return;
    79      }
    80 
    81      double rhobar = alfa;
    82      double phibar = beta;
    83      double bnorm  = beta;
    84      double rnorm  = beta;
    85      sprintf(head1, "   Itn      x(1)      Function");
    86      sprintf(head2, " Compatible   LS      Norm A   Cond A");
    87 
    88      if (show) {
    89           printf(" %s%s\n", head1, head2);
    90           double test1  = 1;
    91           double test2  = alfa / beta;
    92           sprintf(str1, "%6d %12.5e %10.3e",   *itn, x[0], rnorm );
    93           sprintf(str2, "  %8.1e  %8.1e",       test1, test2 );
    94           printf("%s%s\n", str1, str2);
    95      }
    96 
    97      //----------------------------------------------------------------
    98      // Main iteration loop.
    99      //----------------------------------------------------------------
    100      while (*itn < itnlim) {
    101           *itn += 1;
    102           // Perform the next step of the bidiagonalization to obtain the
    103           // next beta, u, alfa, v.  These satisfy the relations
    104           // beta*u  =  a*v   -  alfa*u,
    105           // alfa*v  =  A'*u  -  beta*v.
    106 
    107           u.scale((-alfa));
    108           if (precon) {
    109                CoinDenseVector<double> pv(v * Pr);
    110                matVecMult( 1, u, pv);
    111           } else {
    112                matVecMult( 1, u, v);
     21  //     Initialize.
     22  static char term_msg[8][80] = {
     23    "The exact solution is x = 0",
     24    "The residual Ax - b is small enough, given ATOL and BTOL",
     25    "The least squares error is small enough, given ATOL",
     26    "The estimated condition number has exceeded CONLIM",
     27    "The residual Ax - b is small enough, given machine precision",
     28    "The least squares error is small enough, given machine precision",
     29    "The estimated condition number has exceeded machine precision",
     30    "The iteration limit has been reached"
     31  };
     32
     33  //  printf("***************** Entering LSQR *************\n");
     34  assert(model_);
     35
     36  char str1[100], str2[100], str3[100], str4[100], head1[100], head2[100];
     37
     38  int n = ncols_; // set  m,n from lsqr object
     39
     40  *itn = 0;
     41  *istop = 0;
     42  double ctol = 0;
     43  if (conlim > 0)
     44    ctol = 1 / conlim;
     45
     46  double anorm = 0;
     47  double acond = 0;
     48  double ddnorm = 0;
     49  double xnorm = 0;
     50  double xxnorm = 0;
     51  double z = 0;
     52  double cs2 = -1;
     53  double sn2 = 0;
     54
     55  // Set up the first vectors u and v for the bidiagonalization.
     56  // These satisfy  beta*u = b,  alfa*v = A'u.
     57
     58  CoinDenseVector< double > u(b);
     59  CoinDenseVector< double > v(n, 0.0);
     60  x.clear();
     61  double alfa = 0;
     62  double beta = u.twoNorm();
     63  if (beta > 0) {
     64    u = (1 / beta) * u;
     65    matVecMult(2, v, u);
     66    if (precon)
     67      v = v * Pr;
     68    alfa = v.twoNorm();
     69  }
     70  if (alfa > 0) {
     71    v.scale(1 / alfa);
     72  }
     73  CoinDenseVector< double > w(v);
     74
     75  double arnorm = alfa * beta;
     76  if (arnorm == 0) {
     77    printf("  %s\n\n", term_msg[0]);
     78    return;
     79  }
     80
     81  double rhobar = alfa;
     82  double phibar = beta;
     83  double bnorm = beta;
     84  double rnorm = beta;
     85  sprintf(head1, "   Itn      x(1)      Function");
     86  sprintf(head2, " Compatible   LS      Norm A   Cond A");
     87
     88  if (show) {
     89    printf(" %s%s\n", head1, head2);
     90    double test1 = 1;
     91    double test2 = alfa / beta;
     92    sprintf(str1, "%6d %12.5e %10.3e", *itn, x[0], rnorm);
     93    sprintf(str2, "  %8.1e  %8.1e", test1, test2);
     94    printf("%s%s\n", str1, str2);
     95  }
     96
     97  //----------------------------------------------------------------
     98  // Main iteration loop.
     99  //----------------------------------------------------------------
     100  while (*itn < itnlim) {
     101    *itn += 1;
     102    // Perform the next step of the bidiagonalization to obtain the
     103    // next beta, u, alfa, v.  These satisfy the relations
     104    // beta*u  =  a*v   -  alfa*u,
     105    // alfa*v  =  A'*u  -  beta*v.
     106
     107    u.scale((-alfa));
     108    if (precon) {
     109      CoinDenseVector< double > pv(v * Pr);
     110      matVecMult(1, u, pv);
     111    } else {
     112      matVecMult(1, u, v);
     113    }
     114    beta = u.twoNorm();
     115    if (beta > 0) {
     116      u.scale((1 / beta));
     117      anorm = sqrt(anorm * anorm + alfa * alfa + beta * beta + damp * damp);
     118      v.scale((-beta));
     119      CoinDenseVector< double > vv(n);
     120      vv.clear();
     121      matVecMult(2, vv, u);
     122      if (precon)
     123        vv = vv * Pr;
     124      v = v + vv;
     125      alfa = v.twoNorm();
     126      if (alfa > 0)
     127        v.scale((1 / alfa));
     128    }
     129
     130    // Use a plane rotation to eliminate the damping parameter.
     131    // This alters the diagonal (rhobar) of the lower-bidiagonal matrix.
     132
     133    double rhobar1 = sqrt(rhobar * rhobar + damp * damp);
     134    double cs1 = rhobar / rhobar1;
     135    double sn1 = damp / rhobar1;
     136    double psi = sn1 * phibar;
     137    phibar *= cs1;
     138
     139    // Use a plane rotation to eliminate the subdiagonal element (beta)
     140    // of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.
     141
     142    double rho = sqrt(rhobar1 * rhobar1 + beta * beta);
     143    double cs = rhobar1 / rho;
     144    double sn = beta / rho;
     145    double theta = sn * alfa;
     146    rhobar = -cs * alfa;
     147    double phi = cs * phibar;
     148    phibar = sn * phibar;
     149    double tau = sn * phi;
     150
     151    // Update x and w.
     152
     153    double t1 = phi / rho;
     154    double t2 = -theta / rho;
     155    //    dk           =   ((1/rho)*w);
     156
     157    double w_norm = w.twoNorm();
     158    x = x + t1 * w;
     159    w = v + t2 * w;
     160    ddnorm = ddnorm + (w_norm / rho) * (w_norm / rho);
     161    // if wantvar, var = var  +  dk.*dk; end
     162
     163    // Use a plane rotation on the right to eliminate the
     164    // super-diagonal element (theta) of the upper-bidiagonal matrix.
     165    // Then use the result to estimate  norm(x).
     166
     167    double delta = sn2 * rho;
     168    double gambar = -cs2 * rho;
     169    double rhs = phi - delta * z;
     170    double zbar = rhs / gambar;
     171    xnorm = sqrt(xxnorm + zbar * zbar);
     172    double gamma = sqrt(gambar * gambar + theta * theta);
     173    cs2 = gambar / gamma;
     174    sn2 = theta / gamma;
     175    z = rhs / gamma;
     176    xxnorm = xxnorm + z * z;
     177
     178    // Test for convergence.
     179    // First, estimate the condition of the matrix  Abar,
     180    // and the norms of  rbar  and  Abar'rbar.
     181
     182    acond = anorm * sqrt(ddnorm);
     183    double res1 = phibar * phibar;
     184    double res2 = res1 + psi * psi;
     185    rnorm = sqrt(res1 + res2);
     186    arnorm = alfa * fabs(tau);
     187
     188    // Now use these norms to estimate certain other quantities,
     189    // some of which will be small near a solution.
     190
     191    double test1 = rnorm / bnorm;
     192    double test2 = arnorm / (anorm * rnorm);
     193    double test3 = 1 / acond;
     194    t1 = test1 / (1 + anorm * xnorm / bnorm);
     195    double rtol = btol + atol * anorm * xnorm / bnorm;
     196
     197    // The following tests guard against extremely small values of
     198    // atol, btol  or  ctol.  (The user may have set any or all of
     199    // the parameters  atol, btol, conlim  to 0.)
     200    // The effect is equivalent to the normal tests using
     201    // atol = eps,  btol = eps,  conlim = 1/eps.
     202
     203    if (*itn >= itnlim)
     204      *istop = 7;
     205    if (1 + test3 <= 1)
     206      *istop = 6;
     207    if (1 + test2 <= 1)
     208      *istop = 5;
     209    if (1 + t1 <= 1)
     210      *istop = 4;
     211
     212    // Allow for tolerances set by the user.
     213
     214    if (test3 <= ctol)
     215      *istop = 3;
     216    if (test2 <= atol)
     217      *istop = 2;
     218    if (test1 <= rtol)
     219      *istop = 1;
     220
     221    //-------------------------------------------------------------------
     222    // SPECIAL TEST THAT DEPENDS ON pdco.m.
     223    // Aname in pdco   is  iw in lsqr.
     224    // dy              is  x
     225    // Other stuff     is in info.
     226
     227    // We allow for diagonal preconditioning in pdDDD3.
     228    //-------------------------------------------------------------------
     229    if (*istop > 0) {
     230      double r3new = arnorm;
     231      double r3ratio = r3new / info.r3norm;
     232      double atolold = atol;
     233      double atolnew = atol;
     234
     235      if (atol > info.atolmin) {
     236        if (r3ratio <= 0.1) { // dy seems good
     237          // Relax
     238        } else if (r3ratio <= 0.5) { // Accept dy but make next one more accurate.
     239          atolnew = atolnew * 0.1;
     240        } else { // Recompute dy more accurately
     241          if (show) {
     242            printf("\n                                ");
     243            printf("                                \n");
     244            printf(" %5.1f%7d%7.3f", log10(atolold), *itn, r3ratio);
    113245          }
    114           beta = u.twoNorm();
    115           if (beta > 0) {
    116                u.scale((1 / beta));
    117                anorm = sqrt(anorm * anorm + alfa * alfa + beta * beta +  damp * damp);
    118                v.scale((-beta));
    119                CoinDenseVector<double> vv(n);
    120                vv.clear();
    121                matVecMult( 2, vv, u );
    122                if (precon)
    123                     vv = vv * Pr;
    124                v = v + vv;
    125                alfa  = v.twoNorm();
    126                if (alfa > 0)
    127                     v.scale((1 / alfa));
    128           }
    129 
    130           // Use a plane rotation to eliminate the damping parameter.
    131           // This alters the diagonal (rhobar) of the lower-bidiagonal matrix.
    132 
    133           double rhobar1 = sqrt(rhobar * rhobar + damp * damp);
    134           double cs1     = rhobar / rhobar1;
    135           double sn1     = damp   / rhobar1;
    136           double psi     = sn1 * phibar;
    137           phibar       *= cs1;
    138 
    139           // Use a plane rotation to eliminate the subdiagonal element (beta)
    140           // of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.
    141 
    142           double rho     = sqrt(rhobar1 * rhobar1 +  beta * beta);
    143           double cs      =   rhobar1 / rho;
    144           double sn      =   beta   / rho;
    145           double theta   =   sn * alfa;
    146           rhobar  = - cs * alfa;
    147           double phi     =   cs * phibar;
    148           phibar  =   sn * phibar;
    149           double tau     =   sn * phi;
    150 
    151           // Update x and w.
    152 
    153           double t1      =   phi  / rho;
    154           double t2      = - theta / rho;
    155           //    dk           =   ((1/rho)*w);
    156 
    157           double w_norm = w.twoNorm();
    158           x       = x      +  t1 * w;
    159           w       = v      +  t2 * w;
    160           ddnorm  = ddnorm +  (w_norm / rho) * (w_norm / rho);
    161           // if wantvar, var = var  +  dk.*dk; end
    162 
    163           // Use a plane rotation on the right to eliminate the
    164           // super-diagonal element (theta) of the upper-bidiagonal matrix.
    165           // Then use the result to estimate  norm(x).
    166 
    167           double delta   =   sn2 * rho;
    168           double gambar  = - cs2 * rho;
    169           double rhs     =   phi  -  delta * z;
    170           double zbar    =   rhs / gambar;
    171           xnorm   =   sqrt(xxnorm + zbar * zbar);
    172           double gamma   =   sqrt(gambar * gambar + theta * theta);
    173           cs2     =   gambar / gamma;
    174           sn2     =   theta  / gamma;
    175           z       =   rhs    / gamma;
    176           xxnorm  =   xxnorm  +  z * z;
    177 
    178           // Test for convergence.
    179           // First, estimate the condition of the matrix  Abar,
    180           // and the norms of  rbar  and  Abar'rbar.
    181 
    182           acond   =   anorm * sqrt( ddnorm );
    183           double res1    =   phibar * phibar;
    184           double res2    =   res1  +  psi * psi;
    185           rnorm   =   sqrt( res1 + res2 );
    186           arnorm  =   alfa * fabs( tau );
    187 
    188           // Now use these norms to estimate certain other quantities,
    189           // some of which will be small near a solution.
    190 
    191           double test1   =   rnorm / bnorm;
    192           double test2   =   arnorm / ( anorm * rnorm );
    193           double test3   =       1 / acond;
    194           t1      =   test1 / (1    +  anorm * xnorm / bnorm);
    195           double rtol    =   btol  +  atol *  anorm * xnorm / bnorm;
    196 
    197           // The following tests guard against extremely small values of
    198           // atol, btol  or  ctol.  (The user may have set any or all of
    199           // the parameters  atol, btol, conlim  to 0.)
    200           // The effect is equivalent to the normal tests using
    201           // atol = eps,  btol = eps,  conlim = 1/eps.
    202 
    203           if (*itn >= itnlim)
    204                *istop = 7;
    205           if (1 + test3  <= 1)
    206                *istop = 6;
    207           if (1 + test2  <= 1)
    208                *istop = 5;
    209           if (1 + t1     <= 1)
    210                *istop = 4;
    211 
    212           // Allow for tolerances set by the user.
    213 
    214           if  (test3 <= ctol)
    215                *istop = 3;
    216           if  (test2 <= atol)
    217                *istop = 2;
    218           if  (test1 <= rtol)
    219                *istop = 1;
    220 
    221           //-------------------------------------------------------------------
    222           // SPECIAL TEST THAT DEPENDS ON pdco.m.
    223           // Aname in pdco   is  iw in lsqr.
    224           // dy              is  x
    225           // Other stuff     is in info.
    226 
    227           // We allow for diagonal preconditioning in pdDDD3.
    228           //-------------------------------------------------------------------
    229           if (*istop > 0) {
    230                double r3new     = arnorm;
    231                double r3ratio   = r3new / info.r3norm;
    232                double atolold   = atol;
    233                double atolnew   = atol;
    234 
    235                if (atol > info.atolmin) {
    236                     if     (r3ratio <= 0.1) {  // dy seems good
    237                          // Relax
    238                     } else if (r3ratio <= 0.5) { // Accept dy but make next one more accurate.
    239                          atolnew = atolnew * 0.1;
    240                     } else {                    // Recompute dy more accurately
    241                          if (show) {
    242                               printf("\n                                ");
    243                               printf("                                \n");
    244                               printf(" %5.1f%7d%7.3f", log10(atolold), *itn, r3ratio);
    245                          }
    246                          atol    = atol * 0.1;
    247                          atolnew = atol;
    248                          *istop   = 0;
    249                     }
    250 
    251                     outfo->atolold = atolold;
    252                     outfo->atolnew = atolnew;
    253                     outfo->r3ratio = r3ratio;
    254                }
    255 
    256                //-------------------------------------------------------------------
    257                // See if it is time to print something.
    258                //-------------------------------------------------------------------
    259                int prnt = 0;
    260                if (n     <= 40       ) prnt = 1;
    261                if (*itn   <= 10       ) prnt = 1;
    262                if (*itn   >= itnlim - 10) prnt = 1;
    263                if (*itn % 10 == 0  ) prnt = 1;
    264                if (test3 <=  2 * ctol  ) prnt = 1;
    265                if (test2 <= 10 * atol  ) prnt = 1;
    266                if (test1 <= 10 * rtol  ) prnt = 1;
    267                if (*istop !=  0       ) prnt = 1;
    268 
    269                if (prnt == 1) {
    270                     if (show) {
    271                          sprintf(str1, "   %6d %12.5e %10.3e",   *itn, x[0], rnorm );
    272                          sprintf(str2, "  %8.1e %8.1e",       test1, test2 );
    273                          sprintf(str3, " %8.1e %8.1e",        anorm, acond );
    274                          printf("%s%s%s\n", str1, str2, str3);
    275                     }
    276                }
    277                if (*istop > 0)
    278                     break;
    279           }
    280      }
    281      // End of iteration loop.
    282      // Print the stopping condition.
    283 
    284      if (show) {
    285           printf("\n LSQR finished\n");
    286           //    disp(msg(istop+1,:))
    287           //    disp(' ')
    288           printf("%s\n", term_msg[*istop]);
    289           sprintf(str1, "istop  =%8d     itn    =%8d",      *istop, *itn    );
    290           sprintf(str2, "anorm  =%8.1e   acond  =%8.1e",  anorm, acond  );
    291           sprintf(str3, "rnorm  =%8.1e   arnorm =%8.1e",  rnorm, arnorm );
    292           sprintf(str4, "bnorm  =%8.1e   xnorm  =%8.1e",  bnorm, xnorm  );
    293           printf("%s %s\n", str1, str2);
    294           printf("%s %s\n", str3, str4);
    295      }
    296 }
    297 
    298 
    299 
    300 
    301 void ClpLsqr::matVecMult( int mode, CoinDenseVector<double> *x, CoinDenseVector<double> *y)
    302 {
    303      int n = model_->numberColumns();
    304      int m = model_->numberRows();
    305      CoinDenseVector<double> *temp = new CoinDenseVector<double>(n, 0.0);
    306      double *t_elts = temp->getElements();
    307      double *x_elts = x->getElements();
    308      double *y_elts = y->getElements();
    309      ClpPdco * pdcoModel = (ClpPdco *) model_;
    310      if (mode == 1) {
    311           pdcoModel->matVecMult( 2, temp, y);
    312           for (int k = 0; k < n; k++)
    313                x_elts[k] += (diag1_[k] * t_elts[k]);
    314           for (int k = 0; k < m; k++)
    315                x_elts[n+k] += (diag2_ * y_elts[k]);
    316      } else {
    317           for (int k = 0; k < n; k++)
    318                t_elts[k] = diag1_[k] * y_elts[k];
    319           pdcoModel->matVecMult( 1, x, temp);
    320           for (int k = 0; k < m; k++)
    321                x_elts[k] += diag2_ * y_elts[n+k];
    322      }
    323      delete temp;
    324      return;
    325 }
    326 
    327 void ClpLsqr::matVecMult( int mode, CoinDenseVector<double> &x, CoinDenseVector<double> &y)
    328 {
    329      matVecMult( mode, &x, &y);
    330      return;
     246          atol = atol * 0.1;
     247          atolnew = atol;
     248          *istop = 0;
     249        }
     250
     251        outfo->atolold = atolold;
     252        outfo->atolnew = atolnew;
     253        outfo->r3ratio = r3ratio;
     254      }
     255
     256      //-------------------------------------------------------------------
     257      // See if it is time to print something.
     258      //-------------------------------------------------------------------
     259      int prnt = 0;
     260      if (n <= 40)
     261        prnt = 1;
     262      if (*itn <= 10)
     263        prnt = 1;
     264      if (*itn >= itnlim - 10)
     265        prnt = 1;
     266      if (*itn % 10 == 0)
     267        prnt = 1;
     268      if (test3 <= 2 * ctol)
     269        prnt = 1;
     270      if (test2 <= 10 * atol)
     271        prnt = 1;
     272      if (test1 <= 10 * rtol)
     273        prnt = 1;
     274      if (*istop != 0)
     275        prnt = 1;
     276
     277      if (prnt == 1) {
     278        if (show) {
     279          sprintf(str1, "   %6d %12.5e %10.3e", *itn, x[0], rnorm);
     280          sprintf(str2, "  %8.1e %8.1e", test1, test2);
     281          sprintf(str3, " %8.1e %8.1e", anorm, acond);
     282          printf("%s%s%s\n", str1, str2, str3);
     283        }
     284      }
     285      if (*istop > 0)
     286        break;
     287    }
     288  }
     289  // End of iteration loop.
     290  // Print the stopping condition.
     291
     292  if (show) {
     293    printf("\n LSQR finished\n");
     294    //    disp(msg(istop+1,:))
     295    //    disp(' ')
     296    printf("%s\n", term_msg[*istop]);
     297    sprintf(str1, "istop  =%8d     itn    =%8d", *istop, *itn);
     298    sprintf(str2, "anorm  =%8.1e   acond  =%8.1e", anorm, acond);
     299    sprintf(str3, "rnorm  =%8.1e   arnorm =%8.1e", rnorm, arnorm);
     300    sprintf(str4, "bnorm  =%8.1e   xnorm  =%8.1e", bnorm, xnorm);
     301    printf("%s %s\n", str1, str2);
     302    printf("%s %s\n", str3, str4);
     303  }
     304}
     305
     306void ClpLsqr::matVecMult(int mode, CoinDenseVector< double > *x, CoinDenseVector< double > *y)
     307{
     308  int n = model_->numberColumns();
     309  int m = model_->numberRows();
     310  CoinDenseVector< double > *temp = new CoinDenseVector< double >(n, 0.0);
     311  double *t_elts = temp->getElements();
     312  double *x_elts = x->getElements();
     313  double *y_elts = y->getElements();
     314  ClpPdco *pdcoModel = (ClpPdco *)model_;
     315  if (mode == 1) {
     316    pdcoModel->matVecMult(2, temp, y);
     317    for (int k = 0; k < n; k++)
     318      x_elts[k] += (diag1_[k] * t_elts[k]);
     319    for (int k = 0; k < m; k++)
     320      x_elts[n + k] += (diag2_ * y_elts[k]);
     321  } else {
     322    for (int k = 0; k < n; k++)
     323      t_elts[k] = diag1_[k] * y_elts[k];
     324    pdcoModel->matVecMult(1, x, temp);
     325    for (int k = 0; k < m; k++)
     326      x_elts[k] += diag2_ * y_elts[n + k];
     327  }
     328  delete temp;
     329  return;
     330}
     331
     332void ClpLsqr::matVecMult(int mode, CoinDenseVector< double > &x, CoinDenseVector< double > &y)
     333{
     334  matVecMult(mode, &x, &y);
     335  return;
    331336}
    332337/* Default constructor */
    333 ClpLsqr::ClpLsqr() :
    334      nrows_(0),
    335      ncols_(0),
    336      model_(NULL),
    337      diag1_(NULL),
    338      diag2_(0.0)
    339 {}
     338ClpLsqr::ClpLsqr()
     339  : nrows_(0)
     340  , ncols_(0)
     341  , model_(NULL)
     342  , diag1_(NULL)
     343  , diag2_(0.0)
     344{
     345}
    340346
    341347/* Constructor for use with Pdco model (note modified for pdco!!!!) */
    342 ClpLsqr::ClpLsqr(ClpInterior *model) :
    343      diag1_(NULL),
    344     diag2_(0.0)
    345 {
    346      model_ = model;
    347      nrows_ = model->numberRows() + model->numberColumns();
    348      ncols_ = model->numberRows();
     348ClpLsqr::ClpLsqr(ClpInterior *model)
     349  : diag1_(NULL)
     350  , diag2_(0.0)
     351{
     352  model_ = model;
     353  nrows_ = model->numberRows() + model->numberColumns();
     354  ncols_ = model->numberRows();
    349355}
    350356/** Destructor */
    351357ClpLsqr::~ClpLsqr()
    352358{
    353      // delete [] diag1_; no as we just borrowed it
    354 }
    355 bool
    356 ClpLsqr::setParam(char *parmName, int parmValue)
    357 {
    358      std::cout << "Set lsqr integer parameter " << parmName << "to " << parmValue
    359                << std::endl;
    360      if ( strcmp(parmName, "nrows") == 0) {
    361           nrows_ = parmValue;
    362           return 1;
    363      } else if ( strcmp(parmName, "ncols") == 0) {
    364           ncols_ = parmValue;
    365           return 1;
    366      }
    367      std::cout << "Attempt to set unknown integer parameter name " << parmName << std::endl;
    368      return 0;
    369 }
    370 ClpLsqr::ClpLsqr(const ClpLsqr &rhs) :
    371      nrows_(rhs.nrows_),
    372      ncols_(rhs.ncols_),
    373      model_(rhs.model_),
    374      diag2_(rhs.diag2_)
    375 {
    376      diag1_ = ClpCopyOfArray(rhs.diag1_, nrows_);
     359  // delete [] diag1_; no as we just borrowed it
     360}
     361bool ClpLsqr::setParam(char *parmName, int parmValue)
     362{
     363  std::cout << "Set lsqr integer parameter " << parmName << "to " << parmValue
     364            << std::endl;
     365  if (strcmp(parmName, "nrows") == 0) {
     366    nrows_ = parmValue;
     367    return 1;
     368  } else if (strcmp(parmName, "ncols") == 0) {
     369    ncols_ = parmValue;
     370    return 1;
     371  }
     372  std::cout << "Attempt to set unknown integer parameter name " << parmName << std::endl;
     373  return 0;
     374}
     375ClpLsqr::ClpLsqr(const ClpLsqr &rhs)
     376  : nrows_(rhs.nrows_)
     377  , ncols_(rhs.ncols_)
     378  , model_(rhs.model_)
     379  , diag2_(rhs.diag2_)
     380{
     381  diag1_ = ClpCopyOfArray(rhs.diag1_, nrows_);
    377382}
    378383// Assignment operator. This copies the data
    379384ClpLsqr &
    380 ClpLsqr::operator=(const ClpLsqr & rhs)
    381 {
    382      if (this != &rhs) {
    383           delete [] diag1_;
    384           diag1_ = ClpCopyOfArray(rhs.diag1_, nrows_);
    385           nrows_ = rhs.nrows_;
    386           ncols_ = rhs.ncols_;
    387           model_ = rhs.model_;
    388           diag2_ = rhs.diag2_;
    389      }
    390      return *this;
    391 }
     385ClpLsqr::operator=(const ClpLsqr &rhs)
     386{
     387  if (this != &rhs) {
     388    delete[] diag1_;
     389    diag1_ = ClpCopyOfArray(rhs.diag1_, nrows_);
     390    nrows_ = rhs.nrows_;
     391    ncols_ = rhs.ncols_;
     392    model_ = rhs.model_;
     393    diag2_ = rhs.diag2_;
     394  }
     395  return *this;
     396}
     397
     398/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     399*/
Note: See TracChangeset for help on using the changeset viewer.