Changeset 2385 for trunk/Clp/src/ClpCholeskyPardiso.cpp
 Timestamp:
 Jan 6, 2019 2:43:06 PM (4 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/ClpCholeskyPardiso.cpp
r2236 r2385 21 21 // Default Constructor 22 22 // 23 ClpCholeskyPardiso::ClpCholeskyPardiso 24 25 { 26 23 ClpCholeskyPardiso::ClpCholeskyPardiso(int denseThreshold) 24 : ClpCholeskyBase(denseThreshold) 25 { 26 type_ = 17; 27 27 } 28 28 … … 30 30 // Copy constructor 31 31 // 32 ClpCholeskyPardiso::ClpCholeskyPardiso (const ClpCholeskyPardiso & rhs) 33 : ClpCholeskyBase(rhs) 34 { 35 } 36 32 ClpCholeskyPardiso::ClpCholeskyPardiso(const ClpCholeskyPardiso &rhs) 33 : ClpCholeskyBase(rhs) 34 { 35 } 37 36 38 37 // 39 38 // Destructor 40 39 // 41 ClpCholeskyPardiso::~ClpCholeskyPardiso 40 ClpCholeskyPardiso::~ClpCholeskyPardiso() 42 41 { 43 42 /*  */ 44 43 /* .. Termination and release of memory. */ 45 44 /*  */ 46 int phase = 1; 47 void * voidParameters=reinterpret_cast<void *>(doubleParameters_);48 MKL_INT mtype = 2; 49 MKL_INT nrhs = 1; 50 memset(integerParameters_, 0,sizeof(integerParameters_));51 MKL_INT maxfct, mnum, error, msglvl =0;45 int phase = 1; /* Release internal memory. */ 46 void *voidParameters = reinterpret_cast< void * >(doubleParameters_); 47 MKL_INT mtype = 2; /* Real symmetric matrix */ 48 MKL_INT nrhs = 1; /* Number of right hand sides. */ 49 memset(integerParameters_, 0, sizeof(integerParameters_)); 50 MKL_INT maxfct, mnum, error, msglvl = 0; 52 51 /* Auxiliary variables. */ 53 double ddum; 54 MKL_INT idum; 55 PARDISO 56 57 58 52 double ddum; /* Double dummy */ 53 MKL_INT idum; /* Integer dummy. */ 54 PARDISO(voidParameters, &maxfct, &mnum, &mtype, &phase, 55 &numberRows_, sparseFactor_, 56 choleskyStart_, choleskyRow_, &idum, 57 &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error); 59 58 } 60 59 … … 63 62 // 64 63 ClpCholeskyPardiso & 65 ClpCholeskyPardiso::operator=(const ClpCholeskyPardiso &rhs)66 { 67 68 69 70 64 ClpCholeskyPardiso::operator=(const ClpCholeskyPardiso &rhs) 65 { 66 if (this != &rhs) { 67 ClpCholeskyBase::operator=(rhs); 68 } 69 return *this; 71 70 } 72 71 // 73 72 // Clone 74 73 // 75 ClpCholeskyBase * 76 { 77 74 ClpCholeskyBase *ClpCholeskyPardiso::clone() const 75 { 76 return new ClpCholeskyPardiso(*this); 78 77 } 79 78 /* Orders rows and saves pointer to matrix.and model */ 80 int 81 ClpCholeskyPardiso::order(ClpInterior * model) 82 { 83 numberRows_ = model>numberRows(); 84 rowsDropped_ = new char [numberRows_]; 85 memset(rowsDropped_, 0, numberRows_); 86 numberRowsDropped_ = 0; 87 model_ = model; 88 rowCopy_ = model>clpMatrix()>reverseOrderedCopy(); 89 // Space for starts 90 choleskyStart_ = new CoinBigIndex[numberRows_+1]; 91 const CoinBigIndex * columnStart = model_>clpMatrix()>getVectorStarts(); 92 const int * columnLength = model_>clpMatrix()>getVectorLengths(); 93 const int * row = model_>clpMatrix()>getIndices(); 94 const CoinBigIndex * rowStart = rowCopy_>getVectorStarts(); 95 const int * rowLength = rowCopy_>getVectorLengths(); 96 const int * column = rowCopy_>getIndices(); 97 // We need two arrays for counts 98 int * which = new int [numberRows_]; 99 int * used = new int[numberRows_+1]; 100 CoinZeroN(used, numberRows_); 101 int iRow; 102 sizeFactor_ = 0; 103 int numberColumns = model>numberColumns(); 104 int numberDense = 0; 105 denseThreshold_=0; 106 if (denseThreshold_ > 0) { 107 delete [] whichDense_; 108 delete [] denseColumn_; 109 delete dense_; 110 whichDense_ = new char[numberColumns]; 111 int iColumn; 112 used[numberRows_] = 0; 113 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 114 int length = columnLength[iColumn]; 115 used[length] += 1; 79 int ClpCholeskyPardiso::order(ClpInterior *model) 80 { 81 numberRows_ = model>numberRows(); 82 rowsDropped_ = new char[numberRows_]; 83 memset(rowsDropped_, 0, numberRows_); 84 numberRowsDropped_ = 0; 85 model_ = model; 86 rowCopy_ = model>clpMatrix()>reverseOrderedCopy(); 87 // Space for starts 88 choleskyStart_ = new CoinBigIndex[numberRows_ + 1]; 89 const CoinBigIndex *columnStart = model_>clpMatrix()>getVectorStarts(); 90 const int *columnLength = model_>clpMatrix()>getVectorLengths(); 91 const int *row = model_>clpMatrix()>getIndices(); 92 const CoinBigIndex *rowStart = rowCopy_>getVectorStarts(); 93 const int *rowLength = rowCopy_>getVectorLengths(); 94 const int *column = rowCopy_>getIndices(); 95 // We need two arrays for counts 96 int *which = new int[numberRows_]; 97 int *used = new int[numberRows_ + 1]; 98 CoinZeroN(used, numberRows_); 99 int iRow; 100 sizeFactor_ = 0; 101 int numberColumns = model>numberColumns(); 102 int numberDense = 0; 103 denseThreshold_ = 0; 104 if (denseThreshold_ > 0) { 105 delete[] whichDense_; 106 delete[] denseColumn_; 107 delete dense_; 108 whichDense_ = new char[numberColumns]; 109 int iColumn; 110 used[numberRows_] = 0; 111 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 112 int length = columnLength[iColumn]; 113 used[length] += 1; 114 } 115 int nLong = 0; 116 int stop = CoinMax(denseThreshold_ / 2, 100); 117 for (iRow = numberRows_; iRow >= stop; iRow) { 118 if (used[iRow]) 119 COIN_DETAIL_PRINT(printf("%d columns are of length %d\n", used[iRow], iRow)); 120 nLong += used[iRow]; 121 if (nLong > 50  nLong > (numberColumns >> 2)) 122 break; 123 } 124 CoinZeroN(used, numberRows_); 125 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 126 if (columnLength[iColumn] < denseThreshold_) { 127 whichDense_[iColumn] = 0; 128 } else { 129 whichDense_[iColumn] = 1; 130 numberDense++; 131 } 132 } 133 if (!numberDense  numberDense > 100) { 134 // free 135 delete[] whichDense_; 136 whichDense_ = NULL; 137 denseColumn_ = NULL; 138 dense_ = NULL; 139 } else { 140 // space for dense columns 141 denseColumn_ = new double[numberDense * numberRows_]; 142 // dense cholesky 143 dense_ = new ClpCholeskyDense(); 144 dense_>reserveSpace(NULL, numberDense); 145 COIN_DETAIL_PRINT(printf("Taking %d columns as dense\n", numberDense)); 146 } 147 } 148 for (iRow = 0; iRow < numberRows_; iRow++) { 149 int number = 1; 150 // make sure diagonal exists 151 which[0] = iRow; 152 used[iRow] = 1; 153 if (!rowsDropped_[iRow]) { 154 CoinBigIndex startRow = rowStart[iRow]; 155 CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; 156 for (CoinBigIndex k = startRow; k < endRow; k++) { 157 int iColumn = column[k]; 158 if (!whichDense_  !whichDense_[iColumn]) { 159 CoinBigIndex start = columnStart[iColumn]; 160 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 161 for (CoinBigIndex j = start; j < end; j++) { 162 int jRow = row[j]; 163 if (jRow >= iRow && !rowsDropped_[jRow]) { 164 if (!used[jRow]) { 165 used[jRow] = 1; 166 which[number++] = jRow; 167 } 168 } 116 169 } 117 int nLong = 0; 118 int stop = CoinMax(denseThreshold_ / 2, 100); 119 for (iRow = numberRows_; iRow >= stop; iRow) { 120 if (used[iRow]) 121 COIN_DETAIL_PRINT(printf("%d columns are of length %d\n", used[iRow], iRow)); 122 nLong += used[iRow]; 123 if (nLong > 50  nLong > (numberColumns >> 2)) 124 break; 170 } 171 } 172 sizeFactor_ += number; 173 int j; 174 for (j = 0; j < number; j++) 175 used[which[j]] = 0; 176 } 177 } 178 delete[] which; 179 // Now we have size  create arrays and fill in 180 try { 181 choleskyRow_ = new int[sizeFactor_]; 182 } catch (...) { 183 // no memory 184 delete[] choleskyStart_; 185 choleskyStart_ = NULL; 186 return 1; 187 } 188 try { 189 sparseFactor_ = new double[sizeFactor_]; 190 } catch (...) { 191 // no memory 192 delete[] choleskyRow_; 193 choleskyRow_ = NULL; 194 delete[] choleskyStart_; 195 choleskyStart_ = NULL; 196 return 1; 197 } 198 199 sizeFactor_ = 0; 200 which = choleskyRow_; 201 for (iRow = 0; iRow < numberRows_; iRow++) { 202 int number = 1; 203 // make sure diagonal exists 204 which[0] = iRow; 205 used[iRow] = 1; 206 choleskyStart_[iRow] = sizeFactor_; 207 if (!rowsDropped_[iRow]) { 208 CoinBigIndex startRow = rowStart[iRow]; 209 CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; 210 for (CoinBigIndex k = startRow; k < endRow; k++) { 211 int iColumn = column[k]; 212 if (!whichDense_  !whichDense_[iColumn]) { 213 CoinBigIndex start = columnStart[iColumn]; 214 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 215 for (CoinBigIndex j = start; j < end; j++) { 216 int jRow = row[j]; 217 if (jRow >= iRow && !rowsDropped_[jRow]) { 218 if (!used[jRow]) { 219 used[jRow] = 1; 220 which[number++] = jRow; 221 } 222 } 125 223 } 126 CoinZeroN(used, numberRows_); 127 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 128 if (columnLength[iColumn] < denseThreshold_) { 129 whichDense_[iColumn] = 0; 130 } else { 131 whichDense_[iColumn] = 1; 132 numberDense++; 133 } 134 } 135 if (!numberDense  numberDense > 100) { 136 // free 137 delete [] whichDense_; 138 whichDense_ = NULL; 139 denseColumn_ = NULL; 140 dense_ = NULL; 141 } else { 142 // space for dense columns 143 denseColumn_ = new double [numberDense*numberRows_]; 144 // dense cholesky 145 dense_ = new ClpCholeskyDense(); 146 dense_>reserveSpace(NULL, numberDense); 147 COIN_DETAIL_PRINT(printf("Taking %d columns as dense\n", numberDense)); 148 } 149 } 150 for (iRow = 0; iRow < numberRows_; iRow++) { 151 int number = 1; 152 // make sure diagonal exists 153 which[0] = iRow; 154 used[iRow] = 1; 155 if (!rowsDropped_[iRow]) { 156 CoinBigIndex startRow = rowStart[iRow]; 157 CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; 158 for (CoinBigIndex k = startRow; k < endRow; k++) { 159 int iColumn = column[k]; 160 if (!whichDense_  !whichDense_[iColumn]) { 161 CoinBigIndex start = columnStart[iColumn]; 162 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 163 for (CoinBigIndex j = start; j < end; j++) { 164 int jRow = row[j]; 165 if (jRow >= iRow && !rowsDropped_[jRow]) { 166 if (!used[jRow]) { 167 used[jRow] = 1; 168 which[number++] = jRow; 169 } 170 } 171 } 172 } 173 } 174 sizeFactor_ += number; 175 int j; 176 for (j = 0; j < number; j++) 177 used[which[j]] = 0; 178 } 179 } 180 delete [] which; 181 // Now we have size  create arrays and fill in 182 try { 183 choleskyRow_ = new int [sizeFactor_]; 184 } catch (...) { 185 // no memory 186 delete [] choleskyStart_; 187 choleskyStart_ = NULL; 188 return 1; 189 } 190 try { 191 sparseFactor_ = new double[sizeFactor_]; 192 } catch (...) { 193 // no memory 194 delete [] choleskyRow_; 195 choleskyRow_ = NULL; 196 delete [] choleskyStart_; 197 choleskyStart_ = NULL; 198 return 1; 199 } 200 201 sizeFactor_ = 0; 202 which = choleskyRow_; 203 for (iRow = 0; iRow < numberRows_; iRow++) { 204 int number = 1; 205 // make sure diagonal exists 206 which[0] = iRow; 207 used[iRow] = 1; 208 choleskyStart_[iRow] = sizeFactor_; 209 if (!rowsDropped_[iRow]) { 210 CoinBigIndex startRow = rowStart[iRow]; 211 CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; 212 for (CoinBigIndex k = startRow; k < endRow; k++) { 213 int iColumn = column[k]; 214 if (!whichDense_  !whichDense_[iColumn]) { 215 CoinBigIndex start = columnStart[iColumn]; 216 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 217 for (CoinBigIndex j = start; j < end; j++) { 218 int jRow = row[j]; 219 if (jRow >= iRow && !rowsDropped_[jRow]) { 220 if (!used[jRow]) { 221 used[jRow] = 1; 222 which[number++] = jRow; 223 } 224 } 225 } 226 } 227 } 228 sizeFactor_ += number; 229 int j; 230 for (j = 0; j < number; j++) 231 used[which[j]] = 0; 232 // Sort 233 std::sort(which, which + number); 234 // move which on 235 which += number; 236 } 237 } 238 choleskyStart_[numberRows_] = sizeFactor_; 239 delete [] used; 240 permuteInverse_ = new int [numberRows_]; 241 permute_ = new int[numberRows_]; 242 // Assume void * same size as double 243 void * voidParameters=reinterpret_cast<void *>(doubleParameters_); 244 MKL_INT mtype = 2; /* Real symmetric matrix */ 245 MKL_INT nrhs = 1; /* Number of right hand sides. */ 246 memset(integerParameters_,0,sizeof(integerParameters_)); 247 MKL_INT maxfct, mnum, phase, error, msglvl; 248 /* Auxiliary variables. */ 249 double ddum; /* Double dummy */ 250 MKL_INT idum; /* Integer dummy. */ 251 /*  */ 252 /* .. Setup Pardiso control parameters. */ 253 /*  */ 254 integerParameters_[0] = 1; /* No solver default */ 255 integerParameters_[1] = 2; /* Fillin reordering from METIS */ 256 integerParameters_[3] = 0; /* No iterativedirect algorithm */ 257 integerParameters_[4] = 0; /* No user fillin reducing permutation */ 258 integerParameters_[5] = 0; /* Write solution into x */ 259 integerParameters_[6] = 0; /* Number of refinements done */ 260 integerParameters_[7] = 2; /* Max numbers of iterative refinement steps */ 261 integerParameters_[8] = 0; /* Not in use */ 262 integerParameters_[9] = 13; /* Perturb the pivot elements with 1E13 */ 263 integerParameters_[10] = 1; /* Use nonsymmetric permutation and scaling MPS */ 264 integerParameters_[11] = 0; /* Not in use */ 265 integerParameters_[12] = 1; /* Maximum weighted matching algorithm is switchedoff (default for symmetric). Try integerParameters_[12] = 1 in case of inappropriate accuracy */ 266 integerParameters_[13] = 0; /* Output: Number of perturbed pivots */ 267 integerParameters_[14] = 0; /* Not in use */ 268 integerParameters_[15] = 0; /* Not in use */ 269 integerParameters_[16] = 0; /* Not in use */ 270 integerParameters_[17] = 1; /* Output: Number of nonzeros in the factor LU */ 271 integerParameters_[18] = 1; /* Output: Mflops for LU factorization */ 272 integerParameters_[19] = 0; /* Output: Numbers of CG Iterations */ 273 integerParameters_[19] = 1; // for 2x2 pivoting 274 integerParameters_[34] = 1; /* PARDISO use Cstyle indexing for ia and ja arrays */ 275 integerParameters_[55] = 1; /* Pivoting control */ 276 maxfct = 1; /* Maximum number of numerical factorizations. */ 277 mnum = 1; /* Which factorization to use. */ 278 msglvl = 0; /* Print statistical information in file */ 279 error = 0; /* Initialize error flag */ 280 /*  */ 281 /* .. Initialize the internal solver memory pointer. This is only */ 282 /* necessary for the FIRST call of the PARDISO solver. */ 283 /*  */ 284 memset(voidParameters,0,sizeof(doubleParameters_)); 285 /*  */ 286 /* .. Reordering and Symbolic Factorization. This step also allocates */ 287 /* all memory that is necessary for the factorization. */ 288 /*  */ 289 // temp  fill elements 290 for (int i=0;i<sizeFactor_;i++) 291 sparseFactor_[i]=1.0; 224 } 225 } 226 sizeFactor_ += number; 227 int j; 228 for (j = 0; j < number; j++) 229 used[which[j]] = 0; 230 // Sort 231 std::sort(which, which + number); 232 // move which on 233 which += number; 234 } 235 } 236 choleskyStart_[numberRows_] = sizeFactor_; 237 delete[] used; 238 permuteInverse_ = new int[numberRows_]; 239 permute_ = new int[numberRows_]; 240 // Assume void * same size as double 241 void *voidParameters = reinterpret_cast< void * >(doubleParameters_); 242 MKL_INT mtype = 2; /* Real symmetric matrix */ 243 MKL_INT nrhs = 1; /* Number of right hand sides. */ 244 memset(integerParameters_, 0, sizeof(integerParameters_)); 245 MKL_INT maxfct, mnum, phase, error, msglvl; 246 /* Auxiliary variables. */ 247 double ddum; /* Double dummy */ 248 MKL_INT idum; /* Integer dummy. */ 249 /*  */ 250 /* .. Setup Pardiso control parameters. */ 251 /*  */ 252 integerParameters_[0] = 1; /* No solver default */ 253 integerParameters_[1] = 2; /* Fillin reordering from METIS */ 254 integerParameters_[3] = 0; /* No iterativedirect algorithm */ 255 integerParameters_[4] = 0; /* No user fillin reducing permutation */ 256 integerParameters_[5] = 0; /* Write solution into x */ 257 integerParameters_[6] = 0; /* Number of refinements done */ 258 integerParameters_[7] = 2; /* Max numbers of iterative refinement steps */ 259 integerParameters_[8] = 0; /* Not in use */ 260 integerParameters_[9] = 13; /* Perturb the pivot elements with 1E13 */ 261 integerParameters_[10] = 1; /* Use nonsymmetric permutation and scaling MPS */ 262 integerParameters_[11] = 0; /* Not in use */ 263 integerParameters_[12] = 1; /* Maximum weighted matching algorithm is switchedoff (default for symmetric). Try integerParameters_[12] = 1 in case of inappropriate accuracy */ 264 integerParameters_[13] = 0; /* Output: Number of perturbed pivots */ 265 integerParameters_[14] = 0; /* Not in use */ 266 integerParameters_[15] = 0; /* Not in use */ 267 integerParameters_[16] = 0; /* Not in use */ 268 integerParameters_[17] = 1; /* Output: Number of nonzeros in the factor LU */ 269 integerParameters_[18] = 1; /* Output: Mflops for LU factorization */ 270 integerParameters_[19] = 0; /* Output: Numbers of CG Iterations */ 271 integerParameters_[19] = 1; // for 2x2 pivoting 272 integerParameters_[34] = 1; /* PARDISO use Cstyle indexing for ia and ja arrays */ 273 integerParameters_[55] = 1; /* Pivoting control */ 274 maxfct = 1; /* Maximum number of numerical factorizations. */ 275 mnum = 1; /* Which factorization to use. */ 276 msglvl = 0; /* Print statistical information in file */ 277 error = 0; /* Initialize error flag */ 278 /*  */ 279 /* .. Initialize the internal solver memory pointer. This is only */ 280 /* necessary for the FIRST call of the PARDISO solver. */ 281 /*  */ 282 memset(voidParameters, 0, sizeof(doubleParameters_)); 283 /*  */ 284 /* .. Reordering and Symbolic Factorization. This step also allocates */ 285 /* all memory that is necessary for the factorization. */ 286 /*  */ 287 // temp  fill elements 288 for (int i = 0; i < sizeFactor_; i++) 289 sparseFactor_[i] = 1.0; 292 290 #define ALWAYS_ORDER 293 291 #ifdef ALWAYS_ORDER 294 error=0;292 error = 0; 295 293 #else 296 297 PARDISO(voidParameters, &maxfct, &mnum, &mtype, &phase,298 299 300 294 phase = 11; 295 PARDISO(voidParameters, &maxfct, &mnum, &mtype, &phase, 296 &numberRows_, sparseFactor_, 297 choleskyStart_, choleskyRow_, &idum, 298 &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error); 301 299 #endif 302 if ( error != 0 ) 303 { 304 printf ("ERROR during symbolic factorization: %d\n", error); 305 return 1; 306 } 307 printf ("Number of nonzeros in factors = %d\n", integerParameters_[17]); 308 printf ("Number of factorization MFLOPS = %d\n", integerParameters_[18]); 309 lastNumberDropped_=0; 310 // Modify gamma and delta if 0.0 311 if (!model>gamma() && !model>delta()) { 312 model>setGamma(5.0e5); 313 model>setDelta(5.0e5); 314 } 315 return 0; 300 if (error != 0) { 301 printf("ERROR during symbolic factorization: %d\n", error); 302 return 1; 303 } 304 printf("Number of nonzeros in factors = %d\n", integerParameters_[17]); 305 printf("Number of factorization MFLOPS = %d\n", integerParameters_[18]); 306 lastNumberDropped_ = 0; 307 // Modify gamma and delta if 0.0 308 if (!model>gamma() && !model>delta()) { 309 model>setGamma(5.0e5); 310 model>setDelta(5.0e5); 311 } 312 return 0; 316 313 } 317 314 /* Does Symbolic factorization given permutation. … … 319 316 user must provide factorize and solve. Otherwise the default factorization is used 320 317 returns nonzero if not enough memory */ 321 int 322 ClpCholeskyPardiso::symbolic() 323 { 324 return 0; 318 int ClpCholeskyPardiso::symbolic() 319 { 320 return 0; 325 321 } 326 322 /* Factorize  filling in rowsDropped and returning number dropped */ 327 int 328 ClpCholeskyPardiso::factorize(const double * diagonal, int * rowsDropped) 329 { 330 const CoinBigIndex * columnStart = model_>clpMatrix()>getVectorStarts(); 331 const int * columnLength = model_>clpMatrix()>getVectorLengths(); 332 const int * row = model_>clpMatrix()>getIndices(); 333 const double * element = model_>clpMatrix()>getElements(); 334 const CoinBigIndex * rowStart = rowCopy_>getVectorStarts(); 335 const int * rowLength = rowCopy_>getVectorLengths(); 336 const int * column = rowCopy_>getIndices(); 337 const double * elementByRow = rowCopy_>getElements(); 338 int numberColumns = model_>clpMatrix()>getNumCols(); 339 int iRow; 340 double * work = new double[numberRows_]; 341 CoinZeroN(work, numberRows_); 342 const double * diagonalSlack = diagonal + numberColumns; 343 double largest; 344 //int numberDense = 0; 345 //if (dense_) 346 // numberDense = dense_>numberRows(); 347 //perturbation 348 double perturbation = model_>diagonalPerturbation() * model_>diagonalNorm(); 349 perturbation = perturbation * perturbation; 350 if (perturbation > 1.0) { 323 int ClpCholeskyPardiso::factorize(const double *diagonal, int *rowsDropped) 324 { 325 const CoinBigIndex *columnStart = model_>clpMatrix()>getVectorStarts(); 326 const int *columnLength = model_>clpMatrix()>getVectorLengths(); 327 const int *row = model_>clpMatrix()>getIndices(); 328 const double *element = model_>clpMatrix()>getElements(); 329 const CoinBigIndex *rowStart = rowCopy_>getVectorStarts(); 330 const int *rowLength = rowCopy_>getVectorLengths(); 331 const int *column = rowCopy_>getIndices(); 332 const double *elementByRow = rowCopy_>getElements(); 333 int numberColumns = model_>clpMatrix()>getNumCols(); 334 int iRow; 335 double *work = new double[numberRows_]; 336 CoinZeroN(work, numberRows_); 337 const double *diagonalSlack = diagonal + numberColumns; 338 double largest; 339 //int numberDense = 0; 340 //if (dense_) 341 // numberDense = dense_>numberRows(); 342 //perturbation 343 double perturbation = model_>diagonalPerturbation() * model_>diagonalNorm(); 344 perturbation = perturbation * perturbation; 345 if (perturbation > 1.0) { 351 346 #ifdef COIN_DEVELOP 352 353 347 //if (model_>model()>logLevel()&4) 348 std::cout << "large perturbation " << perturbation << std::endl; 354 349 #endif 355 perturbation = sqrt(perturbation);; 356 perturbation = 1.0; 357 } 358 if (whichDense_) { 359 double * denseDiagonal = dense_>diagonal(); 360 double * dense = denseColumn_; 361 int iDense = 0; 362 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 363 if (whichDense_[iColumn]) { 364 denseDiagonal[iDense++] = 1.0 / diagonal[iColumn]; 365 CoinZeroN(dense, numberRows_); 366 CoinBigIndex start = columnStart[iColumn]; 367 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 368 for (CoinBigIndex j = start; j < end; j++) { 369 int jRow = row[j]; 370 dense[jRow] = element[j]; 371 } 372 dense += numberRows_; 373 } 350 perturbation = sqrt(perturbation); 351 ; 352 perturbation = 1.0; 353 } 354 if (whichDense_) { 355 double *denseDiagonal = dense_>diagonal(); 356 double *dense = denseColumn_; 357 int iDense = 0; 358 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 359 if (whichDense_[iColumn]) { 360 denseDiagonal[iDense++] = 1.0 / diagonal[iColumn]; 361 CoinZeroN(dense, numberRows_); 362 CoinBigIndex start = columnStart[iColumn]; 363 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 364 for (CoinBigIndex j = start; j < end; j++) { 365 int jRow = row[j]; 366 dense[jRow] = element[j]; 367 } 368 dense += numberRows_; 369 } 370 } 371 } 372 double delta2 = model_>delta(); // add delta*delta to diagonal 373 delta2 *= delta2; 374 int *whichNon = new int[numberRows_]; 375 int nNon = 0; 376 for (iRow = 0; iRow < numberRows_; iRow++) { 377 //for (int i=0;i<numberRows_;i++) 378 //assert(!work[i]); 379 for (int i = 0; i < nNon; i++) 380 work[whichNon[i]] = 0.0; 381 nNon = 1; 382 whichNon[0] = iRow; 383 double *put = sparseFactor_ + choleskyStart_[iRow]; 384 int *which = choleskyRow_ + choleskyStart_[iRow]; 385 int number = choleskyStart_[iRow + 1]  choleskyStart_[iRow]; 386 if (!rowLength[iRow]) 387 rowsDropped_[iRow] = 1; 388 if (!rowsDropped_[iRow]) { 389 CoinBigIndex startRow = rowStart[iRow]; 390 CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; 391 work[iRow] = diagonalSlack[iRow] + delta2; 392 for (CoinBigIndex k = startRow; k < endRow; k++) { 393 int iColumn = column[k]; 394 if (!whichDense_  !whichDense_[iColumn]) { 395 CoinBigIndex start = columnStart[iColumn]; 396 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 397 double multiplier = diagonal[iColumn] * elementByRow[k]; 398 for (CoinBigIndex j = start; j < end; j++) { 399 int jRow = row[j]; 400 if (jRow >= iRow && !rowsDropped_[jRow]) { 401 double value = element[j] * multiplier; 402 if (!work[jRow]) 403 whichNon[nNon++] = jRow; 404 work[jRow] += value; 405 } 374 406 } 375 } 376 double delta2 = model_>delta(); // add delta*delta to diagonal 377 delta2 *= delta2; 378 int * whichNon = new int [numberRows_]; 379 int nNon=0; 380 for (iRow = 0; iRow < numberRows_; iRow++) { 381 //for (int i=0;i<numberRows_;i++) 382 //assert(!work[i]); 383 for (int i=0;i<nNon;i++) 384 work[whichNon[i]]=0.0; 385 nNon=1; 386 whichNon[0]=iRow; 387 double * put = sparseFactor_ + choleskyStart_[iRow]; 388 int * which = choleskyRow_ + choleskyStart_[iRow]; 389 int number = choleskyStart_[iRow+1]  choleskyStart_[iRow]; 390 if (!rowLength[iRow]) 391 rowsDropped_[iRow] = 1; 392 if (!rowsDropped_[iRow]) { 393 CoinBigIndex startRow = rowStart[iRow]; 394 CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; 395 work[iRow] = diagonalSlack[iRow] + delta2; 396 for (CoinBigIndex k = startRow; k < endRow; k++) { 397 int iColumn = column[k]; 398 if (!whichDense_  !whichDense_[iColumn]) { 399 CoinBigIndex start = columnStart[iColumn]; 400 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 401 double multiplier = diagonal[iColumn] * elementByRow[k]; 402 for (CoinBigIndex j = start; j < end; j++) { 403 int jRow = row[j]; 404 if (jRow >= iRow && !rowsDropped_[jRow]) { 405 double value = element[j] * multiplier; 406 if (!work[jRow]) 407 whichNon[nNon++]=jRow; 408 work[jRow] += value; 409 } 410 } 411 } 412 } 413 int j; 414 for (j = 0; j < number; j++) { 415 int jRow = which[j]; 416 put[j] = work[jRow]; 417 work[jRow] = 0.0; 418 } 419 } else { 420 // dropped 421 int j; 422 for (j = 1; j < number; j++) { 423 put[j] = 0.0; 424 } 425 put[0] = 1.0e12; 407 } 408 } 409 int j; 410 for (j = 0; j < number; j++) { 411 int jRow = which[j]; 412 put[j] = work[jRow]; 413 work[jRow] = 0.0; 414 } 415 } else { 416 // dropped 417 int j; 418 for (j = 1; j < number; j++) { 419 put[j] = 0.0; 420 } 421 put[0] = 1.0e12; 422 } 423 } 424 delete[] whichNon; 425 //check sizes 426 double largest2 = maximumAbsElement(sparseFactor_, sizeFactor_); 427 largest2 *= 1.0e20; 428 largest = CoinMin(largest2, 1.0e11); 429 int numberDroppedBefore = 0; 430 for (iRow = 0; iRow < numberRows_; iRow++) { 431 int dropped = rowsDropped_[iRow]; 432 // Move to int array 433 rowsDropped[iRow] = dropped; 434 if (!dropped) { 435 CoinBigIndex start = choleskyStart_[iRow]; 436 double diagonal = sparseFactor_[start]; 437 if (diagonal > largest2) { 438 sparseFactor_[start] = diagonal + perturbation; 439 } else { 440 sparseFactor_[start] = diagonal + perturbation; 441 rowsDropped[iRow] = 2; 442 numberDroppedBefore++; 443 } 444 } 445 } 446 //integerParameters_[0] = 0; 447 int maxfct = 1; /* Maximum number of numerical factorizations. */ 448 int mnum = 1; /* Which factorization to use. */ 449 int msglvl = 0; /* Print statistical information in file */ 450 int error = 0; /* Initialize error flag */ 451 int phase = 22; 452 #ifdef ALWAYS_ORDER 453 { 454 int phase = 1; /* Release internal memory. */ 455 void *voidParameters = reinterpret_cast< void * >(doubleParameters_); 456 MKL_INT mtype = 2; /* Real symmetric matrix */ 457 MKL_INT nrhs = 1; /* Number of right hand sides. */ 458 memset(integerParameters_, 0, sizeof(integerParameters_)); 459 MKL_INT maxfct, mnum, error, msglvl = 0; 460 /* Auxiliary variables. */ 461 double ddum; /* Double dummy */ 462 MKL_INT idum; /* Integer dummy. */ 463 PARDISO(voidParameters, &maxfct, &mnum, &mtype, &phase, 464 &numberRows_, sparseFactor_, 465 choleskyStart_, choleskyRow_, &idum, 466 &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error); 467 memset(integerParameters_, 0, sizeof(integerParameters_)); 468 /*  */ 469 /* .. Setup Pardiso control parameters. */ 470 /*  */ 471 integerParameters_[0] = 1; /* No solver default */ 472 integerParameters_[1] = 2; /* Fillin reordering from METIS */ 473 integerParameters_[3] = 0; /* No iterativedirect algorithm */ 474 integerParameters_[4] = 0; /* No user fillin reducing permutation */ 475 integerParameters_[5] = 0; /* Write solution into x */ 476 integerParameters_[6] = 0; /* Number of refinements done */ 477 integerParameters_[7] = 2; /* Max numbers of iterative refinement steps */ 478 integerParameters_[8] = 0; /* Not in use */ 479 integerParameters_[9] = 13; /* Perturb the pivot elements with 1E13 */ 480 integerParameters_[10] = 1; /* Use nonsymmetric permutation and scaling MPS */ 481 integerParameters_[11] = 0; /* Not in use */ 482 integerParameters_[12] = 1; /* Maximum weighted matching algorithm is switchedoff (default for symmetric). Try integerParameters_[12] = 1 in case of inappropriate accuracy */ 483 integerParameters_[13] = 0; /* Output: Number of perturbed pivots */ 484 integerParameters_[14] = 0; /* Not in use */ 485 integerParameters_[15] = 0; /* Not in use */ 486 integerParameters_[16] = 0; /* Not in use */ 487 integerParameters_[17] = 1; /* Output: Number of nonzeros in the factor LU */ 488 integerParameters_[18] = 1; /* Output: Mflops for LU factorization */ 489 integerParameters_[19] = 0; /* Output: Numbers of CG Iterations */ 490 integerParameters_[19] = 1; // for 2x2 pivoting 491 integerParameters_[34] = 1; /* PARDISO use Cstyle indexing for ia and ja arrays */ 492 integerParameters_[55] = 1; /* Pivoting control */ 493 maxfct = 1; /* Maximum number of numerical factorizations. */ 494 mnum = 1; /* Which factorization to use. */ 495 msglvl = 0; /* Print statistical information in file */ 496 error = 0; /* Initialize error flag */ 497 /*  */ 498 /* .. Initialize the internal solver memory pointer. This is only */ 499 /* necessary for the FIRST call of the PARDISO solver. */ 500 /*  */ 501 memset(voidParameters, 0, sizeof(doubleParameters_)); 502 /*  */ 503 /* .. Reordering and Symbolic Factorization. This step also allocates */ 504 /* all memory that is necessary for the factorization. */ 505 /*  */ 506 // pack down 507 CoinBigIndex numberElements = 0; 508 for (int iRow = 0; iRow < numberRows_; iRow++) { 509 CoinBigIndex start = choleskyStart_[iRow]; 510 CoinBigIndex end = choleskyStart_[iRow + 1]; 511 choleskyStart_[iRow] = numberElements; 512 for (CoinBigIndex j = start; j < end; j++) { 513 if (sparseFactor_[j]) { 514 choleskyRow_[numberElements] = choleskyRow_[j]; 515 sparseFactor_[numberElements++] = sparseFactor_[j]; 516 } 517 } 518 } 519 choleskyStart_[numberRows_] = numberElements; 520 phase = 11; 521 PARDISO(voidParameters, &maxfct, &mnum, &mtype, &phase, 522 &numberRows_, sparseFactor_, 523 choleskyStart_, choleskyRow_, &idum, 524 &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error); 525 if (error != 0) { 526 printf("ERROR during symbolic factorization: %d\n", error); 527 return 1; 528 } 529 printf("Number of nonzeros in factors = %d\n", integerParameters_[17]); 530 printf("Number of factorization MFLOPS = %d\n", integerParameters_[18]); 531 } 532 #endif 533 if (numberRowsDropped_ > lastNumberDropped_ && false) { 534 phase = 12; //reorder 535 CoinBigIndex numberElements = 0; 536 for (int iRow = 0; iRow < numberRows_; iRow++) { 537 CoinBigIndex start = choleskyStart_[iRow]; 538 CoinBigIndex end = choleskyStart_[iRow + 1]; 539 choleskyStart_[iRow] = numberElements; 540 if (rowsDropped_[iRow] != 2) { 541 for (CoinBigIndex j = start; j < end; j++) { 542 choleskyRow_[numberElements] = choleskyRow_[j]; 543 sparseFactor_[numberElements++] = sparseFactor_[j]; 544 } 545 } else { 546 rowsDropped_[iRow] = 1; 547 for (CoinBigIndex j = start; j < end; j++) { 548 if (sparseFactor_[j]) { 549 choleskyRow_[numberElements] = choleskyRow_[j]; 550 sparseFactor_[numberElements++] = sparseFactor_[j]; 426 551 } 427 } 428 delete [] whichNon; 429 //check sizes 430 double largest2 = maximumAbsElement(sparseFactor_, sizeFactor_); 431 largest2 *= 1.0e20; 432 largest = CoinMin(largest2, 1.0e11); 433 int numberDroppedBefore = 0; 434 for (iRow = 0; iRow < numberRows_; iRow++) { 435 int dropped = rowsDropped_[iRow]; 436 // Move to int array 437 rowsDropped[iRow] = dropped; 438 if (!dropped) { 439 CoinBigIndex start = choleskyStart_[iRow]; 440 double diagonal = sparseFactor_[start]; 441 if (diagonal > largest2) { 442 sparseFactor_[start] = diagonal + perturbation; 443 } else { 444 sparseFactor_[start] = diagonal + perturbation; 445 rowsDropped[iRow] = 2; 446 numberDroppedBefore++; 447 } 448 } 449 } 450 //integerParameters_[0] = 0; 451 int maxfct = 1; /* Maximum number of numerical factorizations. */ 452 int mnum = 1; /* Which factorization to use. */ 453 int msglvl = 0; /* Print statistical information in file */ 454 int error = 0; /* Initialize error flag */ 455 int phase = 22; 456 #ifdef ALWAYS_ORDER 457 { 458 int phase = 1; /* Release internal memory. */ 459 void * voidParameters=reinterpret_cast<void *>(doubleParameters_); 460 MKL_INT mtype = 2; /* Real symmetric matrix */ 461 MKL_INT nrhs = 1; /* Number of right hand sides. */ 462 memset(integerParameters_,0,sizeof(integerParameters_)); 463 MKL_INT maxfct, mnum, error, msglvl=0; 464 /* Auxiliary variables. */ 465 double ddum; /* Double dummy */ 466 MKL_INT idum; /* Integer dummy. */ 467 PARDISO (voidParameters, &maxfct, &mnum, &mtype, &phase, 468 &numberRows_, sparseFactor_, 469 choleskyStart_, choleskyRow_, &idum, 470 &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error); 471 memset(integerParameters_,0,sizeof(integerParameters_)); 472 /*  */ 473 /* .. Setup Pardiso control parameters. */ 474 /*  */ 475 integerParameters_[0] = 1; /* No solver default */ 476 integerParameters_[1] = 2; /* Fillin reordering from METIS */ 477 integerParameters_[3] = 0; /* No iterativedirect algorithm */ 478 integerParameters_[4] = 0; /* No user fillin reducing permutation */ 479 integerParameters_[5] = 0; /* Write solution into x */ 480 integerParameters_[6] = 0; /* Number of refinements done */ 481 integerParameters_[7] = 2; /* Max numbers of iterative refinement steps */ 482 integerParameters_[8] = 0; /* Not in use */ 483 integerParameters_[9] = 13; /* Perturb the pivot elements with 1E13 */ 484 integerParameters_[10] = 1; /* Use nonsymmetric permutation and scaling MPS */ 485 integerParameters_[11] = 0; /* Not in use */ 486 integerParameters_[12] = 1; /* Maximum weighted matching algorithm is switchedoff (default for symmetric). Try integerParameters_[12] = 1 in case of inappropriate accuracy */ 487 integerParameters_[13] = 0; /* Output: Number of perturbed pivots */ 488 integerParameters_[14] = 0; /* Not in use */ 489 integerParameters_[15] = 0; /* Not in use */ 490 integerParameters_[16] = 0; /* Not in use */ 491 integerParameters_[17] = 1; /* Output: Number of nonzeros in the factor LU */ 492 integerParameters_[18] = 1; /* Output: Mflops for LU factorization */ 493 integerParameters_[19] = 0; /* Output: Numbers of CG Iterations */ 494 integerParameters_[19] = 1; // for 2x2 pivoting 495 integerParameters_[34] = 1; /* PARDISO use Cstyle indexing for ia and ja arrays */ 496 integerParameters_[55] = 1; /* Pivoting control */ 497 maxfct = 1; /* Maximum number of numerical factorizations. */ 498 mnum = 1; /* Which factorization to use. */ 499 msglvl = 0; /* Print statistical information in file */ 500 error = 0; /* Initialize error flag */ 501 /*  */ 502 /* .. Initialize the internal solver memory pointer. This is only */ 503 /* necessary for the FIRST call of the PARDISO solver. */ 504 /*  */ 505 memset(voidParameters,0,sizeof(doubleParameters_)); 506 /*  */ 507 /* .. Reordering and Symbolic Factorization. This step also allocates */ 508 /* all memory that is necessary for the factorization. */ 509 /*  */ 510 // pack down 511 CoinBigIndex numberElements=0; 512 for (int iRow = 0; iRow < numberRows_; iRow++) { 513 CoinBigIndex start = choleskyStart_[iRow]; 514 CoinBigIndex end = choleskyStart_[iRow+1]; 515 choleskyStart_[iRow]=numberElements; 516 for (CoinBigIndex j=start;j<end;j++) { 517 if (sparseFactor_[j]) { 518 choleskyRow_[numberElements]=choleskyRow_[j]; 519 sparseFactor_[numberElements++]=sparseFactor_[j]; 520 } 521 } 522 } 523 choleskyStart_[numberRows_]=numberElements; 524 phase = 11; 525 PARDISO (voidParameters, &maxfct, &mnum, &mtype, &phase, 526 &numberRows_, sparseFactor_, 527 choleskyStart_, choleskyRow_, &idum, 528 &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error); 529 if ( error != 0 ) 530 { 531 printf ("ERROR during symbolic factorization: %d\n", error); 532 return 1; 533 } 534 printf ("Number of nonzeros in factors = %d\n", integerParameters_[17]); 535 printf ("Number of factorization MFLOPS = %d\n", integerParameters_[18]); 536 } 537 #endif 538 if (numberRowsDropped_>lastNumberDropped_&&false) { 539 phase=12; //reorder 540 CoinBigIndex numberElements=0; 541 for (int iRow = 0; iRow < numberRows_; iRow++) { 542 CoinBigIndex start = choleskyStart_[iRow]; 543 CoinBigIndex end = choleskyStart_[iRow+1]; 544 choleskyStart_[iRow]=numberElements; 545 if (rowsDropped_[iRow]!=2) { 546 for (CoinBigIndex j=start;j<end;j++) { 547 choleskyRow_[numberElements]=choleskyRow_[j]; 548 sparseFactor_[numberElements++]=sparseFactor_[j]; 549 } 550 } else { 551 rowsDropped_[iRow]=1; 552 for (CoinBigIndex j=start;j<end;j++) { 553 if (sparseFactor_[j]) { 554 choleskyRow_[numberElements]=choleskyRow_[j]; 555 sparseFactor_[numberElements++]=sparseFactor_[j]; 556 } 557 } 558 } 559 } 560 choleskyStart_[numberRows_]=numberElements; 561 lastNumberDropped_=numberRowsDropped_; 562 } 563 MKL_INT mtype = 2; /* Real symmetric matrix */ 564 void * voidParameters=reinterpret_cast<void *>(doubleParameters_); 565 MKL_INT nrhs = 1; /* Number of right hand sides. */ 566 /* Auxiliary variables. */ 567 double ddum; /* Double dummy */ 568 MKL_INT idum; /* Integer dummy. */ 569 PARDISO (voidParameters, &maxfct, &mnum, &mtype, &phase, 570 &numberRows_, sparseFactor_, 571 choleskyStart_, choleskyRow_, &idum, 572 &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error); 573 printf("%d positive, %d negative\n",integerParameters_[21], 574 integerParameters_[22]); 575 // use some other array? 576 double * originalDiagonal=new double[numberRows_]; 577 pardiso_getdiag(voidParameters,work,originalDiagonal,&mnum,&error); 578 delete [] originalDiagonal; 579 largest=1.0e11; 580 double smallest=COIN_DBL_MAX; 581 int newDropped = 0; 582 for (int i=0;i<numberRows_;i++) { 583 double value=work[i]; 584 if (value!=1.0e100) { 585 //if (value>1.0e13) 586 //printf("large diagonal %g at %d\n",value,i); 587 largest=CoinMax(largest,value); 588 smallest=CoinMin(smallest,value); 589 } else if (!rowsDropped_[i]) { 590 rowsDropped_[i]=2; 591 rowsDropped[i]=2; 592 printf("%d dropped\n",i); 593 newDropped++; 594 numberRowsDropped_++; 595 } 596 } 597 if(smallest<=0.0) 598 printf("small negative diagonal %g\n",smallest); 599 delete [] work; 600 if (model_>messageHandler()>logLevel() > 1) 601 std::cout << "Cholesky  largest " << largest << " smallest " << smallest << std::endl; 602 choleskyCondition_ = largest / smallest; 603 status_ = 0; 604 return newDropped; 552 } 553 } 554 } 555 choleskyStart_[numberRows_] = numberElements; 556 lastNumberDropped_ = numberRowsDropped_; 557 } 558 MKL_INT mtype = 2; /* Real symmetric matrix */ 559 void *voidParameters = reinterpret_cast< void * >(doubleParameters_); 560 MKL_INT nrhs = 1; /* Number of right hand sides. */ 561 /* Auxiliary variables. */ 562 double ddum; /* Double dummy */ 563 MKL_INT idum; /* Integer dummy. */ 564 PARDISO(voidParameters, &maxfct, &mnum, &mtype, &phase, 565 &numberRows_, sparseFactor_, 566 choleskyStart_, choleskyRow_, &idum, 567 &nrhs, integerParameters_, &msglvl, &ddum, &ddum, &error); 568 printf("%d positive, %d negative\n", integerParameters_[21], 569 integerParameters_[22]); 570 // use some other array? 571 double *originalDiagonal = new double[numberRows_]; 572 pardiso_getdiag(voidParameters, work, originalDiagonal, &mnum, &error); 573 delete[] originalDiagonal; 574 largest = 1.0e11; 575 double smallest = COIN_DBL_MAX; 576 int newDropped = 0; 577 for (int i = 0; i < numberRows_; i++) { 578 double value = work[i]; 579 if (value != 1.0e100) { 580 //if (value>1.0e13) 581 //printf("large diagonal %g at %d\n",value,i); 582 largest = CoinMax(largest, value); 583 smallest = CoinMin(smallest, value); 584 } else if (!rowsDropped_[i]) { 585 rowsDropped_[i] = 2; 586 rowsDropped[i] = 2; 587 printf("%d dropped\n", i); 588 newDropped++; 589 numberRowsDropped_++; 590 } 591 } 592 if (smallest <= 0.0) 593 printf("small negative diagonal %g\n", smallest); 594 delete[] work; 595 if (model_>messageHandler()>logLevel() > 1) 596 std::cout << "Cholesky  largest " << largest << " smallest " << smallest << std::endl; 597 choleskyCondition_ = largest / smallest; 598 status_ = 0; 599 return newDropped; 605 600 } 606 601 /* Uses factorization to solve. */ 607 void 608 ClpCholeskyPardiso::solve (double * region) 602 void ClpCholeskyPardiso::solve(double *region) 609 603 { 610 604 /*  */ … … 612 606 /*  */ 613 607 int phase = 33; 614 integerParameters_[7] = 2; 615 double * 616 memcpy(x, region,numberRows_*sizeof(double));617 MKL_INT mtype = 2; 618 void * voidParameters=reinterpret_cast<void *>(doubleParameters_);619 MKL_INT nrhs = 1; 620 MKL_INT maxfct =1, mnum=1, error, msglvl=0;608 integerParameters_[7] = 2; /* Max numbers of iterative refinement steps. */ 609 double *x = new double[numberRows_]; 610 memcpy(x, region, numberRows_ * sizeof(double)); 611 MKL_INT mtype = 2; /* Real symmetric matrix */ 612 void *voidParameters = reinterpret_cast< void * >(doubleParameters_); 613 MKL_INT nrhs = 1; /* Number of right hand sides. */ 614 MKL_INT maxfct = 1, mnum = 1, error, msglvl = 0; 621 615 /* Auxiliary variables. */ 622 MKL_INT idum; /* Integer dummy. */ 623 PARDISO (voidParameters, &maxfct, &mnum, &mtype, &phase, 624 &numberRows_, sparseFactor_, 625 choleskyStart_, choleskyRow_, &idum, 626 &nrhs, integerParameters_, &msglvl, x,region, &error); 627 if ( error != 0 ) 628 { 629 printf ("ERROR during solution: %d\n", error); 630 exit (3); 631 } 632 delete [] x; 616 MKL_INT idum; /* Integer dummy. */ 617 PARDISO(voidParameters, &maxfct, &mnum, &mtype, &phase, 618 &numberRows_, sparseFactor_, 619 choleskyStart_, choleskyRow_, &idum, 620 &nrhs, integerParameters_, &msglvl, x, region, &error); 621 if (error != 0) { 622 printf("ERROR during solution: %d\n", error); 623 exit(3); 624 } 625 delete[] x; 633 626 } 634 627 /*  */ … … 636 629 /* .. in pardiso solver */ 637 630 /*  */ 638 631 // up rows dropped 639 632 extern "C" { 640 int mkl_pardiso_pivot( const double*aii, double*bii, const double*eps )633 int mkl_pardiso_pivot(const double *aii, double *bii, const double *eps) 641 634 { 642 635 //if ( (*bii > *eps)  ( *bii < *eps ) ) 643 if (*bii >*eps)644 645 if ( *bii > 0)646 647 648 649 650 636 if (*bii > *eps) 637 return 0; 638 if (*bii > 0) 639 *bii = *eps; 640 else 641 *bii = *eps; 642 *bii = 1.0e100; 643 return 1; 651 644 } 652 645 } 653 646 #endif 647 648 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 649 */
Note: See TracChangeset
for help on using the changeset viewer.