Changeset 425


Ignore:
Timestamp:
Sep 5, 2004 3:52:01 AM (15 years ago)
Author:
forrest
Message:

Improve cholesky

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpCholeskyBase.cpp

    r424 r425  
    15291529  //perturbation
    15301530  longDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
    1531   perturbation=perturbation*perturbation;
     1531  //perturbation=perturbation*perturbation*100000000.0;
    15321532  if (perturbation>1.0) {
    15331533    //if (model_->model()->logLevel()&4)
    1534       std::cout <<"large perturbation "<<perturbation<<std::endl;
     1534    //std::cout <<"large perturbation "<<perturbation<<std::endl;
    15351535    perturbation=sqrt(perturbation);;
    15361536    perturbation=1.0;
     
    16361636          rowsDropped[iRow]=2;
    16371637          numberDroppedBefore++;
     1638          //printf("dropped - small diagonal %g\n",diagonal);
    16381639        }
    16391640      }
     
    23852386    if (start<end) {
    23862387      CoinBigIndex offset = indexStart_[iRow]-choleskyStart_[iRow];
    2387       longDouble dValue=d[iRow];
    2388       for (CoinBigIndex k=start;k<end;k++) {
    2389         int kRow = choleskyRow_[k+offset];
    2390         assert(kRow>=firstDense_);
    2391         longDouble a_ik=sparseFactor_[k];
    2392         longDouble value1=dValue*a_ik;
    2393         diagonal_[kRow] -= value1*a_ik;
    2394         CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
    2395         for (CoinBigIndex j=k+1;j<end;j++) {
    2396           int jRow=choleskyRow_[j+offset];
    2397           longDouble a_jk = sparseFactor_[j];
    2398           sparseFactor_[base+jRow] -= a_jk*value1;
    2399         }
     2388      if (clique_[iRow]<2) {
     2389        longDouble dValue=d[iRow];
     2390        for (CoinBigIndex k=start;k<end;k++) {
     2391          int kRow = choleskyRow_[k+offset];
     2392          assert(kRow>=firstDense_);
     2393          longDouble a_ik=sparseFactor_[k];
     2394          longDouble value1=dValue*a_ik;
     2395          diagonal_[kRow] -= value1*a_ik;
     2396          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
     2397          for (CoinBigIndex j=k+1;j<end;j++) {
     2398            int jRow=choleskyRow_[j+offset];
     2399            longDouble a_jk = sparseFactor_[j];
     2400            sparseFactor_[base+jRow] -= a_jk*value1;
     2401          }
     2402        }
     2403      } else if (clique_[iRow]<3) {
     2404        // do as pair
     2405        longDouble dValue0=d[iRow];
     2406        longDouble dValue1=d[iRow+1];
     2407        int offset1 = first[iRow+1]-start;
     2408        // skip row
     2409        iRow++;
     2410        for (CoinBigIndex k=start;k<end;k++) {
     2411          int kRow = choleskyRow_[k+offset];
     2412          assert(kRow>=firstDense_);
     2413          longDouble a_ik0=sparseFactor_[k];
     2414          longDouble value0=dValue0*a_ik0;
     2415          longDouble a_ik1=sparseFactor_[k+offset1];
     2416          longDouble value1=dValue1*a_ik1;
     2417          diagonal_[kRow] -= value0*a_ik0 + value1*a_ik1;
     2418          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
     2419          for (CoinBigIndex j=k+1;j<end;j++) {
     2420            int jRow=choleskyRow_[j+offset];
     2421            longDouble a_jk0 = sparseFactor_[j];
     2422            longDouble a_jk1 = sparseFactor_[j+offset1];
     2423            sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1;
     2424          }
     2425        }
     2426#define MANY_REGISTERS
     2427#ifdef MANY_REGISTERS
     2428      } else if (clique_[iRow]==3) {
     2429#else
     2430      } else {
     2431#endif
     2432        // do as clique
     2433        // maybe later get fancy on big cliques and do transpose copy
     2434        // seems only worth going to 3 on Intel
     2435        longDouble dValue0=d[iRow];
     2436        longDouble dValue1=d[iRow+1];
     2437        longDouble dValue2=d[iRow+2];
     2438        // get offsets and skip rows
     2439        int offset1 = first[++iRow]-start;
     2440        int offset2 = first[++iRow]-start;
     2441        for (CoinBigIndex k=start;k<end;k++) {
     2442          int kRow = choleskyRow_[k+offset];
     2443          assert(kRow>=firstDense_);
     2444          double diagonalValue=diagonal_[kRow];
     2445          longDouble a_ik0=sparseFactor_[k];
     2446          longDouble value0=dValue0*a_ik0;
     2447          longDouble a_ik1=sparseFactor_[k+offset1];
     2448          longDouble value1=dValue1*a_ik1;
     2449          longDouble a_ik2=sparseFactor_[k+offset2];
     2450          longDouble value2=dValue2*a_ik2;
     2451          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
     2452          diagonal_[kRow] = diagonalValue - value0*a_ik0 - value1*a_ik1 - value2*a_ik2;
     2453          for (CoinBigIndex j=k+1;j<end;j++) {
     2454            int jRow=choleskyRow_[j+offset];
     2455            longDouble a_jk0 = sparseFactor_[j];
     2456            longDouble a_jk1 = sparseFactor_[j+offset1];
     2457            longDouble a_jk2 = sparseFactor_[j+offset2];
     2458            sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1+a_jk2*value2;
     2459          }
     2460        }
     2461#ifdef MANY_REGISTERS
     2462      } else {
     2463        // do as clique
     2464        // maybe later get fancy on big cliques and do transpose copy
     2465        // maybe only worth going to 3 on Intel (but may have hidden registers)
     2466        longDouble dValue0=d[iRow];
     2467        longDouble dValue1=d[iRow+1];
     2468        longDouble dValue2=d[iRow+2];
     2469        longDouble dValue3=d[iRow+3];
     2470        // get offsets and skip rows
     2471        int offset1 = first[++iRow]-start;
     2472        int offset2 = first[++iRow]-start;
     2473        int offset3 = first[++iRow]-start;
     2474        for (CoinBigIndex k=start;k<end;k++) {
     2475          int kRow = choleskyRow_[k+offset];
     2476          assert(kRow>=firstDense_);
     2477          double diagonalValue=diagonal_[kRow];
     2478          longDouble a_ik0=sparseFactor_[k];
     2479          longDouble value0=dValue0*a_ik0;
     2480          longDouble a_ik1=sparseFactor_[k+offset1];
     2481          longDouble value1=dValue1*a_ik1;
     2482          longDouble a_ik2=sparseFactor_[k+offset2];
     2483          longDouble value2=dValue2*a_ik2;
     2484          longDouble a_ik3=sparseFactor_[k+offset3];
     2485          longDouble value3=dValue3*a_ik3;
     2486          CoinBigIndex base = choleskyStart_[kRow]-kRow-1;
     2487          diagonal_[kRow] = diagonalValue - (value0*a_ik0 + value1*a_ik1 + value2*a_ik2+value3*a_ik3);
     2488          for (CoinBigIndex j=k+1;j<end;j++) {
     2489            int jRow=choleskyRow_[j+offset];
     2490            longDouble a_jk0 = sparseFactor_[j];
     2491            longDouble a_jk1 = sparseFactor_[j+offset1];
     2492            longDouble a_jk2 = sparseFactor_[j+offset2];
     2493            longDouble a_jk3 = sparseFactor_[j+offset3];
     2494            sparseFactor_[base+jRow] -= a_jk0*value0+a_jk1*value1+a_jk2*value2+a_jk3*value3;
     2495          }
     2496        }
     2497#endif
    24002498      }
    24012499    }
  • trunk/ClpCholeskyDense.cpp

    r423 r425  
    10211021  aa = aOther-4*BLOCK;
    10221022  if (nUnder==BLOCK) {
     1023    //#define INTEL
     1024#ifdef INTEL
     1025    aa+=2*BLOCK;
     1026    for (j = 0; j < BLOCK; j +=2) {
     1027      aa+=2*BLOCK;
     1028      for (i=0;i< BLOCK;i+=2) {
     1029        longDouble t00=aa[i+0*BLOCK];
     1030        longDouble t10=aa[i+1*BLOCK];
     1031        longDouble t01=aa[i+1+0*BLOCK];
     1032        longDouble t11=aa[i+1+1*BLOCK];
     1033        for (k=0;k<BLOCK;k++) {
     1034          longDouble multiplier = work[k];
     1035          longDouble a00=aUnder[i+k*BLOCK]*multiplier;
     1036          longDouble a01=aUnder[i+1+k*BLOCK]*multiplier;
     1037          t00 -= a00 * above[j + 0 + k * BLOCK];
     1038          t10 -= a00 * above[j + 1 + k * BLOCK];
     1039          t01 -= a01 * above[j + 0 + k * BLOCK];
     1040          t11 -= a01 * above[j + 1 + k * BLOCK];
     1041        }
     1042        aa[i+0*BLOCK]=t00;
     1043        aa[i+1*BLOCK]=t10;
     1044        aa[i+1+0*BLOCK]=t01;
     1045        aa[i+1+1*BLOCK]=t11;
     1046      }
     1047    }
     1048#else
    10231049    for (j = 0; j < BLOCK; j +=4) {
    10241050      aa+=4*BLOCK;
     
    10811107      }
    10821108    }
     1109#endif
    10831110  } else {
    10841111    int odd=nUnder&1;
  • trunk/ClpPredictorCorrector.cpp

    r424 r425  
    2424#include <iostream>
    2525
    26 static double eScale=1.0e57;
     26static double eScale=1.0e27;
    2727static double eBaseCaution=1.0e-12;
    2828static double eBase=1.0e-12;
  • trunk/ClpSolve.cpp

    r423 r425  
    10371037    //barrier.setGamma(1.0e-8);
    10381038    //barrier.setDelta(1.0e-8);
     1039    barrier.setDiagonalPerturbation(1.0e-14);
    10391040    if (aggressiveGamma) {
    10401041      barrier.setGamma(1.0e-3);
Note: See TracChangeset for help on using the changeset viewer.