Changeset 1368 for trunk/Clp


Ignore:
Timestamp:
May 30, 2009 4:52:38 AM (11 years ago)
Author:
forrest
Message:

changes for cholesky including code from Anshul Gupta

Location:
trunk/Clp/src
Files:
9 edited

Legend:

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

    r1367 r1368  
    11// Copyright (C) 2002, International Business Machines
    22// Corporation and others.  All Rights Reserved.
     3/*----------------------------------------------------------------------------*/
     4/*      Ordering code - courtesy of Anshul Gupta                              */
     5/*      (C) Copyright IBM Corporation 1997, 2009.  All Rights Reserved.       */
     6/*----------------------------------------------------------------------------*/
     7
     8/*  A compact no-frills Approximate Minimum Local Fill ordering code.
     9
     10    References:
     11
     12[1] Ordering Sparse Matrices Using Approximate Minimum Local Fill.
     13    Edward Rothberg, SGI Manuscript, April 1996.
     14[2] An Approximate Minimum Degree Ordering Algorithm.
     15    T. Davis, P. Amestoy, and I. Duff, TR-94-039, CIS Department,
     16    University of Florida, December 1994.
     17*/
     18/*----------------------------------------------------------------------------*/
     19
    320
    421#include "CoinPragma.hpp"
    522
    6 #include <iostream>
     23#include <iostream> 
    724
    825#include "ClpCholeskyBase.hpp"
     
    621638{
    622639  model_=model;
     640#define BASE_ORDER 2
     641#if BASE_ORDER>0
     642  if (!doKKT_&&model_->numberRows()>6) {
     643    if (preOrder(false,true,false))
     644      return -1;
     645    //rowsDropped_ = new char [numberRows_];
     646    numberRowsDropped_=0;
     647    memset(rowsDropped_,0,numberRows_);
     648    //rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
     649    // approximate minimum degree
     650    return orderAMD();
     651  }
     652#endif
    623653  int numberRowsModel = model_->numberRows();
    624654  int numberColumns = model_->numberColumns();
     
    636666  rowsDropped_ = new char [numberRows_];
    637667  numberRowsDropped_=0;
     668  memset(rowsDropped_,0,numberRows_);
    638669  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
    639670  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
     
    749780  delete [] count;
    750781  permuteInverse_ = new int [numberRows_];
    751   memset(rowsDropped_,0,numberRows_);
    752782  for (iRow=0;iRow<numberRows_;iRow++) {
    753783    //permute_[iRow]=iRow; // force no permute
     
    757787  }
    758788  return 0;
     789}
     790#if BASE_ORDER==1
     791/* Orders rows and saves pointer to matrix.and model */
     792int
     793ClpCholeskyBase::orderAMD()
     794{
     795  permuteInverse_ = new int [numberRows_];
     796  permute_ = new int[numberRows_];
     797  // Do ordering
     798  int returnCode=0;
     799  // get more space and full matrix
     800  int space = 6*sizeFactor_+100000;
     801  int * temp = new int [space];
     802  int * which = new int[2*numberRows_];
     803  CoinBigIndex * tempStart = new CoinBigIndex [numberRows_+1];
     804  memset(which,0,numberRows_*sizeof(int));
     805  for (int iRow=0;iRow<numberRows_;iRow++) {
     806    which[iRow]+= choleskyStart_[iRow+1]-choleskyStart_[iRow]-1;
     807    for (CoinBigIndex j=choleskyStart_[iRow]+1;j<choleskyStart_[iRow+1];j++) {
     808      int jRow=choleskyRow_[j];
     809      which[jRow]++;
     810    }
     811  }
     812  CoinBigIndex sizeFactor =0;
     813  for (int iRow=0;iRow<numberRows_;iRow++) {
     814    int length = which[iRow];
     815    permute_[iRow]=length;
     816    tempStart[iRow]=sizeFactor;
     817    which[iRow]=sizeFactor;
     818    sizeFactor += length;
     819  }
     820  for (int iRow=0;iRow<numberRows_;iRow++) {
     821    assert (choleskyRow_[choleskyStart_[iRow]]==iRow);
     822    for (CoinBigIndex j=choleskyStart_[iRow]+1;j<choleskyStart_[iRow+1];j++) {
     823      int jRow=choleskyRow_[j];
     824      int put = which[iRow];
     825      temp[put]=jRow;
     826      which[iRow]++;
     827      put = which[jRow];
     828      temp[put]=iRow;
     829      which[jRow]++;
     830    }
     831  }
     832  for (int iRow=1;iRow<numberRows_;iRow++)
     833    assert (which[iRow-1]==tempStart[iRow]);
     834  CoinBigIndex lastSpace=sizeFactor;
     835  delete [] choleskyRow_;
     836  choleskyRow_=temp;
     837  delete [] choleskyStart_;
     838  choleskyStart_=tempStart;
     839  // linked lists of sizes and lengths
     840  int * first = new int [numberRows_];
     841  int * next = new int [numberRows_];
     842  int * previous = new int [numberRows_];
     843  char * mark = new char[numberRows_];
     844  memset(mark,0,numberRows_);
     845  CoinBigIndex * sort = new CoinBigIndex [numberRows_];
     846  int * order = new int [numberRows_];
     847  for (int iRow=0;iRow<numberRows_;iRow++) {
     848    first[iRow]=-1;
     849    next[iRow]=-1;
     850    previous[iRow]=-1;
     851    permuteInverse_[iRow]=-1;
     852  }
     853  int large = 1000+2*numberRows_;
     854  for (int iRow=0;iRow<numberRows_;iRow++) {
     855    int n = permute_[iRow];
     856    if (first[n]<0) {
     857      first[n]=iRow;
     858      previous[iRow]=n+large;
     859      next[iRow]=n+2*large;
     860    } else {
     861      int k=first[n];
     862      assert (k<numberRows_);
     863      first[n]=iRow;
     864      previous[iRow]=n+large;
     865      assert (previous[k]==n+large);
     866      next[iRow]=k;
     867      previous[k]=iRow;
     868    }
     869  }
     870  int smallest=0;
     871  int done=0;
     872  int numberCompressions=0;
     873  int numberExpansions=0;
     874  while (smallest<numberRows_) {
     875    if (first[smallest]<0||first[smallest]>numberRows_) {
     876      smallest++;
     877      continue;
     878    }
     879    int iRow=first[smallest];
     880    if (false) {
     881      int kRow=-1;
     882      int ss=999999;
     883      for (int jRow=numberRows_-1;jRow>=0;jRow--) {
     884        if (permuteInverse_[jRow]<0) {
     885          int length = permute_[jRow];
     886          assert (length>0);
     887          for (CoinBigIndex j=choleskyStart_[jRow];
     888               j<choleskyStart_[jRow]+length;j++) {
     889            int jjRow=choleskyRow_[j];
     890            assert (permuteInverse_[jjRow]<0);
     891          }
     892          if (length<ss) {
     893            ss=length;
     894            kRow=jRow;
     895          }
     896        }
     897      }
     898      assert (smallest==ss);
     899      printf("list chose %d - full chose %d - length %d\n",
     900             iRow,kRow,ss);
     901    }
     902    int kNext = next[iRow];
     903    first[smallest]=kNext;
     904    if (kNext<numberRows_)
     905      previous[kNext]=previous[iRow];
     906    previous[iRow]=-1;
     907    next[iRow]=-1;
     908    permuteInverse_[iRow]=done;
     909    done++;
     910    // Now add edges
     911    CoinBigIndex start=choleskyStart_[iRow];
     912    CoinBigIndex end=choleskyStart_[iRow]+permute_[iRow];
     913    int nSave=0;
     914    for (CoinBigIndex k=start;k<end;k++) {
     915      int kRow=choleskyRow_[k];
     916      assert (permuteInverse_[kRow]<0);
     917      //if (permuteInverse_[kRow]<0)
     918      which[nSave++]=kRow;
     919    }
     920    for (int i=0;i<nSave;i++) {
     921      int kRow = which[i];
     922      int length = permute_[kRow];
     923      CoinBigIndex start=choleskyStart_[kRow];
     924      CoinBigIndex end=choleskyStart_[kRow]+length;
     925      for (CoinBigIndex j=start;j<end;j++) {
     926        int jRow=choleskyRow_[j];
     927        mark[jRow]=1;
     928      }
     929      mark[kRow]=1;
     930      int n=nSave;
     931      for (int j=0;j<nSave;j++) {
     932        int jRow=which[j];
     933        if (!mark[jRow]) {
     934          which[n++]=jRow;
     935        }
     936      }
     937      for (CoinBigIndex j=start;j<end;j++) {
     938        int jRow=choleskyRow_[j];
     939        mark[jRow]=0;
     940      }
     941      mark[kRow]=0;
     942      if (n>nSave) {
     943        bool inPlace = (n-nSave==1);
     944        if (!inPlace) {
     945          // extend
     946          int length = n-nSave+end-start;
     947          if (length+lastSpace>space) {
     948            // need to compress
     949            numberCompressions++;
     950            int newN=0;
     951            for (int iRow=0;iRow<numberRows_;iRow++) {
     952              CoinBigIndex start=choleskyStart_[iRow];
     953              if (permuteInverse_[iRow]<0) {
     954                sort[newN]=start;
     955                order[newN++]=iRow;
     956              } else {
     957                choleskyStart_[iRow]=-1;
     958                permute_[iRow]=0;
     959              }
     960            }
     961            CoinSort_2(sort,sort+newN,order);
     962            sizeFactor=0;
     963            for (int k=0;k<newN;k++) {
     964              int iRow=order[k];
     965              int length = permute_[iRow];
     966              CoinBigIndex start=choleskyStart_[iRow];
     967              choleskyStart_[iRow]=sizeFactor;
     968              for (int j=0;j<length;j++)
     969                choleskyRow_[sizeFactor+j]=choleskyRow_[start+j];
     970              sizeFactor += length;
     971            }
     972            lastSpace=sizeFactor;
     973            if ((sizeFactor+numberRows_+1000)*4>3*space) {
     974              // need to expand
     975              numberExpansions++;
     976              space = (3*space)/2;
     977              int * temp = new int [space];
     978              memcpy(temp,choleskyRow_,sizeFactor*sizeof(int));
     979              delete [] choleskyRow_;
     980              choleskyRow_=temp;
     981            }
     982          }
     983        }
     984        // Now add
     985        start=choleskyStart_[kRow];
     986        end=choleskyStart_[kRow]+permute_[kRow];
     987        if (!inPlace)
     988          choleskyStart_[kRow]=lastSpace;
     989        CoinBigIndex put = choleskyStart_[kRow];
     990        for (CoinBigIndex j=start;j<end;j++) {
     991          int jRow=choleskyRow_[j];
     992          if (permuteInverse_[jRow]<0)
     993            choleskyRow_[put++]=jRow;
     994        }
     995        for (int j=nSave;j<n;j++) {
     996          int jRow=which[j];
     997          choleskyRow_[put++]=jRow;
     998        }
     999        if (!inPlace) {
     1000          permute_[kRow]=put-lastSpace;
     1001          lastSpace=put;
     1002        } else {
     1003          permute_[kRow]=put-choleskyStart_[kRow];
     1004        }
     1005      } else {
     1006        // pack down for new counts
     1007        CoinBigIndex put=start;
     1008        for (CoinBigIndex j=start;j<end;j++) {
     1009          int jRow=choleskyRow_[j];
     1010          if (permuteInverse_[jRow]<0)
     1011            choleskyRow_[put++]=jRow;
     1012        }
     1013        permute_[kRow]=put-start;
     1014      }
     1015      // take out
     1016      int iNext = next[kRow];
     1017      int iPrevious = previous[kRow];
     1018      if (iPrevious<numberRows_) {
     1019        next[iPrevious]=iNext;
     1020      } else {
     1021        assert (iPrevious==length+large);
     1022        first[length]=iNext;
     1023      }
     1024      if (iNext<numberRows_) {
     1025        previous[iNext]=iPrevious;
     1026      } else {
     1027        assert (iNext==length+2*large);
     1028      }
     1029      // put in
     1030      length=permute_[kRow];
     1031      smallest = CoinMin(smallest,length);
     1032      if (first[length]<0||first[length]>numberRows_) {
     1033        first[length]=kRow;
     1034        previous[kRow]=length+large;
     1035        next[kRow]=length+2*large;
     1036      } else {
     1037        int k=first[length];
     1038        assert (k<numberRows_);
     1039        first[length]=kRow;
     1040        previous[kRow]=length+large;
     1041        assert (previous[k]==length+large);
     1042        next[kRow]=k;
     1043        previous[k]=kRow;
     1044      }
     1045    }
     1046  }
     1047  // do rest
     1048  for (int iRow=0;iRow<numberRows_;iRow++) {
     1049    if (permuteInverse_[iRow]<0)
     1050      permuteInverse_[iRow]=done++;
     1051  }
     1052  printf("%d compressions, %d expansions\n",
     1053         numberCompressions,numberExpansions);
     1054  assert (done==numberRows_);
     1055  delete [] sort;
     1056  delete [] order;
     1057  delete [] which;
     1058  delete [] mark;
     1059  delete [] first;
     1060  delete [] next;
     1061  delete [] previous;
     1062  delete [] choleskyRow_;
     1063  choleskyRow_=NULL;
     1064  delete [] choleskyStart_;
     1065  choleskyStart_=NULL;
     1066#ifndef NDEBUG
     1067  for (int iRow=0;iRow<numberRows_;iRow++) {
     1068    permute_[iRow]=-1;
     1069  }
     1070#endif
     1071  for (int iRow=0;iRow<numberRows_;iRow++) {
     1072    permute_[permuteInverse_[iRow]]=iRow;
     1073  }
     1074#ifndef NDEBUG
     1075  for (int iRow=0;iRow<numberRows_;iRow++) {
     1076    assert(permute_[iRow]>=0&&permute_[iRow]<numberRows_);
     1077  }
     1078#endif
     1079  return returnCode;
     1080}
     1081#elif BASE_ORDER==2
     1082/*----------------------------------------------------------------------------*/
     1083/*      (C) Copyright IBM Corporation 1997, 2009.  All Rights Reserved.       */
     1084/*----------------------------------------------------------------------------*/
     1085
     1086/*  A compact no-frills Approximate Minimum Local Fill ordering code.
     1087
     1088    References:
     1089
     1090[1] Ordering Sparse Matrices Using Approximate Minimum Local Fill.
     1091    Edward Rothberg, SGI Manuscript, April 1996.
     1092[2] An Approximate Minimum Degree Ordering Algorithm.
     1093    T. Davis, P. Amestoy, and I. Duff, TR-94-039, CIS Department,
     1094    University of Florida, December 1994.
     1095*/
     1096
     1097#include <math.h>
     1098#include <stdlib.h>
     1099
     1100typedef int             WSI;
     1101
     1102#define NORDTHRESH      7
     1103#define MAXIW           2147000000
     1104
     1105//#define WSSMP_DBG
     1106#ifdef WSSMP_DBG
     1107static void chk (WSI ind, WSI i, WSI lim) {
     1108
     1109  if (i <= 0 || i > lim) {
     1110    printf("%d: chk: myamlf: WAR**: i, lim = %d, %d\n",ind,i,lim);
     1111  }
     1112}
     1113#endif
     1114
     1115static void
     1116myamlf(WSI n, WSI xadj[], WSI adjncy[], WSI dgree[], WSI varbl[],
     1117       WSI snxt[], WSI perm[], WSI invp[], WSI head[], WSI lsize[],
     1118       WSI flag[], WSI erscore[], WSI locaux, WSI adjln, WSI speed)
     1119{
     1120WSI l, i, j, k;
     1121double x, y;
     1122WSI maxmum, fltag, nodeg, scln, nm1, deg, tn,
     1123    locatns, ipp, jpp, nnode, numpiv, node,
     1124    nodeln, currloc, counter, numii, mindeg,
     1125    i0, i1, i2, i4, i5, i6, i7, i9,
     1126    j0, j1, j2, j3, j4, j5, j6, j7, j8, j9;
     1127/*                                             n cannot be less than NORDTHRESH
     1128 if (n < 3) {
     1129    if (n > 0) perm[0] = invp[0] = 1;
     1130    if (n > 1) perm[1] = invp[1] = 2;
     1131    return;
     1132 }
     1133*/
     1134#ifdef WSSMP_DBG
     1135 printf("Myamlf: n, locaux, adjln, speed = %d,%d,%d,%d\n",n,locaux,adjln,speed);
     1136 for (i = 0; i < n; i++) flag[i] = 0;
     1137 k = xadj[0];
     1138 for (i = 1; i <= n; i++) {
     1139   j = k + dgree[i-1];
     1140   if (j < k || k < 1) printf("ERR**: myamlf: i, j, k = %d,%d,%d\n",i,j,k);
     1141   for (l = k; l < j; l++) {
     1142      if (adjncy[l - 1] < 1 || adjncy[l - 1] > n || adjncy[l - 1] == i)
     1143        printf("ERR**: myamlf: i, l, adjj[l] = %d,%d,%d\n",i,l,adjncy[l - 1]);
     1144      if (flag[adjncy[l - 1] - 1] == i)
     1145        printf("ERR**: myamlf: duplicate entry: %d - %d\n",i,adjncy[l - 1]);
     1146      flag[adjncy[l - 1] - 1] = i;
     1147      nm1 = adjncy[l - 1] - 1;
     1148      for (fltag = xadj[nm1]; fltag < xadj[nm1] + dgree[nm1]; fltag++) {
     1149        if (adjncy[fltag - 1] == i) goto L100;
     1150      }
     1151      printf("ERR**: Unsymmetric %d %d\n",i,fltag);
     1152L100:;
     1153   }
     1154   k = j;
     1155   if (k > locaux) printf("ERR**: i, j, k, locaux = %d, %d, %d, %d\n",
     1156                                  i, j, k, locaux);
     1157 }
     1158 if (n*(n-1) < locaux-1) printf("ERR**: myamlf: too many nnz, n, ne = %d, %d\n",
     1159                                n,adjln);
     1160 for (i = 0; i < n; i++) flag[i] = 1;
     1161 if (n > 10000) printf ("Finished error checking in AMF\n");
     1162 for (i = locaux; i <= adjln; i++) adjncy[i - 1] = -22;
     1163#endif
     1164
     1165 maxmum = 0;
     1166 mindeg = 1;
     1167 fltag = 2;
     1168 locatns = locaux - 1;
     1169 nm1 = n - 1;
     1170 counter = 1;
     1171 for (l = 0; l < n; l++) {
     1172     j = erscore[l];
     1173#ifdef WSSMP_DBG
     1174chk(1,j,n);
     1175#endif
     1176     if (j > 0) {
     1177         nnode = head[j - 1];
     1178         if (nnode) perm[nnode - 1] = l + 1;
     1179         snxt[l] = nnode;
     1180         head[j - 1] = l + 1;
     1181     } else {
     1182         invp[l] = -(counter++);
     1183         flag[l] = xadj[l] = 0;
     1184     }
     1185 }
     1186 while (counter <= n) {
     1187     for (deg = mindeg; head[deg - 1] < 1; deg++);
     1188     nodeg = 0;
     1189#ifdef WSSMP_DBG
     1190chk(2,deg,n-1);
     1191#endif
     1192     node = head[deg - 1];
     1193#ifdef WSSMP_DBG
     1194chk(3,node,n);
     1195#endif
     1196     mindeg = deg;
     1197     nnode = snxt[node - 1];
     1198     if (nnode) perm[nnode - 1] = 0;
     1199     head[deg - 1] = nnode;
     1200     nodeln = invp[node - 1];
     1201     numpiv = varbl[node - 1];
     1202     invp[node - 1] = -counter;
     1203     counter += numpiv;
     1204     varbl[node - 1] = -numpiv;
     1205     if (nodeln) {
     1206         j4 = locaux;
     1207         i5 = lsize[node - 1] - nodeln;
     1208         i2 = nodeln + 1;
     1209         l = xadj[node - 1];
     1210         for (i6 = 1; i6 <= i2; ++i6) {
     1211             if (i6 == i2) {
     1212                 tn = node, i0 = l, scln = i5;
     1213             } else {
     1214#ifdef WSSMP_DBG
     1215chk(4,l,adjln);
     1216#endif
     1217                 tn = adjncy[l-1];
     1218                 l++;
     1219#ifdef WSSMP_DBG
     1220chk(5,tn,n);
     1221#endif
     1222                 i0 = xadj[tn - 1];
     1223                 scln = lsize[tn - 1];
     1224             }
     1225             for (i7 = 1; i7 <= scln; ++i7) {
     1226#ifdef WSSMP_DBG
     1227chk(6,i0,adjln);
     1228#endif
     1229               i = adjncy[i0 - 1];
     1230               i0++;
     1231#ifdef WSSMP_DBG
     1232chk(7,i,n);
     1233#endif
     1234               numii = varbl[i - 1];
     1235               if (numii > 0) {
     1236                 if (locaux > adjln) {
     1237                   lsize[node - 1] -= i6;
     1238                   xadj[node - 1] = l;
     1239                   if (!lsize[node - 1]) xadj[node - 1] = 0;
     1240                   xadj[tn - 1] = i0;
     1241                   lsize[tn - 1] = scln - i7;
     1242                   if (!lsize[tn - 1]) xadj[tn - 1] = 0;
     1243                   for (j = 1; j <= n; j++) {
     1244                     i4 = xadj[j - 1];
     1245                     if (i4 > 0) {
     1246                       xadj[j - 1] = adjncy[i4 - 1];
     1247#ifdef WSSMP_DBG
     1248chk(8,i4,adjln);
     1249#endif
     1250                       adjncy[i4 - 1] = -j;
     1251                     }
     1252                   }
     1253                   i9 = j4 - 1;
     1254                   j6 = j7 = 1;
     1255                   while (j6 <= i9) {
     1256#ifdef WSSMP_DBG
     1257chk(9,j6,adjln);
     1258#endif
     1259                     j = -adjncy[j6 - 1];
     1260                     j6++;
     1261                     if (j > 0) {
     1262#ifdef WSSMP_DBG
     1263chk(10,j7,adjln);
     1264chk(11,j,n);
     1265#endif
     1266                       adjncy[j7 - 1] = xadj[j - 1];
     1267                       xadj[j - 1] = j7++;
     1268                       j8 = lsize[j - 1] - 1 + j7;
     1269                       for (; j7 < j8; j7++,j6++) adjncy[j7-1] = adjncy[j6-1];
     1270                     }
     1271                   }
     1272                   for (j0=j7; j4 < locaux; j4++,j7++) {
     1273#ifdef WSSMP_DBG
     1274chk(12,j4,adjln);
     1275#endif
     1276                     adjncy[j7 - 1] = adjncy[j4 - 1];
     1277                   }
     1278                   j4 = j0;
     1279                   locaux = j7;
     1280                   i0 = xadj[tn - 1];
     1281                   l = xadj[node - 1];
     1282                 }
     1283                 adjncy[locaux-1] = i;
     1284                 locaux++;
     1285                 varbl[i - 1] = -numii;
     1286                 nodeg += numii;
     1287                 ipp = perm[i - 1];
     1288                 nnode = snxt[i - 1];
     1289#ifdef WSSMP_DBG
     1290if (ipp) chk(13,ipp,n);
     1291if (nnode) chk(14,nnode,n);
     1292#endif
     1293                 if (ipp) snxt[ipp - 1] = nnode;
     1294                 else head[erscore[i - 1] - 1] = nnode;
     1295                 if (nnode) perm[nnode - 1] = ipp;
     1296               }
     1297             }
     1298             if (tn != node) {
     1299                 flag[tn - 1] = 0, xadj[tn - 1] = -node;
     1300             }
     1301         }
     1302         currloc = (j5 = locaux) - j4;
     1303         locatns += currloc;
     1304     } else {
     1305         i1 = (j4 = xadj[node - 1]) + lsize[node - 1];
     1306         for (j = j5 = j4; j < i1; j++) {
     1307#ifdef WSSMP_DBG
     1308chk(15,j,adjln);
     1309#endif
     1310             i = adjncy[j - 1];
     1311#ifdef WSSMP_DBG
     1312chk(16,i,n);
     1313#endif
     1314             numii = varbl[i - 1];
     1315             if (numii > 0) {
     1316                 nodeg += numii;
     1317                 varbl[i - 1] = -numii;
     1318#ifdef WSSMP_DBG
     1319chk(17,j5,adjln);
     1320#endif
     1321                 adjncy[j5-1] = i;
     1322                 ipp = perm[i - 1];
     1323                 nnode = snxt[i - 1];
     1324                 j5++;
     1325#ifdef WSSMP_DBG
     1326if (ipp) chk(18,ipp,n);
     1327if (nnode) chk(19,nnode,n);
     1328#endif
     1329                 if (ipp) snxt[ipp - 1] = nnode;
     1330                 else head[erscore[i - 1] - 1] = nnode;
     1331                 if (nnode) perm[nnode - 1] = ipp;
     1332             }
     1333         }
     1334         currloc = 0;
     1335     }
     1336#ifdef WSSMP_DBG
     1337chk(20,node,n);
     1338#endif
     1339     lsize[node - 1] = j5 - (xadj[node - 1] = j4);
     1340     dgree[node - 1] = nodeg;
     1341     if (fltag+n < 0 || fltag+n > MAXIW) {
     1342       for (i = 1; i <= n; i++) if (flag[i - 1]) flag[i - 1] = 1;
     1343       fltag = 2;
     1344     }
     1345     for (j3 = j4; j3 < j5; j3++) {
     1346#ifdef WSSMP_DBG
     1347chk(21,j3,adjln);
     1348#endif
     1349       i = adjncy[j3 - 1];
     1350#ifdef WSSMP_DBG
     1351chk(22,i,n);
     1352#endif
     1353       j = invp[i - 1];
     1354       if (j > 0) {
     1355         i4 = fltag - (numii = -varbl[i - 1]);
     1356         i2 = xadj[i - 1] + j;
     1357         for (l = xadj[i - 1]; l < i2; l++) {
     1358#ifdef WSSMP_DBG
     1359chk(23,l,adjln);
     1360#endif
     1361           tn = adjncy[l - 1];
     1362#ifdef WSSMP_DBG
     1363chk(24,tn,n);
     1364#endif
     1365           j9 = flag[tn - 1];
     1366           if (j9 >= fltag) j9 -= numii; else if (j9) j9 = dgree[tn - 1] + i4;
     1367           flag[tn - 1] = j9;
     1368         }
     1369       }
     1370     }
     1371     for (j3 = j4; j3 < j5; j3++) {
     1372#ifdef WSSMP_DBG
     1373chk(25,j3,adjln);
     1374#endif
     1375         i = adjncy[j3 - 1];
     1376         i5 = deg = 0;
     1377#ifdef WSSMP_DBG
     1378chk(26,i,n);
     1379#endif
     1380         j1 = (i4 = xadj[i - 1]) + invp[i - 1];
     1381         for (l = j0 = i4; l < j1; l++) {
     1382#ifdef WSSMP_DBG
     1383chk(27,l,adjln);
     1384#endif
     1385             tn = adjncy[l - 1];
     1386#ifdef WSSMP_DBG
     1387chk(70,tn,n);
     1388#endif
     1389             j8 = flag[tn - 1];
     1390             if (j8) {
     1391                 deg += (j8 - fltag);
     1392#ifdef WSSMP_DBG
     1393chk(28,i4,adjln);
     1394#endif
     1395                 adjncy[i4-1] = tn;
     1396                 i5 += tn;
     1397                 i4++;
     1398                 while (i5 >= nm1) i5 -= nm1;
     1399             }
     1400         }
     1401         invp[i - 1] = (j2 = i4) - j0 + 1;
     1402         i2 = j0 + lsize[i - 1];
     1403         for (l = j1; l < i2; l++) {
     1404#ifdef WSSMP_DBG
     1405chk(29,l,adjln);
     1406#endif
     1407             j = adjncy[l - 1];
     1408#ifdef WSSMP_DBG
     1409chk(30,j,n);
     1410#endif
     1411             numii = varbl[j - 1];
     1412             if (numii > 0) {
     1413                 deg += numii;
     1414#ifdef WSSMP_DBG
     1415chk(31,i4,adjln);
     1416#endif
     1417                 adjncy[i4-1] = j;
     1418                 i5 += j;
     1419                 i4++;
     1420                 while (i5 >= nm1) i5 -= nm1;
     1421             }
     1422         }
     1423         if (invp[i - 1] == 1 && j2 == i4) {
     1424             numii = -varbl[i - 1];
     1425             xadj[i - 1] = -node;
     1426             nodeg -= numii;
     1427             counter += numii;
     1428             numpiv += numii;
     1429             invp[i - 1] = varbl[i - 1] = 0;
     1430         } else {
     1431             if (dgree[i - 1] > deg) dgree[i - 1] = deg;
     1432#ifdef WSSMP_DBG
     1433chk(32,j2,adjln);
     1434chk(33,j0,adjln);
     1435#endif
     1436             adjncy[i4 - 1] = adjncy[j2 - 1];
     1437             adjncy[j2 - 1] = adjncy[j0 - 1];
     1438             adjncy[j0 - 1] = node;
     1439             lsize[i - 1] = i4 - j0 + 1;
     1440             i5++;
     1441#ifdef WSSMP_DBG
     1442chk(35,i5,n);
     1443#endif
     1444             j = head[i5 - 1];
     1445             if (j > 0) {
     1446                 snxt[i - 1] = perm[j - 1];
     1447                 perm[j - 1] = i;
     1448             } else {
     1449                 snxt[i - 1] = -j;
     1450                 head[i5 - 1] = -i;
     1451             }
     1452             perm[i - 1] = i5;
     1453         }
     1454     }
     1455#ifdef WSSMP_DBG
     1456chk(36,node,n);
     1457#endif
     1458     dgree[node - 1] = nodeg;
     1459     if (maxmum < nodeg) maxmum = nodeg;
     1460     fltag += maxmum;
     1461#ifdef WSSMP_DBG
     1462     if (fltag+n+1 < 0) printf("ERR**: fltag = %d (A)\n",fltag);
     1463#endif
     1464     for (j3 = j4; j3 < j5; ++j3) {
     1465#ifdef WSSMP_DBG
     1466chk(37,j3,adjln);
     1467#endif
     1468       i = adjncy[j3 - 1];
     1469#ifdef WSSMP_DBG
     1470chk(38,i,n);
     1471#endif
     1472       if (varbl[i - 1] < 0) {
     1473         i5 = perm[i - 1];
     1474#ifdef WSSMP_DBG
     1475chk(39,i5,n);
     1476#endif
     1477         j = head[i5 - 1];
     1478         if (j) {
     1479           if (j < 0) {
     1480                head[i5 - 1] = 0, i = -j;
     1481           } else {
     1482#ifdef WSSMP_DBG
     1483chk(40,j,n);
     1484#endif
     1485                i = perm[j - 1];
     1486                perm[j - 1] = 0;
     1487           }
     1488           while (i) {
     1489#ifdef WSSMP_DBG
     1490chk(41,i,n);
     1491#endif
     1492             if (!snxt[i - 1]) {
     1493                 i = 0;
     1494             } else {
     1495               k = invp[i - 1];
     1496               i2 = xadj[i - 1] + (scln = lsize[i - 1]);
     1497               for (l = xadj[i - 1] + 1; l < i2; l++) {
     1498#ifdef WSSMP_DBG
     1499chk(42,l,adjln);
     1500chk(43,adjncy[l - 1],n);
     1501#endif
     1502                 flag[adjncy[l - 1] - 1] = fltag;
     1503               }
     1504               jpp = i;
     1505               j = snxt[i - 1];
     1506               while (j) {
     1507#ifdef WSSMP_DBG
     1508chk(44,j,n);
     1509#endif
     1510                 if (lsize[j - 1] == scln && invp[j - 1] == k) {
     1511                   i2 = xadj[j - 1] + scln;
     1512                   j8 = 1;
     1513                   for (l = xadj[j - 1] + 1; l < i2; l++) {
     1514#ifdef WSSMP_DBG
     1515chk(45,l,adjln);
     1516chk(46,adjncy[l - 1],n);
     1517#endif
     1518                       if (flag[adjncy[l - 1] - 1] != fltag) {
     1519                         j8 = 0;
     1520                         break;
     1521                       }
     1522                   }
     1523                   if (j8) {
     1524                     xadj[j - 1] = -i;
     1525                     varbl[i - 1] += varbl[j - 1];
     1526                     varbl[j - 1] = invp[j - 1] = 0;
     1527#ifdef WSSMP_DBG
     1528chk(47,j,n);
     1529chk(48,jpp,n);
     1530#endif
     1531                     j = snxt[j - 1];
     1532                     snxt[jpp - 1] = j;
     1533                   } else {
     1534                     jpp = j;
     1535#ifdef WSSMP_DBG
     1536chk(49,j,n);
     1537#endif
     1538                     j = snxt[j - 1];
     1539                   }
     1540                 } else {
     1541                   jpp = j;
     1542#ifdef WSSMP_DBG
     1543chk(50,j,n);
     1544#endif
     1545                   j = snxt[j - 1];
     1546                 }
     1547               }
     1548               fltag++;
     1549#ifdef WSSMP_DBG
     1550               if (fltag+n+1 < 0) printf("ERR**: fltag = %d (B)\n",fltag);
     1551#endif
     1552#ifdef WSSMP_DBG
     1553chk(51,i,n);
     1554#endif
     1555               i = snxt[i - 1];
     1556             }
     1557           }
     1558         }
     1559       }
     1560     }
     1561     j8 = nm1 - counter;
     1562     switch (speed) {
     1563case 1:
     1564       for (j = j3 = j4; j3 < j5; j3++) {
     1565#ifdef WSSMP_DBG
     1566chk(52,j3,adjln);
     1567#endif
     1568         i = adjncy[j3 - 1];
     1569#ifdef WSSMP_DBG
     1570chk(53,i,n);
     1571#endif
     1572         numii = varbl[i - 1];
     1573         if (numii < 0) {
     1574           k = 0;
     1575           i4 = (l = xadj[i - 1]) + invp[i - 1];
     1576           for (; l < i4; l++) {
     1577#ifdef WSSMP_DBG
     1578chk(54,l,adjln);
     1579chk(55,adjncy[l - 1],n);
     1580#endif
     1581              i5 = dgree[adjncy[l - 1] - 1];
     1582              if (k < i5) k = i5;
     1583           }
     1584           x = static_cast<double> (k - 1);
     1585           varbl[i - 1] = abs(numii);
     1586           j9 = dgree[i - 1] + nodeg;
     1587           deg = (j8 > j9 ? j9 : j8) + numii;
     1588           if (deg < 1) deg = 1;
     1589           y = static_cast<double> (deg);
     1590           dgree[i - 1] = deg;
     1591/*
     1592           printf("%i %f should >= %i %f\n",deg,y,k-1,x);
     1593           if (y < x) printf("** \n"); else printf("\n");
     1594*/
     1595           y = y*y - y;
     1596           x = y - (x*x - x);
     1597           if (x < 1.1) x = 1.1;
     1598           deg = static_cast<WSI>(sqrt(x));
     1599/*         if (deg < 1) deg = 1; */
     1600           erscore[i - 1] = deg;
     1601#ifdef WSSMP_DBG
     1602chk(56,deg,n-1);
     1603#endif
     1604           nnode = head[deg - 1];
     1605           if (nnode) perm[nnode - 1] = i;
     1606           snxt[i - 1] = nnode;
     1607           perm[i - 1] = 0;
     1608#ifdef WSSMP_DBG
     1609chk(57,j,adjln);
     1610#endif
     1611           head[deg - 1] = adjncy[j-1] = i;
     1612           j++;
     1613           if (deg < mindeg) mindeg = deg;
     1614         }
     1615       }
     1616       break;
     1617case 2:
     1618       for (j = j3 = j4; j3 < j5; j3++) {
     1619#ifdef WSSMP_DBG
     1620chk(58,j3,adjln);
     1621#endif
     1622         i = adjncy[j3 - 1];
     1623#ifdef WSSMP_DBG
     1624chk(59,i,n);
     1625#endif
     1626         numii = varbl[i - 1];
     1627         if (numii < 0) {
     1628           x = static_cast<double> (dgree[adjncy[xadj[i - 1] - 1] - 1] - 1);
     1629           varbl[i - 1] = abs(numii);
     1630           j9 = dgree[i - 1] + nodeg;
     1631           deg = (j8 > j9 ? j9 : j8) + numii;
     1632           if (deg < 1) deg = 1;
     1633           y = static_cast<double> (deg);
     1634           dgree[i - 1] = deg;
     1635/*
     1636           printf("%i %f should >= %i %f",deg,y,dgree[adjncy[xadj[i - 1] - 1] - 1]-1,x);
     1637           if (y < x) printf("** \n"); else printf("\n");
     1638*/
     1639           y = y*y - y;
     1640           x = y - (x*x - x);
     1641           if (x < 1.1) x = 1.1;
     1642           deg = static_cast<WSI>(sqrt(x));
     1643/*         if (deg < 1) deg = 1; */
     1644           erscore[i - 1] = deg;
     1645#ifdef WSSMP_DBG
     1646chk(60,deg,n-1);
     1647#endif
     1648           nnode = head[deg - 1];
     1649           if (nnode) perm[nnode - 1] = i;
     1650           snxt[i - 1] = nnode;
     1651           perm[i - 1] = 0;
     1652#ifdef WSSMP_DBG
     1653chk(61,j,adjln);
     1654#endif
     1655           head[deg - 1] = adjncy[j-1] = i;
     1656           j++;
     1657           if (deg < mindeg) mindeg = deg;
     1658         }
     1659       }
     1660       break;
     1661default:
     1662       for (j = j3 = j4; j3 < j5; j3++) {
     1663#ifdef WSSMP_DBG
     1664chk(62,j3,adjln);
     1665#endif
     1666         i = adjncy[j3 - 1];
     1667#ifdef WSSMP_DBG
     1668chk(63,i,n);
     1669#endif
     1670         numii = varbl[i - 1];
     1671         if (numii < 0) {
     1672           varbl[i - 1] = abs(numii);
     1673           j9 = dgree[i - 1] + nodeg;
     1674           deg = (j8 > j9 ? j9 : j8) + numii;
     1675           if (deg < 1) deg = 1;
     1676           dgree[i - 1] = deg;
     1677#ifdef WSSMP_DBG
     1678chk(64,deg,n-1);
     1679#endif
     1680           nnode = head[deg - 1];
     1681           if (nnode) perm[nnode - 1] = i;
     1682           snxt[i - 1] = nnode;
     1683           perm[i - 1] = 0;
     1684#ifdef WSSMP_DBG
     1685chk(65,j,adjln);
     1686#endif
     1687           head[deg - 1] = adjncy[j-1] = i;
     1688           j++;
     1689           if (deg < mindeg) mindeg = deg;
     1690         }
     1691       }
     1692       break;
     1693     }
     1694     if (currloc) {
     1695#ifdef WSSMP_DBG
     1696chk(66,node,n);
     1697#endif
     1698       locatns += (lsize[node - 1] - currloc), locaux = j;
     1699     }
     1700     varbl[node - 1] = numpiv + nodeg;
     1701     lsize[node - 1] = j - j4;
     1702     if (!lsize[node - 1]) flag[node - 1] = xadj[node - 1] = 0;
     1703 }
     1704 for (l = 1; l <= n; l++) {
     1705   if (!invp[l - 1]) {
     1706     for (i = -xadj[l - 1]; invp[i - 1] >= 0; i = -xadj[i - 1]);
     1707       tn = i;
     1708#ifdef WSSMP_DBG
     1709chk(67,tn,n);
     1710#endif
     1711       k = -invp[tn - 1];
     1712       i = l;
     1713#ifdef WSSMP_DBG
     1714chk(68,i,n);
     1715#endif
     1716       while (invp[i - 1] >= 0) {
     1717         nnode = -xadj[i - 1];
     1718         xadj[i - 1] = -tn;
     1719         if (!invp[i - 1]) invp[i - 1] = k++;
     1720         i = nnode;
     1721       }
     1722       invp[tn - 1] = -k;
     1723   }
     1724 }
     1725 for (l = 0; l < n; l++) {
     1726   i = abs(invp[l]);
     1727#ifdef WSSMP_DBG
     1728chk(69,i,n);
     1729#endif
     1730   invp[l] = i;
     1731   perm[i - 1] = l + 1;
     1732 }
     1733 return;
     1734}
     1735// This code is not needed, but left in in case needed sometime
     1736#if 0
     1737/*C--------------------------------------------------------------------------*/
     1738void amlfdr(WSI *n, WSI xadj[], WSI adjncy[], WSI dgree[], WSI *adjln,
     1739            WSI *locaux, WSI varbl[], WSI snxt[], WSI perm[],
     1740            WSI head[], WSI invp[], WSI lsize[], WSI flag[], WSI *ispeed)
     1741{
     1742WSI nn,nlocaux,nadjln,speed,i,j,mx,mxj,*erscore;
     1743
     1744#ifdef WSSMP_DBG
     1745  printf("Calling amlfdr for n, speed = %d, %d\n",*n,*ispeed);
     1746#endif
     1747
     1748  if ((nn = *n) == 0) return;
     1749
     1750#ifdef WSSMP_DBG
     1751  if (nn == 31) {
     1752    printf("n = %d; adjln = %d; locaux = %d; ispeed = %d\n",
     1753            *n,*adjln,*locaux,*ispeed);
     1754  }
     1755#endif
     1756
     1757  if (nn < NORDTHRESH) {
     1758    for (i = 0; i < nn; i++) lsize[i] = i;
     1759    for (i = nn; i > 0; i--) {
     1760      mx = dgree[0];
     1761      mxj = 0;
     1762      for (j = 1; j < i; j++)
     1763        if (dgree[j] > mx) {
     1764          mx = dgree[j];
     1765          mxj = j;
     1766        }
     1767      invp[lsize[mxj]] = i;
     1768      dgree[mxj] = dgree[i-1];
     1769      lsize[mxj] = lsize[i-1];
     1770    }
     1771    return;
     1772  }
     1773  speed = *ispeed;
     1774  if (speed < 3) {
     1775/* 
     1776    erscore = (WSI *)malloc(nn * sizeof(WSI));
     1777    if (erscore == NULL) speed = 3;
     1778*/
     1779    wscmal (&nn, &i, &erscore);
     1780    if (i != 0) speed = 3;
     1781  }
     1782  if (speed > 2) erscore = dgree;
     1783  if (speed < 3) {
     1784    for (i = 0; i < nn; i++) {
     1785      perm[i] = 0;
     1786      invp[i] = 0;
     1787      head[i] = 0;
     1788      flag[i] = 1;
     1789      varbl[i] = 1;
     1790      lsize[i] = dgree[i];
     1791      erscore[i] = dgree[i];
     1792    }
     1793  } else {
     1794    for (i = 0; i < nn; i++) {
     1795      perm[i] = 0;
     1796      invp[i] = 0;
     1797      head[i] = 0;
     1798      flag[i] = 1;
     1799      varbl[i] = 1;
     1800      lsize[i] = dgree[i];
     1801    }
     1802  }
     1803  nlocaux = *locaux;
     1804  nadjln = *adjln;
     1805
     1806  myamlf(nn, xadj, adjncy, dgree, varbl, snxt, perm, invp,
     1807         head, lsize, flag, erscore, nlocaux, nadjln, speed);
     1808/*
     1809  if (speed < 3) free(erscore);
     1810*/
     1811  if (speed < 3) wscfr(&erscore);
     1812  return;
     1813}
     1814#endif // end of taking out amlfdr
     1815/*C--------------------------------------------------------------------------*/
     1816#endif
     1817// Orders rows
     1818int
     1819ClpCholeskyBase::orderAMD()
     1820{
     1821  permuteInverse_ = new int [numberRows_];
     1822  permute_ = new int[numberRows_];
     1823  // Do ordering
     1824  int returnCode=0;
     1825  // get full matrix
     1826  int space = 2*sizeFactor_+10000+4*numberRows_;
     1827  int * temp = new int [space];
     1828  CoinBigIndex * count = new CoinBigIndex [numberRows_];
     1829  CoinBigIndex * tempStart = new CoinBigIndex [numberRows_+1];
     1830  memset(count,0,numberRows_*sizeof(int));
     1831  for (int iRow=0;iRow<numberRows_;iRow++) {
     1832    count[iRow]+= choleskyStart_[iRow+1]-choleskyStart_[iRow]-1;
     1833    for (CoinBigIndex j=choleskyStart_[iRow]+1;j<choleskyStart_[iRow+1];j++) {
     1834      int jRow=choleskyRow_[j];
     1835      count[jRow]++;
     1836    }
     1837  }
     1838#define OFFSET 1
     1839  CoinBigIndex sizeFactor =0;
     1840  for (int iRow=0;iRow<numberRows_;iRow++) {
     1841    int length = count[iRow];
     1842    permute_[iRow]=length;
     1843    // add 1 to starts
     1844    tempStart[iRow]=sizeFactor+OFFSET;
     1845    count[iRow]=sizeFactor;
     1846    sizeFactor += length;
     1847  }
     1848  tempStart[numberRows_]=sizeFactor+OFFSET;
     1849  // add 1 to rows
     1850  for (int iRow=0;iRow<numberRows_;iRow++) {
     1851    assert (choleskyRow_[choleskyStart_[iRow]]==iRow);
     1852    for (CoinBigIndex j=choleskyStart_[iRow]+1;j<choleskyStart_[iRow+1];j++) {
     1853      int jRow=choleskyRow_[j];
     1854      int put = count[iRow];
     1855      temp[put]=jRow+OFFSET;
     1856      count[iRow]++;
     1857      put = count[jRow];
     1858      temp[put]=iRow+OFFSET;
     1859      count[jRow]++;
     1860    }
     1861  }
     1862  for (int iRow=1;iRow<numberRows_;iRow++)
     1863    assert (count[iRow-1]==tempStart[iRow]-OFFSET);
     1864  delete [] choleskyRow_;
     1865  choleskyRow_=temp;
     1866  delete [] choleskyStart_;
     1867  choleskyStart_=tempStart;
     1868  int locaux=sizeFactor+OFFSET;
     1869  delete [] count;
     1870  int speed=integerParameters_[0];
     1871  if (speed<1||speed>2)
     1872    speed=3;
     1873  int * use = new int [((speed<3) ? 7 : 6)*numberRows_];
     1874  int * dgree = use;
     1875  int * varbl = dgree+numberRows_;
     1876  int * snxt = varbl+numberRows_;
     1877  int * head = snxt+numberRows_;
     1878  int * lsize = head+numberRows_;
     1879  int * flag = lsize+numberRows_;
     1880  int * erscore;
     1881  for (int i=0;i<numberRows_;i++) {
     1882    dgree[i]=choleskyStart_[i+1]-choleskyStart_[i];
     1883    head[i]=dgree[i];
     1884    snxt[i]=0;
     1885    permute_[i]=0;
     1886    permuteInverse_[i]=0;
     1887    head[i]=0;
     1888    flag[i]=1;
     1889    varbl[i]=1;
     1890    lsize[i]=dgree[i];
     1891  }
     1892  if (speed<3) {
     1893    erscore = flag+numberRows_;
     1894    for (int i=0;i<numberRows_;i++)
     1895      erscore[i]=dgree[i];
     1896  } else {
     1897    erscore = dgree;
     1898  }
     1899  myamlf(numberRows_, choleskyStart_, choleskyRow_,
     1900         dgree, varbl, snxt, permute_, permuteInverse_,
     1901         head, lsize, flag, erscore, locaux, space, speed);
     1902  for (int iRow=0;iRow<numberRows_;iRow++) {
     1903    permute_[iRow]--;
     1904  }
     1905  for (int iRow=0;iRow<numberRows_;iRow++) {
     1906    permuteInverse_[permute_[iRow]]=iRow;
     1907  }
     1908  for (int iRow=0;iRow<numberRows_;iRow++) {
     1909    assert (permuteInverse_[iRow]>=0&&permuteInverse_[iRow]<numberRows_);
     1910  }
     1911  delete [] use;
     1912  delete [] choleskyRow_;
     1913  choleskyRow_=NULL;
     1914  delete [] choleskyStart_;
     1915  choleskyStart_=NULL;
     1916  return returnCode;
    7591917}
    7601918/* Does Symbolic factorization given permutation.
     
    13382496      int k=iRow;
    13392497      int linked = link_[iRow];
     2498#ifndef NDEBUG
     2499      int count=0;
     2500#endif
    13402501      while (linked<=kRow) {
    13412502        k=linked;
    13422503        linked = link_[k];
     2504#ifndef NDEBUG
     2505        count++;
     2506        assert (count<numberRows_);
     2507#endif
    13432508      }
    13442509      nz++;
     
    23873552    dense.setIntegerParameter(20,0);
    23883553    dense.setIntegerParameter(34,firstPositive);
     3554    dense.setModel(model_);
    23893555    dense.factorizePart2(dropped);
    23903556    largest=dense.getDoubleParameter(3);
  • trunk/Clp/src/ClpCholeskyBase.hpp

    r1367 r1368  
    3030#if CLP_LONG_CHOLESKY>1
    3131typedef long double longDouble;
    32 #define CHOL_SMALL_VALUE 1.0e-15
     32#define CHOL_SMALL_VALUE 1.0e-15 
    3333#elif CLP_LONG_CHOLESKY==1
    3434typedef double longDouble;
     
    7878  virtual void solve (CoinWorkDouble * region) ;
    7979#endif
     80private:
     81  /// AMD ordering
     82  int orderAMD();
     83public:
    8084  //@}
    8185
     
    166170protected:
    167171  /// Sets type
    168   void setType(int type) {type_=type;}
     172  inline void setType(int type) {type_=type;}
     173  /// model.
     174  inline void setModel(ClpInterior * model)
     175  { model_=model;}
    169176   //@}
    170177   
  • trunk/Clp/src/ClpCholeskyDense.cpp

    r1367 r1368  
    1 // Copyright (C) 2003, International Business Machines
    2 // Corporation and others.  All Rights Reserved.
    3 
    4 
    5 
     1/* Copyright (C) 2003, International Business Machines Corporation
     2   and others.  All Rights Reserved. */
    63#include "CoinPragma.hpp"
    74#include "CoinHelperFunctions.hpp"
     
    1310#include "ClpQuadraticObjective.hpp"
    1411
    15 //#############################################################################
    16 // Constructors / Destructor / Assignment
    17 //#############################################################################
    18 
    19 //-------------------------------------------------------------------
    20 // Default Constructor
    21 //-------------------------------------------------------------------
     12/*#############################################################################*/
     13/* Constructors / Destructor / Assignment*/
     14/*#############################################################################*/
     15
     16/*-------------------------------------------------------------------*/
     17/* Default Constructor */
     18/*-------------------------------------------------------------------*/
    2219ClpCholeskyDense::ClpCholeskyDense ()
    2320  : ClpCholeskyBase(),
     
    2724}
    2825
    29 //-------------------------------------------------------------------
    30 // Copy constructor
    31 //-------------------------------------------------------------------
     26/*-------------------------------------------------------------------*/
     27/* Copy constructor */
     28/*-------------------------------------------------------------------*/
    3229ClpCholeskyDense::ClpCholeskyDense (const ClpCholeskyDense & rhs)
    3330  : ClpCholeskyBase(rhs),
    3431    borrowSpace_(rhs.borrowSpace_)
    3532{
    36   assert(!rhs.borrowSpace_||!rhs.sizeFactor_); // can't do if borrowing space
    37 }
    38 
    39 
    40 //-------------------------------------------------------------------
    41 // Destructor
    42 //-------------------------------------------------------------------
     33  assert(!rhs.borrowSpace_||!rhs.sizeFactor_); /* can't do if borrowing space*/
     34}
     35
     36
     37/*-------------------------------------------------------------------*/
     38/* Destructor */
     39/*-------------------------------------------------------------------*/
    4340ClpCholeskyDense::~ClpCholeskyDense ()
    4441{
    4542  if (borrowSpace_) {
    46     // set NULL
     43    /* set NULL*/
    4744    sparseFactor_=NULL;
    4845    workDouble_=NULL;
     
    5148}
    5249
    53 //----------------------------------------------------------------
    54 // Assignment operator
    55 //-------------------------------------------------------------------
     50/*----------------------------------------------------------------*/
     51/* Assignment operator */
     52/*-------------------------------------------------------------------*/
    5653ClpCholeskyDense &
    5754ClpCholeskyDense::operator=(const ClpCholeskyDense& rhs)
    5855{
    5956  if (this != &rhs) {
    60     assert(!rhs.borrowSpace_||!rhs.sizeFactor_); // can't do if borrowing space
     57    assert(!rhs.borrowSpace_||!rhs.sizeFactor_); /* can't do if borrowing space*/
    6158    ClpCholeskyBase::operator=(rhs);
    6259    borrowSpace_=rhs.borrowSpace_;
     
    6461  return *this;
    6562}
    66 //-------------------------------------------------------------------
    67 // Clone
    68 //-------------------------------------------------------------------
     63/*-------------------------------------------------------------------*/
     64/* Clone*/
     65/*-------------------------------------------------------------------*/
    6966ClpCholeskyBase * ClpCholeskyDense::clone() const
    7067{
    7168  return new ClpCholeskyDense(*this);
    7269}
    73 // If not power of 2 then need to redo a bit
     70/* If not power of 2 then need to redo a bit*/
    7471#define BLOCK 16
    7572#define BLOCKSHIFT 4
    76 // Block unroll if power of 2 and at least 8
     73/* Block unroll if power of 2 and at least 8*/
    7774#define BLOCKUNROLL
    7875
     
    8885  numberRows_ = numberRows;
    8986  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
    90   // allow one stripe extra
     87  /* allow one stripe extra*/
    9188  numberBlocks = numberBlocks + ((numberBlocks*(numberBlocks+1))/2);
    9289  sizeFactor_=numberBlocks*BLOCKSQ;
    93   //#define CHOL_COMPARE
     90  /*#define CHOL_COMPARE*/
    9491#ifdef CHOL_COMPARE 
    9592  sizeFactor_ += 95000;
     
    116113{
    117114 int numberBlocks = (numberRows+BLOCK-1)>>BLOCKSHIFT;
    118   // allow one stripe extra
     115 /* allow one stripe extra*/
    119116  numberBlocks = numberBlocks + ((numberBlocks*(numberBlocks+1))/2);
    120117  CoinBigIndex sizeFactor=numberBlocks*BLOCKSQ;
     
    165162  int numberColumns=model_->clpMatrix()->getNumCols();
    166163  CoinZeroN(sparseFactor_,sizeFactor_);
    167   //perturbation
     164  /*perturbation*/
    168165  CoinWorkDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
    169166  perturbation=perturbation*perturbation;
    170167  if (perturbation>1.0) {
    171168#ifdef COIN_DEVELOP
    172     //if (model_->model()->logLevel()&4)
     169    /*if (model_->model()->logLevel()&4) */
    173170      std::cout <<"large perturbation "<<perturbation<<std::endl;
    174171#endif
     
    180177  CoinWorkDouble largest=1.0;
    181178  CoinWorkDouble smallest=COIN_DBL_MAX;
    182   CoinWorkDouble delta2 = model_->delta(); // add delta*delta to diagonal
     179  CoinWorkDouble delta2 = model_->delta(); /* add delta*delta to diagonal*/
    183180  delta2 *= delta2;
    184181  if (!doKKT_) {
    185182    longDouble * work = sparseFactor_;
    186     work--; // skip diagonal
     183    work--; /* skip diagonal*/
    187184    int addOffset=numberRows_-1;
    188185    const double * diagonalSlack = diagonal + numberColumns;
    189     // largest in initial matrix
     186    /* largest in initial matrix*/
    190187    CoinWorkDouble largest2=1.0e-20;
    191188    for (iRow=0;iRow<numberRows_;iRow++) {
     
    217214        largest2 = CoinMax(largest2,CoinAbs(diagonalValue));
    218215      } else {
    219         // dropped
     216        /* dropped*/
    220217        diagonal_[iRow]=1.0;
    221218      }
     
    223220      work += addOffset;
    224221    }
    225     //check sizes
     222    /*check sizes*/
    226223    largest2*=1.0e-20;
    227224    largest = CoinMin(largest2,CHOL_SMALL_VALUE);
     
    229226    for (iRow=0;iRow<numberRows_;iRow++) {
    230227      int dropped=rowsDropped_[iRow];
    231       // Move to int array
     228      /* Move to int array*/
    232229      rowsDropped[iRow]=dropped;
    233230      if (!dropped) {
     
    246243    doubleParameters_[3]=0.0;
    247244    doubleParameters_[4]=COIN_DBL_MAX;
    248     integerParameters_[34]=0; // say all must be positive
     245    integerParameters_[34]=0; /* say all must be positive*/
    249246#ifdef CHOL_COMPARE 
    250247    if (numberRows_<200)
     
    258255      std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
    259256    choleskyCondition_=largest/smallest;
    260     //drop fresh makes some formADAT easier
     257    /*drop fresh makes some formADAT easier*/
    261258    if (newDropped||numberRowsDropped_) {
    262259      newDropped=0;
     
    265262        rowsDropped_[i]=dropped;
    266263        if (dropped==2) {
    267           //dropped this time
     264          /*dropped this time*/
    268265          rowsDropped[newDropped++]=i;
    269266          rowsDropped_[i]=0;
     
    274271    }
    275272  } else {
    276     // KKT
     273    /* KKT*/
    277274    CoinPackedMatrix * quadratic = NULL;
    278275    ClpQuadraticObjective * quadraticObj =
     
    280277    if (quadraticObj)
    281278      quadratic = quadraticObj->quadraticObjective();
    282     // matrix
     279    /* matrix*/
    283280    int numberRowsModel = model_->numberRows();
    284281    int numberColumns = model_->numberColumns();
    285282    int numberTotal = numberColumns + numberRowsModel;
    286283    longDouble * work = sparseFactor_;
    287     work--; // skip diagonal
     284    work--; /* skip diagonal*/
    288285    int addOffset=numberRows_-1;
    289286    int iColumn;
     
    298295          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
    299296          for (CoinBigIndex j=start;j<end;j++) {
    300             //choleskyRow_[numberElements]=row[j]+numberTotal;
    301             //sparseFactor_[numberElements++]=element[j];
     297            /*choleskyRow_[numberElements]=row[j]+numberTotal;*/
     298            /*sparseFactor_[numberElements++]=element[j];*/
    302299            work[row[j]+numberTotal]=element[j];
    303300            largest = CoinMax(largest,CoinAbs(element[j]));
     
    310307      }
    311308    } else {
    312       // Quadratic
     309      /* Quadratic*/
    313310      const int * columnQuadratic = quadratic->getIndices();
    314311      const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
     
    345342      }
    346343    }
    347     // slacks
     344    /* slacks*/
    348345    for (iColumn=numberColumns;iColumn<numberTotal;iColumn++) {
    349346      CoinWorkDouble value = diagonal[iColumn];
     
    359356      work += addOffset;
    360357    }
    361     // Finish diagonal
     358    /* Finish diagonal*/
    362359    for (iRow=0;iRow<numberRowsModel;iRow++) {
    363360      diagonal_[iRow+numberTotal]=delta2;
    364361    }
    365     //check sizes
     362    /*check sizes*/
    366363    largest*=1.0e-20;
    367364    largest = CoinMin(largest,CHOL_SMALL_VALUE);
     
    370367    doubleParameters_[3]=0.0;
    371368    doubleParameters_[4]=COIN_DBL_MAX;
    372     // Set up LDL cutoff
     369    /* Set up LDL cutoff*/
    373370    integerParameters_[34]=numberTotal;
    374     // KKT
     371    /* KKT*/
    375372    int * rowsDropped2 = new int[numberRows_];
    376373    CoinZeroN(rowsDropped2,numberRows_);
     
    386383      std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
    387384    choleskyCondition_=largest/smallest;
    388     // Should save adjustments in ..R_
     385    /* Should save adjustments in ..R_*/
    389386    int n1=0,n2=0;
    390387    CoinWorkDouble * primalR = model_->primalR();
     
    393390      if (rowsDropped2[iRow]) {
    394391        n1++;
    395         //printf("row region1 %d dropped\n",iRow);
    396         //rowsDropped_[iRow]=1;
     392        /*printf("row region1 %d dropped\n",iRow);*/
     393        /*rowsDropped_[iRow]=1;*/
    397394        rowsDropped_[iRow]=0;
    398395        primalR[iRow]=doubleParameters_[20];
     
    405402      if (rowsDropped2[iRow]) {
    406403        n2++;
    407         //printf("row region2 %d dropped\n",iRow);
    408         //rowsDropped_[iRow]=1;
     404        /*printf("row region2 %d dropped\n",iRow);*/
     405        /*rowsDropped_[iRow]=1;*/
    409406        rowsDropped_[iRow]=0;
    410407        dualR[iRow-numberTotal]=doubleParameters_[34];
     
    426423  diagonal_ = sparseFactor_ + 40000;
    427424  sparseFactor_=diagonal_ + numberRows_;
    428   //memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));
     425  /*memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));*/
    429426  CoinMemcpyN(xx,40000,sparseFactor_);
    430427  CoinMemcpyN(yy,numberRows_,diagonal_);
     
    435432  int firstPositive=integerParameters_[34];
    436433  longDouble * work = sparseFactor_;
    437   // Allow for triangular
     434  /* Allow for triangular*/
    438435  int addOffset=numberRows_-1;
    439436  work--;
     
    451448    bool dropColumn=false;
    452449    if (iColumn<firstPositive) {
    453       // must be negative
     450      /* must be negative*/
    454451      if (diagonalValue<=-dropValue) {
    455452        smallest = CoinMin(smallest,-diagonalValue);
     
    464461      }
    465462    } else {
    466       // must be positive
     463      /* must be positive*/
    467464      if (diagonalValue>=dropValue) {
    468465        smallest = CoinMin(smallest,diagonalValue);
     
    494491      }
    495492    } else {
    496       // drop column
     493      /* drop column*/
    497494      rowsDropped[iColumn]=2;
    498495      numberDropped++;
     
    511508  diagonal_=yy;
    512509}
    513 //#define POS_DEBUG
     510/*#define POS_DEBUG*/
    514511#ifdef POS_DEBUG
    515512static int counter=0;
     
    522519  iCol=0;
    523520  int kk=k>>BLOCKSQSHIFT;
    524   //printf("%d %d %d %d\n",k,kk,BLOCKSQ,BLOCKSQSHIFT);
     521  /*printf("%d %d %d %d\n",k,kk,BLOCKSQ,BLOCKSQSHIFT);*/
    525522  k=kk;
    526523  while (k>=numberBlocks) {
     
    541538{
    542539  int iColumn;
    543   //longDouble * xx = sparseFactor_;
    544   //longDouble * yy = diagonal_;
    545   //diagonal_ = sparseFactor_ + 40000;
    546   //sparseFactor_=diagonal_ + numberRows_;
    547   //memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));
    548   //memcpy(sparseFactor_,xx,40000*sizeof(double));
    549   //memcpy(diagonal_,yy,numberRows_*sizeof(double));
     540  /*longDouble * xx = sparseFactor_;*/
     541  /*longDouble * yy = diagonal_;*/
     542  /*diagonal_ = sparseFactor_ + 40000;*/
     543  /*sparseFactor_=diagonal_ + numberRows_;*/
     544  /*memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));*/
     545  /*memcpy(sparseFactor_,xx,40000*sizeof(double));*/
     546  /*memcpy(diagonal_,yy,numberRows_*sizeof(double));*/
    550547  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
    551   // later align on boundary
     548  /* later align on boundary*/
    552549  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
    553550  int n=numberRows_;
    554551  int nRound = numberRows_&(~(BLOCK-1));
    555   // adjust if exact
     552  /* adjust if exact*/
    556553  if (nRound==n)
    557554    nRound -= BLOCK;
    558555  int sizeLastBlock=n-nRound;
    559   int get = n*(n-1)/2; // ? as no diagonal
     556  int get = n*(n-1)/2; /* ? as no diagonal*/
    560557  int block = numberBlocks*(numberBlocks+1)/2;
    561558  int ifOdd;
     
    566563    ifOdd=1;
    567564    int put = BLOCKSQ;
    568     // do last separately
     565    /* do last separately*/
    569566    put -= (BLOCK-sizeLastBlock)*(BLOCK+1);
    570567    for (iColumn=numberRows_-1;iColumn>=nRound;iColumn--) {
     
    575572        assert (aa+put2>=sparseFactor_+get);
    576573      }
    577       // save diagonal as well
     574      /* save diagonal as well*/
    578575      aa[--put2]=diagonal_[iColumn];
    579576    }
     
    581578    block--;
    582579  } else {
    583     // exact fit
     580    /* exact fit*/
    584581    rowLast=numberRows_-1;
    585582    ifOdd=0;
    586583  }
    587   // Now main loop
     584  /* Now main loop*/
    588585  int nBlock=0;
    589586  for (;n>0;n-=BLOCK) {
     
    592589    int put = BLOCKSQ;
    593590    int putLast=0;
    594     // see if we have small block
     591    /* see if we have small block*/
    595592    if (ifOdd) {
    596593      aaLast = &a[(block-1)*BLOCKSQ];
     
    600597    for (iColumn=n-1;iColumn>=n-BLOCK;iColumn--) {
    601598      if (aaLast) {
    602         // last bit
     599        /* last bit*/
    603600        for (int iRow=numberRows_-1;iRow>rowLast;iRow--) {
    604601          aaLast[--putLast] = sparseFactor_[--get];
     
    617614        }
    618615        if (j-BLOCK<iColumn) {
    619           // save diagonal as well
     616          /* save diagonal as well*/
    620617          aPut[--put2]=diagonal_[iColumn];
    621618        }
     
    628625    block -= nBlock+ifOdd;
    629626  }
    630   factor(a,numberRows_,numberBlocks,
     627  ClpCholeskyDenseC  info;
     628  info.diagonal_=diagonal_;
     629  info.doubleParameters_[0]=doubleParameters_[10];
     630  info.integerParameters_[0]=integerParameters_[34];
     631#ifndef CLP_CILK
     632  ClpCholeskyCfactor(&info,a,numberRows_,numberBlocks,
    631633         diagonal_,workDouble_,rowsDropped);
    632   //sparseFactor_=xx;
    633   //diagonal_=yy;
    634 }
    635 // Non leaf recursive factor
    636 void
    637 ClpCholeskyDense::factor(longDouble * a, int n, int numberBlocks,
     634#else
     635  info.a=a;
     636  info.n=numberRows_;
     637  info.numberBlocks=numberBlocks;
     638  info.work=workDouble_;
     639  info.rowsDropped=rowsDropped;
     640  info.integerParameters_[1]=model_->numberThreads();
     641  ClpCholeskySpawn(&info);
     642#endif
     643  double largest=0.0;
     644  double smallest=COIN_DBL_MAX;
     645  int numberDropped=0;
     646  for (int i=0;i<numberRows_;i++) {
     647    if (diagonal_[i]) {
     648      largest = CoinMax(largest,CoinAbs(diagonal_[i]));
     649      smallest = CoinMin(smallest,CoinAbs(diagonal_[i]));
     650    } else {
     651      numberDropped++;
     652    }
     653  }
     654  doubleParameters_[3]=CoinMax(doubleParameters_[3],1.0/smallest);
     655  doubleParameters_[4]=CoinMin(doubleParameters_[4],1.0/largest);
     656  integerParameters_[20]+= numberDropped;
     657}
     658  /* Non leaf recursive factor*/
     659void 
     660ClpCholeskyCfactor(ClpCholeskyDenseC * thisStruct, longDouble * a, int n, int numberBlocks,
    638661            longDouble * diagonal, longDouble * work, int * rowsDropped)
    639662{
    640663  if (n<=BLOCK) {
    641     factorLeaf(a,n,diagonal,work,rowsDropped);
     664    ClpCholeskyCfactorLeaf(thisStruct, a,n,diagonal,work,rowsDropped);
    642665  } else {
    643666    int nb=number_blocks((n+1)>>1);
     
    647670    int nintri=(nb*(nb+1))>>1;
    648671    int nbelow=(numberBlocks-nb)*nb;
    649     factor(a,nThis,numberBlocks,diagonal,work,rowsDropped);
    650     triRec(a,nThis,a+number_entries(nb),diagonal,work,nLeft,nb,0,numberBlocks);
     672    ClpCholeskyCfactor(thisStruct, a,nThis,numberBlocks,diagonal,work,rowsDropped);
     673    ClpCholeskyCtriRec(thisStruct, a,nThis,a+number_entries(nb),diagonal,work,nLeft,nb,0,numberBlocks);
    651674    aother=a+number_entries(nintri+nbelow);
    652     recTri(a+number_entries(nb),nLeft,nThis,nb,0,aother,diagonal,work,numberBlocks);
    653     factor(aother,nLeft,
     675    ClpCholeskyCrecTri(thisStruct, a+number_entries(nb),nLeft,nThis,nb,0,aother,diagonal,work,numberBlocks);
     676    ClpCholeskyCfactor(thisStruct, aother,nLeft,
    654677        numberBlocks-nb,diagonal+nThis,work+nThis,rowsDropped);
    655678  }
    656679}
    657 // Non leaf recursive triangle rectangle update
     680  /* Non leaf recursive triangle rectangle update*/
    658681void
    659 ClpCholeskyDense::triRec(longDouble * aTri, int nThis, longDouble * aUnder,
     682ClpCholeskyCtriRec(ClpCholeskyDenseC * thisStruct, longDouble * aTri, int nThis, longDouble * aUnder,
    660683                         longDouble * diagonal, longDouble * work,
    661684                         int nLeft, int iBlock, int jBlock,
     
    663686{
    664687  if (nThis<=BLOCK&&nLeft<=BLOCK) {
    665     triRecLeaf(aTri,aUnder,diagonal,work,nLeft);
     688    ClpCholeskyCtriRecLeaf(thisStruct, aTri,aUnder,diagonal,work,nLeft);
    666689  } else if (nThis<nLeft) {
    667690    int nb=number_blocks((nLeft+1)>>1);
    668691    int nLeft2=number_rows(nb);
    669     triRec(aTri,nThis,aUnder,diagonal,work,nLeft2,iBlock,jBlock,numberBlocks);
    670     triRec(aTri,nThis,aUnder+number_entries(nb),diagonal,work,nLeft-nLeft2,
     692    ClpCholeskyCtriRec(thisStruct, aTri,nThis,aUnder,diagonal,work,nLeft2,iBlock,jBlock,numberBlocks);
     693    ClpCholeskyCtriRec(thisStruct, aTri,nThis,aUnder+number_entries(nb),diagonal,work,nLeft-nLeft2,
    671694          iBlock+nb,jBlock,numberBlocks);
    672695  } else {
     
    678701    int nintri=(nb*(nb+1))>>1;
    679702    int nbelow=(numberBlocks-nb)*nb;
    680     triRec(aTri,nThis2,aUnder,diagonal,work,nLeft,iBlock,jBlock,numberBlocks);
     703    ClpCholeskyCtriRec(thisStruct, aTri,nThis2,aUnder,diagonal,work,nLeft,iBlock,jBlock,numberBlocks);
    681704    /* and rectangular update */
    682705    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
    683706      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
    684707    aother=aUnder+number_entries(i);
    685     recRec(aTri+number_entries(nb),nThis-nThis2,nLeft,nThis2,aUnder,aother,
     708    ClpCholeskyCrecRec(thisStruct, aTri+number_entries(nb),nThis-nThis2,nLeft,nThis2,aUnder,aother,
    686709          work,kBlock,jBlock,numberBlocks);
    687     triRec(aTri+number_entries(nintri+nbelow),nThis-nThis2,aother,diagonal+nThis2,
     710    ClpCholeskyCtriRec(thisStruct, aTri+number_entries(nintri+nbelow),nThis-nThis2,aother,diagonal+nThis2,
    688711           work+nThis2,nLeft,
    689712      iBlock-nb,kBlock-nb,numberBlocks-nb);
    690713  }
    691714}
    692 // Non leaf recursive rectangle triangle update
     715  /* Non leaf recursive rectangle triangle update*/
    693716void
    694 ClpCholeskyDense::recTri(longDouble * aUnder, int nTri, int nDo,
     717ClpCholeskyCrecTri(ClpCholeskyDenseC * thisStruct, longDouble * aUnder, int nTri, int nDo,
    695718                         int iBlock, int jBlock,longDouble * aTri,
    696719                         longDouble * diagonal, longDouble * work,
     
    698721{
    699722  if (nTri<=BLOCK&&nDo<=BLOCK) {
    700     recTriLeaf(aUnder,aTri,diagonal,work,nTri);
     723    ClpCholeskyCrecTriLeaf(thisStruct, aUnder,aTri,diagonal,work,nTri);
    701724  } else if (nTri<nDo) {
    702725    int nb=number_blocks((nDo+1)>>1);
     
    704727    longDouble * aother;
    705728    int i;
    706     recTri(aUnder,nTri,nDo2,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
     729    ClpCholeskyCrecTri(thisStruct, aUnder,nTri,nDo2,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
    707730    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
    708731      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
    709732    aother=aUnder+number_entries(i);
    710     recTri(aother,nTri,nDo-nDo2,iBlock-nb,jBlock,aTri,diagonal+nDo2,
     733    ClpCholeskyCrecTri(thisStruct, aother,nTri,nDo-nDo2,iBlock-nb,jBlock,aTri,diagonal+nDo2,
    711734           work+nDo2,numberBlocks-nb);
    712735  } else {
     
    715738    longDouble * aother;
    716739    int i;
    717     recTri(aUnder,nTri2,nDo,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
     740    ClpCholeskyCrecTri(thisStruct, aUnder,nTri2,nDo,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
    718741    /* and rectangular update */
    719742    i=((numberBlocks-iBlock)*(numberBlocks-iBlock+1)-
    720743      (numberBlocks-iBlock-nb)*(numberBlocks-iBlock-nb+1))>>1;
    721744    aother=aTri+number_entries(nb);
    722     recRec(aUnder,nTri2,nTri-nTri2,nDo,aUnder+number_entries(nb),aother,
     745    ClpCholeskyCrecRec(thisStruct, aUnder,nTri2,nTri-nTri2,nDo,aUnder+number_entries(nb),aother,
    723746           work,iBlock,jBlock,numberBlocks);
    724     recTri(aUnder+number_entries(nb),nTri-nTri2,nDo,iBlock+nb,jBlock,
     747    ClpCholeskyCrecTri(thisStruct, aUnder+number_entries(nb),nTri-nTri2,nDo,iBlock+nb,jBlock,
    725748         aTri+number_entries(i),diagonal,work,numberBlocks);
    726749  }
     
    731754*/
    732755void
    733 ClpCholeskyDense::recRec(longDouble * above, int nUnder, int nUnderK,
     756ClpCholeskyCrecRec(ClpCholeskyDenseC * thisStruct, longDouble * above, int nUnder, int nUnderK,
    734757                         int nDo, longDouble * aUnder, longDouble *aOther,
    735758                         longDouble * work,
     
    739762  if (nDo<=BLOCK&&nUnder<=BLOCK&&nUnderK<=BLOCK) {
    740763    assert (nDo==BLOCK&&nUnder==BLOCK);
    741     recRecLeaf(above , aUnder ,  aOther, work, nUnderK);
     764    ClpCholeskyCrecRecLeaf(thisStruct, above , aUnder ,  aOther, work, nUnderK);
    742765  } else if (nDo<=nUnderK&&nUnder<=nUnderK) {
    743766    int nb=number_blocks((nUnderK+1)>>1);
    744767    int nUnder2=number_rows(nb);
    745     recRec(above,nUnder,nUnder2,nDo,aUnder,aOther,work,
     768    ClpCholeskyCrecRec(thisStruct, above,nUnder,nUnder2,nDo,aUnder,aOther,work,
    746769               iBlock,jBlock,numberBlocks);
    747     recRec(above,nUnder,nUnderK-nUnder2,nDo,aUnder+number_entries(nb),
     770    ClpCholeskyCrecRec(thisStruct, above,nUnder,nUnderK-nUnder2,nDo,aUnder+number_entries(nb),
    748771        aOther+number_entries(nb),work,iBlock,jBlock,numberBlocks);
    749772  } else if (nUnderK<=nDo&&nUnder<=nDo) {
     
    751774    int nDo2=number_rows(nb);
    752775    int i;
    753     recRec(above,nUnder,nUnderK,nDo2,aUnder,aOther,work,
     776    ClpCholeskyCrecRec(thisStruct, above,nUnder,nUnderK,nDo2,aUnder,aOther,work,
    754777         iBlock,jBlock,numberBlocks);
    755778    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
    756779      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
    757     recRec(above+number_entries(i),nUnder,nUnderK,nDo-nDo2,
     780    ClpCholeskyCrecRec(thisStruct, above+number_entries(i),nUnder,nUnderK,nDo-nDo2,
    758781           aUnder+number_entries(i),
    759782           aOther,work+nDo2,
     
    763786    int nUnder2=number_rows(nb);
    764787    int i;
    765     recRec(above,nUnder2,nUnderK,nDo,aUnder,aOther,work,
     788    ClpCholeskyCrecRec(thisStruct, above,nUnder2,nUnderK,nDo,aUnder,aOther,work,
    766789               iBlock,jBlock,numberBlocks);
    767790    i=((numberBlocks-iBlock)*(numberBlocks-iBlock-1)-
    768791      (numberBlocks-iBlock-nb)*(numberBlocks-iBlock-nb-1))>>1;
    769     recRec(above+number_entries(nb),nUnder-nUnder2,nUnderK,nDo,aUnder,
     792    ClpCholeskyCrecRec(thisStruct, above+number_entries(nb),nUnder-nUnder2,nUnderK,nDo,aUnder,
    770793        aOther+number_entries(i),work,iBlock+nb,jBlock,numberBlocks);
    771794  }
    772795}
    773 // Leaf recursive factor
     796  /* Leaf recursive factor*/
    774797void
    775 ClpCholeskyDense::factorLeaf(longDouble * a, int n,
     798ClpCholeskyCfactorLeaf(ClpCholeskyDenseC * thisStruct, longDouble * a, int n,
    776799                longDouble * diagonal, longDouble * work, int * rowsDropped)
    777800{
    778   CoinWorkDouble largest=doubleParameters_[3];
    779   CoinWorkDouble smallest=doubleParameters_[4];
    780   double dropValue = doubleParameters_[10];
    781   int firstPositive=integerParameters_[34];
    782   int rowOffset=diagonal-diagonal_;
    783   int numberDropped=0;
     801  double dropValue = thisStruct->doubleParameters_[0];
     802  int firstPositive=thisStruct->integerParameters_[0];
     803  int rowOffset=diagonal-thisStruct->diagonal_;
    784804  int i, j, k;
    785805  CoinWorkDouble t00,temp1;
     
    787807  aa = a-BLOCK;
    788808  for (j = 0; j < n; j ++) {
     809    bool dropColumn;
     810    CoinWorkDouble useT00;
    789811    aa+=BLOCK;
    790812    t00 = aa[j];
     
    793815      t00 -= a[j + k * BLOCK] * a[j + k * BLOCK] * multiplier;
    794816    }
    795     bool dropColumn=false;
    796     CoinWorkDouble useT00=t00;
     817    dropColumn=false;
     818    useT00=t00;
    797819    if (j+rowOffset<firstPositive) {
    798       // must be negative
     820      /* must be negative*/
    799821      if (t00<=-dropValue) {
    800         smallest = CoinMin(smallest,-t00);
    801         largest = CoinMax(largest,-t00);
    802         //aa[j]=t00;
     822        /*aa[j]=t00;*/
    803823        t00 = 1.0/t00;
    804824      } else {
    805825        dropColumn=true;
    806         //aa[j]=-1.0e100;
     826        /*aa[j]=-1.0e100;*/
    807827        useT00=-1.0e-100;
    808828        t00=0.0;
    809         integerParameters_[20]++;
    810829      }
    811830    } else {
    812       // must be positive
     831      /* must be positive*/
    813832      if (t00>=dropValue) {
    814         smallest = CoinMin(smallest,t00);
    815         largest = CoinMax(largest,t00);
    816         //aa[j]=t00;
     833        /*aa[j]=t00;*/
    817834        t00 = 1.0/t00;
    818835      } else {
    819836        dropColumn=true;
    820         //aa[j]=1.0e100;
     837        /*aa[j]=1.0e100;*/
    821838        useT00=1.0e-100;
    822839        t00=0.0;
    823         integerParameters_[20]++;
    824840      }
    825841    }
     
    837853      }
    838854    } else {
    839       // drop column
     855      /* drop column*/
    840856      rowsDropped[j+rowOffset]=2;
    841       numberDropped++;
    842857      diagonal[j]=0.0;
    843       //aa[j]=1.0e100;
     858      /*aa[j]=1.0e100;*/
    844859      work[j]=1.0e100;
    845860      for (i=j+1;i<n;i++) {
     
    848863    }
    849864  }
    850   doubleParameters_[3]=largest;
    851   doubleParameters_[4]=smallest;
    852   integerParameters_[20]+=numberDropped;
    853 }
    854 // Leaf recursive triangle rectangle update
     865}
     866  /* Leaf recursive triangle rectangle update*/
    855867void
    856 ClpCholeskyDense::triRecLeaf(longDouble * aTri, longDouble * aUnder, longDouble * diagonal, longDouble * work,
     868ClpCholeskyCtriRecLeaf(ClpCholeskyDenseC * thisStruct, longDouble * aTri, longDouble * aUnder, longDouble * diagonal, longDouble * work,
    857869                int nUnder)
    858870{
     
    862874  int irt,ict;
    863875  int it=bNumber(aTri,irt,ict);
    864   //printf("%d %d\n",iu,it);
     876  /*printf("%d %d\n",iu,it);*/
    865877  printf("trirecleaf  under (%d,%d), tri (%d,%d)\n",
    866878         iru,icu,irt,ict);
    867   assert (diagonal==diagonal_+ict*BLOCK);
     879  assert (diagonal==thisStruct->diagonal_+ict*BLOCK);
    868880#endif
    869881  int j;
     
    873885    aa = aTri-2*BLOCK;
    874886    for (j = 0; j < BLOCK; j +=2) {
    875       aa+=2*BLOCK;
     887      int i;
    876888      CoinWorkDouble temp0 = diagonal[j];
    877889      CoinWorkDouble temp1 = diagonal[j+1];
    878       for (int i=0;i<BLOCK;i+=2) {
     890      aa+=2*BLOCK;
     891      for ( i=0;i<BLOCK;i+=2) {
     892        CoinWorkDouble at1;
    879893        CoinWorkDouble t00=aUnder[i+j*BLOCK];
    880894        CoinWorkDouble t10=aUnder[i+ BLOCK + j*BLOCK];
    881895        CoinWorkDouble t01=aUnder[i+1+j*BLOCK];
    882896        CoinWorkDouble t11=aUnder[i+1+ BLOCK + j*BLOCK];
    883         for (int k = 0; k < j; ++k) {
     897        int k;
     898        for (k = 0; k < j; ++k) {
    884899          CoinWorkDouble multiplier=work[k];
    885900          CoinWorkDouble au0 = aUnder[i + k * BLOCK]*multiplier;
     
    893908        }
    894909        t00 *= temp0;
    895         CoinWorkDouble at1 = aTri[j + 1 + j*BLOCK]*work[j];
     910        at1 = aTri[j + 1 + j*BLOCK]*work[j];
    896911        t10 -= t00 * at1;
    897912        t01 *= temp0;
     
    907922    aa = aTri-BLOCK;
    908923    for (j = 0; j < BLOCK; j ++) {
     924      int i;
     925      CoinWorkDouble temp1 = diagonal[j];
    909926      aa+=BLOCK;
    910       CoinWorkDouble temp1 = diagonal[j];
    911       for (int i=0;i<nUnder;i++) {
     927      for (i=0;i<nUnder;i++) {
     928        int k;
    912929        CoinWorkDouble t00=aUnder[i+j*BLOCK];
    913         for (int k = 0; k < j; ++k) {
     930        for ( k = 0; k < j; ++k) {
    914931          CoinWorkDouble multiplier=work[k];
    915932          t00 -= aUnder[i + k * BLOCK] * aTri[j + k * BLOCK] * multiplier;
     
    922939#endif
    923940}
    924 // Leaf recursive rectangle triangle update
    925 void ClpCholeskyDense::recTriLeaf(longDouble * aUnder, longDouble * aTri,
     941  /* Leaf recursive rectangle triangle update*/
     942void ClpCholeskyCrecTriLeaf(ClpCholeskyDenseC * thisStruct, longDouble * aUnder, longDouble * aTri,
    926943                                  longDouble * diagonal, longDouble * work, int nUnder)
    927944{
     
    931948  int irt,ict;
    932949  int it=bNumber(aTri,irt,ict);
    933   //printf("%d %d\n",iu,it);
     950  /*printf("%d %d\n",iu,it);*/
    934951  printf("rectrileaf  under (%d,%d), tri (%d,%d)\n",
    935952         iru,icu,irt,ict);
    936   assert (diagonal==diagonal_+icu*BLOCK);
     953  assert (diagonal==thisStruct->diagonal_+icu*BLOCK);
    937954#endif
    938955  int i, j, k;
     
    10061023*/
    10071024void
    1008 ClpCholeskyDense::recRecLeaf(const longDouble * COIN_RESTRICT above,
     1025ClpCholeskyCrecRecLeaf(ClpCholeskyDenseC * thisStruct,
     1026                                         const longDouble * COIN_RESTRICT above,
    10091027                             const longDouble * COIN_RESTRICT aUnder,
    10101028                             longDouble * COIN_RESTRICT aOther,
     
    10191037  int iro,ico;
    10201038  int io=bNumber(aOther,iro,ico);
    1021   //printf("%d %d %d\n",ia,iu,io);
     1039  /*printf("%d %d %d\n",ia,iu,io);*/
    10221040  printf("recrecleaf above (%d,%d), under (%d,%d), other (%d,%d)\n",
    10231041         ira,ica,iru,icu,iro,ico);
     
    10281046  aa = aOther-4*BLOCK;
    10291047  if (nUnder==BLOCK) {
    1030     //#define INTEL
     1048    /*#define INTEL*/
    10311049#ifdef INTEL
    10321050    aa+=2*BLOCK;
     
    12281246  }
    12291247#endif
    1230   //longDouble * xx = sparseFactor_;
    1231   //longDouble * yy = diagonal_;
    1232   //diagonal_ = sparseFactor_ + 40000;
    1233   //sparseFactor_=diagonal_ + numberRows_;
     1248  /*longDouble * xx = sparseFactor_;*/
     1249  /*longDouble * yy = diagonal_;*/
     1250  /*diagonal_ = sparseFactor_ + 40000;*/
     1251  /*sparseFactor_=diagonal_ + numberRows_;*/
    12341252  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
    1235   // later align on boundary
     1253  /* later align on boundary*/
    12361254  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
    12371255  int iBlock;
     
    12611279    aa+=BLOCKSQ;
    12621280  }
    1263   // do diagonal outside
     1281  /* do diagonal outside*/
    12641282  for (iColumn=0;iColumn<numberRows_;iColumn++)
    12651283    region[iColumn] *= diagonal_[iColumn];
     
    12991317#endif
    13001318}
    1301 // Forward part of solve 1
     1319/* Forward part of solve 1*/
    13021320void
    13031321ClpCholeskyDense::solveF1(longDouble * a,int n,double * region)
     
    13101328      t00 -= region[k] * a[j + k * BLOCK];
    13111329    }
    1312     //t00*=a[j + j * BLOCK];
     1330    /*t00*=a[j + j * BLOCK];*/
    13131331    region[j] = t00;
    13141332  }
    13151333}
    1316 // Forward part of solve 2
     1334/* Forward part of solve 2*/
    13171335void
    13181336ClpCholeskyDense::solveF2(longDouble * a,int n,double * region, double * region2)
     
    14261444#endif
    14271445}
    1428 // Backward part of solve 1
     1446/* Backward part of solve 1*/
    14291447void
    14301448ClpCholeskyDense::solveB1(longDouble * a,int n,double * region)
     
    14371455      t00 -= region[k] * a[k + j * BLOCK];
    14381456    }
    1439     //t00*=a[j + j * BLOCK];
     1457    /*t00*=a[j + j * BLOCK];*/
    14401458    region[j] = t00;
    14411459  }
    14421460}
    1443 // Backward part of solve 2
     1461/* Backward part of solve 2*/
    14441462void
    14451463ClpCholeskyDense::solveB2(longDouble * a,int n,double * region, double * region2)
     
    15591577{
    15601578  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
    1561   // later align on boundary
     1579  /* later align on boundary*/
    15621580  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
    15631581  int iBlock;
     
    15871605    aa+=BLOCKSQ;
    15881606  }
    1589   // do diagonal outside
     1607  /* do diagonal outside*/
    15901608  for (iColumn=0;iColumn<numberRows_;iColumn++)
    15911609    region[iColumn] *= diagonal_[iColumn];
     
    16171635  }
    16181636}
    1619 // Forward part of solve 1
     1637/* Forward part of solve 1*/
    16201638void
    16211639ClpCholeskyDense::solveF1Long(longDouble * a,int n,CoinWorkDouble * region)
     
    16281646      t00 -= region[k] * a[j + k * BLOCK];
    16291647    }
    1630     //t00*=a[j + j * BLOCK];
     1648    /*t00*=a[j + j * BLOCK];*/
    16311649    region[j] = t00;
    16321650  }
    16331651}
    1634 // Forward part of solve 2
     1652/* Forward part of solve 2*/
    16351653void
    16361654ClpCholeskyDense::solveF2Long(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
     
    17441762#endif
    17451763}
    1746 // Backward part of solve 1
     1764/* Backward part of solve 1*/
    17471765void
    17481766ClpCholeskyDense::solveB1Long(longDouble * a,int n,CoinWorkDouble * region)
     
    17551773      t00 -= region[k] * a[k + j * BLOCK];
    17561774    }
    1757     //t00*=a[j + j * BLOCK];
     1775    /*t00*=a[j + j * BLOCK];*/
    17581776    region[j] = t00;
    17591777  }
    17601778}
    1761 // Backward part of solve 2
     1779/* Backward part of solve 2*/
    17621780void
    17631781ClpCholeskyDense::solveB2Long(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
     
    18771895{
    18781896  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
    1879   // later align on boundary
     1897  /* later align on boundary*/
    18801898  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
    18811899  int iBlock;
     
    19051923    aa+=BLOCKSQ;
    19061924  }
    1907   // do diagonal outside
     1925  /* do diagonal outside*/
    19081926  for (iColumn=0;iColumn<numberRows_;iColumn++)
    19091927    region[iColumn] *= diagonal_[iColumn];
     
    19351953  }
    19361954}
    1937 // Forward part of solve 1
     1955/* Forward part of solve 1*/
    19381956void
    19391957ClpCholeskyDense::solveF1LongWork(longDouble * a,int n,CoinWorkDouble * region)
     
    19461964      t00 -= region[k] * a[j + k * BLOCK];
    19471965    }
    1948     //t00*=a[j + j * BLOCK];
     1966    /*t00*=a[j + j * BLOCK];*/
    19491967    region[j] = t00;
    19501968  }
    19511969}
    1952 // Forward part of solve 2
     1970/* Forward part of solve 2*/
    19531971void
    19541972ClpCholeskyDense::solveF2LongWork(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
     
    20622080#endif
    20632081}
    2064 // Backward part of solve 1
     2082/* Backward part of solve 1*/
    20652083void
    20662084ClpCholeskyDense::solveB1LongWork(longDouble * a,int n,CoinWorkDouble * region)
     
    20732091      t00 -= region[k] * a[k + j * BLOCK];
    20742092    }
    2075     //t00*=a[j + j * BLOCK];
     2093    /*t00*=a[j + j * BLOCK];*/
    20762094    region[j] = t00;
    20772095  }
    20782096}
    2079 // Backward part of solve 2
     2097/* Backward part of solve 2*/
    20802098void
    20812099ClpCholeskyDense::solveB2LongWork(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
  • trunk/Clp/src/ClpCholeskyDense.hpp

    r1367 r1368  
    1 // Copyright (C) 2003, International Business Machines
    2 // Corporation and others.  All Rights Reserved.
     1/* Copyright (C) 2003, International Business Machines Corporation
     2   and others.  All Rights Reserved. */
    33#ifndef ClpCholeskyDense_H
    44#define ClpCholeskyDense_H
     
    77class ClpMatrixBase;
    88
    9 /** Dense class for Clp Cholesky factorization
    10 
    11 */
    129class ClpCholeskyDense : public ClpCholeskyBase {
    1310 
    1411public:
    1512   /**@name Virtual methods that the derived classes provides  */
    16    //@{
     13   /**@{*/
    1714  /** Orders rows and saves pointer to matrix.and model.
    1815   Returns non-zero if not enough memory */
     
    2825  /** Uses factorization to solve. */
    2926  virtual void solve (double * region) ;
    30   //@}
     27  /**@}*/
    3128
    32    /**@name Non virtual methods for ClpCholeskyDense  */
    33    //@{
     29  /**@name Non virtual methods for ClpCholeskyDense  */
     30  /**@{*/
    3431  /** Reserves space.
    3532      If factor not NULL then just uses passed space
     
    4239  /** part 2 of Factorize - filling in rowsDropped - blocked */
    4340  void factorizePart3(int * rowsDropped) ;
    44   /// Non leaf recursive factor
    45   void factor(longDouble * a, int n, int numberBlocks,
    46               longDouble * diagonal, longDouble * work, int * rowsDropped);
    47   /// Non leaf recursive triangle rectangle update
    48   void triRec(longDouble * aTri, int nThis, longDouble * aUnder, longDouble * diagonal, longDouble * work,
    49               int nLeft, int iBlock, int jBlock,
    50               int numberBlocks);
    51   /// Non leaf recursive rectangle triangle update
    52   void recTri(longDouble * aUnder, int nTri, int nDo,
    53               int iBlock, int jBlock,longDouble * aTri,
    54               longDouble * diagonal, longDouble * work,
    55               int numberBlocks);
     41  /** Forward part of solve */
     42  void solveF1(longDouble * a,int n,double * region);
     43  void solveF2(longDouble * a,int n,double * region,double * region2);
     44  /** Backward part of solve */
     45  void solveB1(longDouble * a,int n,double * region);
     46  void solveB2(longDouble * a,int n,double * region,double * region2);
     47  /** Uses factorization to solve. */
     48  void solveLong (CoinWorkDouble * region) ;
     49  /** Forward part of solve */
     50  void solveF1Long(longDouble * a,int n,CoinWorkDouble * region);
     51  void solveF2Long(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
     52  /** Backward part of solve */
     53  void solveB1Long(longDouble * a,int n,CoinWorkDouble * region);
     54  void solveB2Long(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
     55  /** Uses factorization to solve. */
     56  void solveLongWork (CoinWorkDouble * region) ;
     57  /** Forward part of solve */
     58  void solveF1LongWork(longDouble * a,int n,CoinWorkDouble * region);
     59  void solveF2LongWork(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
     60  /** Backward part of solve */
     61  void solveB1LongWork(longDouble * a,int n,CoinWorkDouble * region);
     62  void solveB2LongWork(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
     63  int bNumber(const longDouble * array,int &, int&);
     64  /** A */
     65  inline longDouble * aMatrix() const
     66  { return sparseFactor_;}
     67  /** Diagonal */
     68  inline longDouble * diagonal() const
     69  { return diagonal_;}
     70  /**@}*/
     71
     72
     73  /**@name Constructors, destructor */
     74  /**@{*/
     75  /** Default constructor. */
     76  ClpCholeskyDense();
     77  /** Destructor  */
     78  virtual ~ClpCholeskyDense();
     79  /** Copy */
     80  ClpCholeskyDense(const ClpCholeskyDense&);
     81  /** Assignment */
     82  ClpCholeskyDense& operator=(const ClpCholeskyDense&);
     83  /** Clone */
     84  virtual ClpCholeskyBase * clone() const ;
     85  /**@}*/
     86   
     87   
     88private:
     89  /**@name Data members */
     90  /**@{*/
     91     /** Just borrowing space */
     92  bool borrowSpace_;
     93  /**@}*/
     94};
     95
     96/* structure for C */
     97typedef struct{
     98  longDouble * diagonal_;
     99  longDouble * a;
     100  longDouble * work;
     101  int * rowsDropped;
     102  double doubleParameters_[1]; /* corresponds to 10 */
     103  int integerParameters_[2]; /* corresponds to 34, nThreads */
     104  int n;
     105  int numberBlocks;
     106} ClpCholeskyDenseC;
     107
     108extern "C" {
     109  void ClpCholeskySpawn(void *);
     110}
     111  /**Non leaf recursive factor */
     112  void
     113  ClpCholeskyCfactor(ClpCholeskyDenseC * thisStruct,
     114                     longDouble * a, int n, int numberBlocks,
     115                     longDouble * diagonal, longDouble * work, int * rowsDropped);
     116 
     117  /**Non leaf recursive triangle rectangle update */
     118  void
     119  ClpCholeskyCtriRec(ClpCholeskyDenseC * thisStruct,
     120                     longDouble * aTri, int nThis,
     121                     longDouble * aUnder, longDouble * diagonal,
     122                     longDouble * work,
     123                     int nLeft, int iBlock, int jBlock,
     124                     int numberBlocks);
     125  /**Non leaf recursive rectangle triangle update */
     126  void
     127  ClpCholeskyCrecTri(ClpCholeskyDenseC * thisStruct,
     128                     longDouble * aUnder, int nTri, int nDo,
     129                     int iBlock, int jBlock,longDouble * aTri,
     130                     longDouble * diagonal, longDouble * work,
     131                     int numberBlocks);
    56132  /** Non leaf recursive rectangle rectangle update,
    57133      nUnder is number of rows in iBlock,
    58134      nUnderK is number of rows in kBlock
    59135  */
    60   void recRec(longDouble * above, int nUnder, int nUnderK,
    61               int nDo, longDouble * aUnder, longDouble *aOther,
    62               longDouble * work,
    63               int iBlock, int jBlock,
    64               int numberBlocks);
    65   /// Leaf recursive factor
    66   void factorLeaf(longDouble * a, int n,
    67               longDouble * diagonal, longDouble * work, int * rowsDropped);
    68   /// Leaf recursive triangle rectangle update
    69   void triRecLeaf(longDouble * aTri, longDouble * aUnder,
    70                   longDouble * diagonal, longDouble * work,
    71                   int nUnder);
    72   /// Leaf recursive rectangle triangle update
    73   void recTriLeaf(longDouble * aUnder, longDouble * aTri,
    74                   longDouble * diagonal, longDouble * work, int nUnder);
     136  void
     137  ClpCholeskyCrecRec(ClpCholeskyDenseC * thisStruct,
     138                     longDouble * above, int nUnder, int nUnderK,
     139                     int nDo, longDouble * aUnder, longDouble *aOther,
     140                     longDouble * work,
     141                     int iBlock, int jBlock,
     142                     int numberBlocks);
     143  /**Leaf recursive factor */
     144  void
     145  ClpCholeskyCfactorLeaf(ClpCholeskyDenseC * thisStruct,
     146                         longDouble * a, int n,
     147                         longDouble * diagonal, longDouble * work,
     148                         int * rowsDropped);
     149  /**Leaf recursive triangle rectangle update */
     150  void
     151  ClpCholeskyCtriRecLeaf(ClpCholeskyDenseC * thisStruct,
     152                         longDouble * aTri, longDouble * aUnder,
     153                         longDouble * diagonal, longDouble * work,
     154                         int nUnder);
     155  /**Leaf recursive rectangle triangle update */
     156  void
     157  ClpCholeskyCrecTriLeaf(ClpCholeskyDenseC * thisStruct,
     158                         longDouble * aUnder, longDouble * aTri,
     159                         longDouble * diagonal, longDouble * work, int nUnder);
    75160  /** Leaf recursive rectangle rectangle update,
    76161      nUnder is number of rows in iBlock,
    77162      nUnderK is number of rows in kBlock
    78163  */
    79   void recRecLeaf(const longDouble * COIN_RESTRICT above,
    80                   const longDouble * COIN_RESTRICT aUnder,
    81                   longDouble * COIN_RESTRICT aOther,
    82                   const longDouble * COIN_RESTRICT work,
    83                   int nUnder);
    84   /// Forward part of solve
    85   void solveF1(longDouble * a,int n,double * region);
    86   void solveF2(longDouble * a,int n,double * region,double * region2);
    87   /// Backward part of solve
    88   void solveB1(longDouble * a,int n,double * region);
    89   void solveB2(longDouble * a,int n,double * region,double * region2);
    90   /** Uses factorization to solve. */
    91   void solveLong (CoinWorkDouble * region) ;
    92   /// Forward part of solve
    93   void solveF1Long(longDouble * a,int n,CoinWorkDouble * region);
    94   void solveF2Long(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
    95   /// Backward part of solve
    96   void solveB1Long(longDouble * a,int n,CoinWorkDouble * region);
    97   void solveB2Long(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
    98   /** Uses factorization to solve. */
    99   void solveLongWork (CoinWorkDouble * region) ;
    100   /// Forward part of solve
    101   void solveF1LongWork(longDouble * a,int n,CoinWorkDouble * region);
    102   void solveF2LongWork(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
    103   /// Backward part of solve
    104   void solveB1LongWork(longDouble * a,int n,CoinWorkDouble * region);
    105   void solveB2LongWork(longDouble * a,int n,CoinWorkDouble * region,CoinWorkDouble * region2);
    106   int bNumber(const longDouble * array,int &, int&);
    107   /// A
    108   inline longDouble * aMatrix() const
    109   { return sparseFactor_;}
    110   /// Diagonal
    111   inline longDouble * diagonal() const
    112   { return diagonal_;}
    113   //@}
    114 
    115 
    116   /**@name Constructors, destructor */
    117   //@{
    118   /** Default constructor. */
    119   ClpCholeskyDense();
    120   /** Destructor  */
    121   virtual ~ClpCholeskyDense();
    122   // Copy
    123   ClpCholeskyDense(const ClpCholeskyDense&);
    124   // Assignment
    125   ClpCholeskyDense& operator=(const ClpCholeskyDense&);
    126   /// Clone
    127   virtual ClpCholeskyBase * clone() const ;
    128   //@}
    129    
    130    
    131 private:
    132   /**@name Data members */
    133   //@{
    134   /// Just borrowing space
    135   bool borrowSpace_;
    136   //@}
    137 };
    138 
     164  void
     165  ClpCholeskyCrecRecLeaf(ClpCholeskyDenseC * thisStruct,
     166                         const longDouble * COIN_RESTRICT above,
     167                         const longDouble * COIN_RESTRICT aUnder,
     168                         longDouble * COIN_RESTRICT aOther,
     169                         const longDouble * COIN_RESTRICT work,
     170                         int nUnder);
    139171#endif
  • trunk/Clp/src/ClpMain.cpp

    r1366 r1368  
    190190    int choleskyType = 0;
    191191    int gamma=0;
    192     int scaleBarrier=0;
     192    parameters[whichParam(BARRIERSCALE,numberParameters,parameters)].setCurrentOption(2);
     193    int scaleBarrier=2;
    193194    int doKKT=0;
    194195    int crossover=2; // do crossover unless quadratic
     
    710711              } else if (method==ClpSolve::useBarrier||method==ClpSolve::useBarrierNoCross) {
    711712                int barrierOptions = choleskyType;
    712                 if (scaleBarrier)
    713                   barrierOptions |= 8;
     713                if (scaleBarrier) {
     714                  if ((scaleBarrier&1)!=0)
     715                    barrierOptions |= 8;
     716                  barrierOptions |= 2048*(scaleBarrier>>1);
     717                }
    714718                if (doKKT)
    715719                  barrierOptions |= 16;
  • trunk/Clp/src/ClpMessage.cpp

    r1367 r1368  
    7272  {CLP_BARRIER_ITERATION,35,1,"%d Primal %g Dual %g Complementarity %g - %d fixed, rank %d"},
    7373  {CLP_BARRIER_OBJECTIVE_GAP,36,3,"Feasible - objective gap %g"},
    74   {CLP_BARRIER_GONE_INFEASIBLE,37,2,"Gone infeasible"},
     74  {CLP_BARRIER_GONE_INFEASIBLE,37,2,"Infeasible"},
    7575  {CLP_BARRIER_CLOSE_TO_OPTIMAL,38,2,"Close to optimal after %d iterations with complementarity %g"},
    7676  {CLP_BARRIER_COMPLEMENTARITY,39,2,"Complementarity %g - %s"},
     
    8282  {CLP_BARRIER_INFO,45,3,"Detail - %s"},
    8383  {CLP_BARRIER_END,46,1,"At end primal/dual infeasibilities %g/%g, complementarity gap %g, objective %g"},
    84   {CLP_BARRIER_ACCURACY,47,2,"Relative error in phase %d, %d refinements reduced from %g to %g"},
     84  {CLP_BARRIER_ACCURACY,47,2,"Relative error in phase %d, %d passes %g => %g"},
    8585  {CLP_BARRIER_SAFE,48,2,"Initial safe primal value %g, objective norm %g"},
    8686  {CLP_BARRIER_NEGATIVE_GAPS,49,3,"%d negative gaps summing to %g"},
  • trunk/Clp/src/ClpPredictorCorrector.cpp

    r1367 r1368  
    881881    if (numberGoodTries&&handler_->logLevel()>1) {
    882882      printf("%d centering steps moved from (gap %.18g, dual %.18g, primal %.18g) to (gap %.18g, dual %.18g, primal %.18g)\n",
    883              numberGoodTries,nextGap,originalDualStep,originalPrimalStep,
    884              nextCenterGap, actualDualStep_,actualPrimalStep_);
     883             numberGoodTries,static_cast<double>(nextGap),static_cast<double>(originalDualStep),
     884             static_cast<double>(originalPrimalStep),
     885             static_cast<double>(nextCenterGap), static_cast<double>(actualDualStep_),
     886             static_cast<double>(actualPrimalStep_));
    885887    }
    886888    // save last gap
     
    991993  CoinWorkDouble maximumDualStep=COIN_DBL_MAX;
    992994  int numberTotal = numberRows_+numberColumns_;
    993 #if CLP_LONG_CHOLESKY<2
    994995  CoinWorkDouble tolerance = 1.0e-12;
    995 #else
    996   CoinWorkDouble tolerance = 1.0e-14;
    997 #endif
    998996  int chosenPrimalSequence=-1;
    999997  int chosenDualSequence=-1;
     
    15811579  CoinWorkDouble saveMaximum=0.0;
    15821580  double firstError=0.0;
    1583   double lastError=0.0;
     1581  double lastError2=0.0;
    15841582  while (!goodSolve&&numberTries<30) {
    15851583    CoinWorkDouble lastError=relativeError;
     
    16901688      firstError=relativeError;
    16911689    if (relativeError<lastError) {
    1692       lastError=relativeError;
     1690      lastError2=relativeError;
    16931691      maximumRHSChange_= maximumRHSChange;
    16941692      if (relativeError>projectionTolerance&&numberTries<=3) {
     
    17411739    }
    17421740  } /* endwhile */
    1743   if (firstError>1.0e-9||numberTries>1) {
     1741  if (firstError>1.0e-8||numberTries>1) {
    17441742    handler_->message(CLP_BARRIER_ACCURACY,messages_)
    17451743      <<phase<<numberTries<<static_cast<double>(firstError)
    1746       <<static_cast<double>(lastError)
     1744      <<static_cast<double>(lastError2)
    17471745      <<CoinMessageEol;
    17481746  }
     
    25072505      }
    25082506    }
    2509     if (0) {
    2510       int i=1324;
    2511       printf("%d %.18g %.18g %.18g %.18g %.18g %.18g %.18g\n",i,solution_[i],
    2512              diagonal_[i],dj_[i],
    2513              lowerSlack_[i],zVec_[i],
    2514              upperSlack_[i],wVec_[i]);
    2515       printf("%d %.18g %.18g %.18g %.18g %.18g\n",i,rhsC_[i],
    2516              rhsZ_[i],rhsL_[i],
    2517              rhsW_[i],rhsU_[i]);
    2518     }
    25192507#if 0
    25202508    for (int i=0;i<3;i++) {
     
    29362924{
    29372925  CoinWorkDouble complementarityMultiplier =1.0/numberComplementarityPairs_;
    2938 #if CLP_LONG_CHOLESKY<2
    29392926  const CoinWorkDouble gamma = 1.0e-8;
    29402927  const CoinWorkDouble gammap = 1.0e-8;
    29412928  CoinWorkDouble gammad = 1.0e-8;
    2942 #else
    2943   const CoinWorkDouble gamma = 1.0e-8;
    2944   const CoinWorkDouble gammap = 1.0e-8;
    2945   CoinWorkDouble gammad = 1.0e-8;
    2946 #endif
    29472929  int nextNumber;
    29482930  int nextNumberItems;
     
    31123094  }
    31133095  CoinWorkDouble offsetObjective=0.0;
    3114   const CoinWorkDouble killTolerance=primalTolerance();
     3096  CoinWorkDouble killTolerance=primalTolerance();
    31153097  CoinWorkDouble qDiagonal;
    3116 #if CLP_LONG_CHOLESKY<2
    31173098  if (mu_<1.0) {
    31183099    qDiagonal=1.0e-8;
     
    31203101    qDiagonal=1.0e-8*mu_;
    31213102  }
    3122 #else
    3123   if (mu_<1.0) {
    3124     qDiagonal=1.0e-8;
    3125   } else {
    3126     qDiagonal=1.0e-8*mu_;
    3127   }
    3128 #endif
    31293103  //CoinWorkDouble nextMu = nextGap/(static_cast<CoinWorkDouble>(2*numberComplementarityPairs_));
    31303104  //printf("using gap of %g\n",nextMu);
     
    31603134    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
    31613135#endif
     3136#ifndef CLP_CAUTION
     3137#define KILL_ITERATION 50
     3138#else
     3139#define KILL_ITERATION 100
     3140#endif
    31623141  if (!quadraticObj||1) {
    3163     if (numberIterations_<50) {
     3142    if (numberIterations_<KILL_ITERATION) {
    31643143      killFactor = 1.0;
    3165     } else if (numberIterations_<100) {
    3166       killFactor = 10.0;
     3144    } else if (numberIterations_<2*KILL_ITERATION) {
     3145      killFactor = 5.0;
    31673146      stepLength_=CoinMax(stepLength_,0.9995);
    3168     } else if (numberIterations_<200) {
    3169       killFactor = 100.0;
     3147    } else if (numberIterations_<4*KILL_ITERATION) {
     3148      killFactor = 20.0;
    31703149      stepLength_=CoinMax(stepLength_,0.99995);
    31713150    } else {
    3172       killFactor = 1.0e5;
     3151      killFactor = 1.0e2;
    31733152      stepLength_=CoinMax(stepLength_,0.999995);
    31743153    }
     
    32213200  CoinWorkDouble gamma2 = gamma_*gamma_; // gamma*gamma will be added to diagonal
    32223201  CoinWorkDouble gammaOffset=0.0;
     3202#if 0
    32233203  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    32243204  const int * columnLength = matrix_->getVectorLengths();
    32253205  const int * row = matrix_->getIndices();
    32263206  const double * element = matrix_->getElements();
     3207#endif
    32273208  for (iColumn=0;iColumn<numberTotal;iColumn++) {
    32283209    if (!flagged(iColumn)) {
     
    32403221      CoinWorkDouble sLower=extra;
    32413222      CoinWorkDouble kill;
    3242 #if CLP_LONG_CHOLESKY>1
    3243       killTolerance *= 1.0e-1;
    3244 #endif
    32453223      if (CoinAbs(newPrimal)>1.0e4) {
    32463224        kill=killTolerance*1.0e-4*newPrimal;
     
    32873265        upperSlack_[iColumn] = largeGap2;
    32883266        // so we can just have one test
    3289         fakeOldBounds=true;
     3267        fakeOldBounds=true; 
    32903268      }
    32913269      CoinWorkDouble lowerBoundInfeasibility=0.0;
    32923270      CoinWorkDouble upperBoundInfeasibility=0.0;
    3293       double saveNewPrimal = newPrimal;
     3271      //double saveNewPrimal = newPrimal;
    32943272      if (lowerBound(iColumn)) {
    32953273        CoinWorkDouble oldSlack = lowerSlack_[iColumn];
     
    33013279          newSlack = lowerSlack_[iColumn];
    33023280        CoinWorkDouble epsilon = CoinAbs(newSlack)*epsilonBase;
    3303 #if CLP_LONG_CHOLESKY<2
    3304         if (epsilon>1.0e-5) {
    3305           //cout<<"bad"<<endl;
    3306           epsilon=1.0e-5;
    3307         }
    3308 #else
    33093281        epsilon=CoinMin(epsilon,1.0e-5);
    3310 #endif
    33113282        //epsilon=1.0e-14;
    33123283        //make sure reasonable
     
    33263297            larger=feasibleSlack;
    33273298          }
    3328 #if CLP_LONG_CHOLESKY<2
    33293299          if (CoinAbs(feasibleSlack-newSlack)<1.0e-6*larger) {
    33303300            newSlack=feasibleSlack;
    33313301          }
    3332 #else
    3333           if (CoinAbs(feasibleSlack-newSlack)<1.0e-6*larger) {
    3334             newSlack=feasibleSlack;
    3335           }
    3336 #endif
    33373302        }
    33383303        if (zVec_[iColumn]>dualTolerance) {
     
    33723337          newSlack = upperSlack_[iColumn];
    33733338        CoinWorkDouble epsilon = CoinAbs(newSlack)*epsilonBase;
    3374 #if CLP_LONG_CHOLESKY<2
    3375         if (epsilon>1.0e-5) {
    3376           //cout<<"bad"<<endl;
    3377           epsilon=1.0e-5;
    3378         }
    3379 #else
    33803339        epsilon=CoinMin(epsilon,1.0e-5);
    3381 #endif
    33823340        //make sure reasonable
    33833341        //epsilon=1.0e-14;
     
    33973355            larger=feasibleSlack;
    33983356          }
    3399 #if CLP_LONG_CHOLESKY<2
    34003357          if (CoinAbs(feasibleSlack-newSlack)<1.0e-6*larger) {
    34013358            newSlack=feasibleSlack;
    34023359          }
    3403 #else
    3404           if (CoinAbs(feasibleSlack-newSlack)<1.0e-6*larger) {
    3405             newSlack=feasibleSlack;
    3406           }
    3407 #endif
    34083360        }
    34093361        if (wVec_[iColumn]>dualTolerance) {
     
    34333385        }
    34343386      }
    3435       if (false&&newPrimal!=saveNewPrimal&&iColumn<numberColumns_) {
     3387#if 0
     3388      if (newPrimal!=saveNewPrimal&&iColumn<numberColumns_) {
    34363389        // adjust slacks
    34373390        double movement = newPrimal-saveNewPrimal;
     
    34433396        }
    34443397      }
     3398#endif
    34453399      solution_[iColumn]=newPrimal;
    34463400      if (CoinAbs(newPrimal)>solutionNorm) {
     
    35813535        if (solution_[iColumn]!=lower_[iColumn]&&
    35823536            solution_[iColumn]!=upper_[iColumn]) {
    3583           printf("%d %g %g %g\n",iColumn,lower_[iColumn],
    3584                  solution_[iColumn],upper_[iColumn]);
     3537          printf("%d %g %g %g\n",iColumn,static_cast<double>(lower_[iColumn]),
     3538                 static_cast<double>(solution_[iColumn]),static_cast<double>(upper_[iColumn]));
    35853539        }
    35863540        diagonal_[iColumn]=0.0;
  • trunk/Clp/src/ClpSimplexPrimal.cpp

    r1367 r1368  
    816816    dummy=4;
    817817    matrix_->generalExpanded(this,9,dummy);
    818 #ifdef KEEP_GOING_IF_FIXED
     818#ifdef CLP_CAUTION
    819819    double lastAverageInfeasibility=sumDualInfeasibilities_/
    820820      static_cast<double>(numberDualInfeasibilities_+10);
     
    823823    double sumInfeasibility =  nonLinearCost_->sumInfeasibilities();
    824824    int reason2=0;
    825 #ifdef KEEP_GOING_IF_FIXED
     825#ifdef CLP_CAUTION
     826#if CLP_CAUTION==2
     827    double test2=1.0e5;
     828#else
     829    double test2=1.0e-1;
     830#endif
    826831    if (!lastSumInfeasibility&&sumInfeasibility&&
    827          lastAverageInfeasibility<1.0e-1&&numberPivots>10)
     832         lastAverageInfeasibility<test2&&numberPivots>10)
    828833      reason2=3;
    829834    if (lastSumInfeasibility<1.0e-6&&sumInfeasibility>1.0e-3&&
  • trunk/Clp/src/ClpSolve.cpp

    r1367 r1368  
    3636#include "ClpCholeskyWssmp.hpp"
    3737#include "ClpCholeskyWssmpKKT.hpp"
    38 #define FAST_BARRIER
    3938#endif
    4039#ifdef UFL_BARRIER
    4140#include "ClpCholeskyUfl.hpp"
    42 #define FAST_BARRIER
    4341#endif
    4442#ifdef TAUCS_BARRIER
    4543#include "ClpCholeskyTaucs.hpp"
    46 #define FAST_BARRIER
    4744#endif
    4845#ifdef MUMPS_BARRIER
    4946#include "ClpCholeskyMumps.hpp"
    50 #define FAST_BARRIER
    51 #endif
    52 #ifdef COIN_DEVELOP
    53 #ifndef FAST_BARRIER
    54 static int numberBarrier=0;
    55 #endif
    5647#endif
    5748#ifdef COIN_HAS_VOL
     
    18761867    bool doKKT=false;
    18771868    bool forceFixing=false;
     1869    int speed=0;
    18781870    if (barrierOptions&16) {
    18791871      barrierOptions &= ~16;
     
    18961888      barrier.setProjectionTolerance(1.0e-9);
    18971889    }
     1890    if (barrierOptions&(2048|4096)) {
     1891      speed = (barrierOptions&(2048|4096))>>11;
     1892      barrierOptions &= ~(2048|4096);
     1893    }
    18981894    if (barrierOptions&8) {
    18991895      barrierOptions &= ~8;
    19001896      scale=true;
    19011897    }
    1902 #ifdef COIN_DEVELOP
    1903 #ifndef FAST_BARRIER
    1904     if (!numberBarrier)
    1905       std::cout<<"Warning - the default ordering is just on row counts! "
    1906                <<"The factorization is being improved"<<std::endl;
    1907     numberBarrier++;
    1908 #endif
    1909 #endif
    19101898    // If quadratic force KKT
    19111899    if (quadraticObj) {
     
    19171905      if (!doKKT) {
    19181906        ClpCholeskyBase * cholesky = new ClpCholeskyBase();
     1907        cholesky->setIntegerParameter(0,speed);
    19191908        barrier.setCholesky(cholesky);
    19201909      } else {
     
    19881977    if (saveMaxIts<1000) {
    19891978      barrier.setMaximumBarrierIterations(saveMaxIts);
    1990       model2->setMaximumIterations(1000000);
     1979      model2->setMaximumIterations(10000000);
    19911980    }
    19921981#ifndef SAVEIT
     
    23942383      }
    23952384    }
    2396     model2->setMaximumIterations(saveMaxIts);
     2385   
     2386    //model2->setMaximumIterations(saveMaxIts);
    23972387#ifdef BORROW
    23982388    delete [] rowPrimal;
     
    25182508    delete model2;
    25192509  }
    2520   setMaximumIterations(saveMaxIterations);
     2510  if (method!=ClpSolve::useBarrierNoCross&&
     2511      method!=ClpSolve::useBarrier)
     2512    setMaximumIterations(saveMaxIterations);
    25212513  std::string statusMessage[]={"Unknown","Optimal","PrimalInfeasible","DualInfeasible","Stopped",
    25222514                               "Errors","User stopped"};
Note: See TracChangeset for help on using the changeset viewer.