289 | | void |
290 | | ClpPackedMatrix::times(double scalar, |
291 | | const double * COIN_RESTRICT x, double * COIN_RESTRICT y) const |
292 | | { |
293 | | int iRow, iColumn; |
294 | | // get matrix data pointers |
295 | | const int * COIN_RESTRICT row = matrix_->getIndices(); |
296 | | const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
297 | | const double * COIN_RESTRICT elementByColumn = matrix_->getElements(); |
298 | | //memset(y,0,matrix_->getNumRows()*sizeof(double)); |
299 | | assert (((flags_&0x02) != 0) == matrix_->hasGaps()); |
300 | | if (!(flags_ & 2)) { |
301 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
302 | | CoinBigIndex j; |
303 | | double value = x[iColumn]; |
304 | | if (value) { |
305 | | CoinBigIndex start = columnStart[iColumn]; |
306 | | CoinBigIndex end = columnStart[iColumn+1]; |
307 | | value *= scalar; |
308 | | for (j = start; j < end; j++) { |
309 | | iRow = row[j]; |
310 | | y[iRow] += value * elementByColumn[j]; |
311 | | } |
312 | | } |
313 | | } |
314 | | } else { |
315 | | const int * columnLength = matrix_->getVectorLengths(); |
316 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
317 | | CoinBigIndex j; |
318 | | double value = x[iColumn]; |
319 | | if (value) { |
320 | | CoinBigIndex start = columnStart[iColumn]; |
321 | | CoinBigIndex end = start + columnLength[iColumn]; |
322 | | value *= scalar; |
323 | | for (j = start; j < end; j++) { |
324 | | iRow = row[j]; |
325 | | y[iRow] += value * elementByColumn[j]; |
326 | | } |
327 | | } |
328 | | } |
329 | | } |
| 294 | void ClpPackedMatrix::times(double scalar, |
| 295 | const double *COIN_RESTRICT x, double *COIN_RESTRICT y) const |
| 296 | { |
| 297 | int iRow, iColumn; |
| 298 | // get matrix data pointers |
| 299 | const int *COIN_RESTRICT row = matrix_->getIndices(); |
| 300 | const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
| 301 | const double *COIN_RESTRICT elementByColumn = matrix_->getElements(); |
| 302 | //memset(y,0,matrix_->getNumRows()*sizeof(double)); |
| 303 | assert(((flags_ & 0x02) != 0) == matrix_->hasGaps()); |
| 304 | if (!(flags_ & 2)) { |
| 305 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 306 | CoinBigIndex j; |
| 307 | double value = x[iColumn]; |
| 308 | if (value) { |
| 309 | CoinBigIndex start = columnStart[iColumn]; |
| 310 | CoinBigIndex end = columnStart[iColumn + 1]; |
| 311 | value *= scalar; |
| 312 | for (j = start; j < end; j++) { |
| 313 | iRow = row[j]; |
| 314 | y[iRow] += value * elementByColumn[j]; |
| 315 | } |
| 316 | } |
| 317 | } |
| 318 | } else { |
| 319 | const int *columnLength = matrix_->getVectorLengths(); |
| 320 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 321 | CoinBigIndex j; |
| 322 | double value = x[iColumn]; |
| 323 | if (value) { |
| 324 | CoinBigIndex start = columnStart[iColumn]; |
| 325 | CoinBigIndex end = start + columnLength[iColumn]; |
| 326 | value *= scalar; |
| 327 | for (j = start; j < end; j++) { |
| 328 | iRow = row[j]; |
| 329 | y[iRow] += value * elementByColumn[j]; |
| 330 | } |
| 331 | } |
| 332 | } |
| 333 | } |
368 | | int numberThreads=abcState(); |
369 | | if (numberThreads) { |
370 | | clpTempInfo info[ABOCA_LITE]; |
371 | | int chunk = (numberActiveColumns_+numberThreads-1)/numberThreads; |
372 | | int n=0; |
373 | | for (int i=0;i<numberThreads;i++) { |
374 | | info[i].spare=y; |
375 | | info[i].work=const_cast<double *>(x); |
376 | | info[i].startColumn=n; |
377 | | info[i].element=elementByColumn; |
378 | | info[i].start=columnStart; |
379 | | info[i].row=row; |
380 | | info[i].numberToDo=CoinMin(chunk,numberActiveColumns_-n); |
381 | | n += chunk; |
382 | | } |
383 | | for (int i=0;i<numberThreads;i++) { |
384 | | cilk_spawn transposeTimesBit(info[i]); |
385 | | } |
386 | | cilk_sync; |
387 | | } else { |
388 | | #endif |
389 | | CoinBigIndex start = columnStart[0]; |
390 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
391 | | CoinBigIndex j; |
392 | | CoinBigIndex next = columnStart[iColumn+1]; |
393 | | double value = y[iColumn]; |
394 | | for (j = start; j < next; j++) { |
395 | | int jRow = row[j]; |
396 | | value -= x[jRow] * elementByColumn[j]; |
397 | | } |
398 | | start = next; |
399 | | y[iColumn] = value; |
400 | | } |
| 371 | int numberThreads = abcState(); |
| 372 | if (numberThreads) { |
| 373 | clpTempInfo info[ABOCA_LITE]; |
| 374 | int chunk = (numberActiveColumns_ + numberThreads - 1) / numberThreads; |
| 375 | int n = 0; |
| 376 | for (int i = 0; i < numberThreads; i++) { |
| 377 | info[i].spare = y; |
| 378 | info[i].work = const_cast< double * >(x); |
| 379 | info[i].startColumn = n; |
| 380 | info[i].element = elementByColumn; |
| 381 | info[i].start = columnStart; |
| 382 | info[i].row = row; |
| 383 | info[i].numberToDo = CoinMin(chunk, numberActiveColumns_ - n); |
| 384 | n += chunk; |
| 385 | } |
| 386 | for (int i = 0; i < numberThreads; i++) { |
| 387 | cilk_spawn transposeTimesBit(info[i]); |
| 388 | } |
| 389 | cilk_sync; |
| 390 | } else { |
| 391 | #endif |
| 392 | CoinBigIndex start = columnStart[0]; |
| 393 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 394 | CoinBigIndex j; |
| 395 | CoinBigIndex next = columnStart[iColumn + 1]; |
| 396 | double value = y[iColumn]; |
| 397 | for (j = start; j < next; j++) { |
| 398 | int jRow = row[j]; |
| 399 | value -= x[jRow] * elementByColumn[j]; |
| 400 | } |
| 401 | start = next; |
| 402 | y[iColumn] = value; |
| 403 | } |
402 | | } |
403 | | #endif |
404 | | } else { |
405 | | CoinBigIndex start = columnStart[0]; |
406 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
407 | | CoinBigIndex j; |
408 | | CoinBigIndex next = columnStart[iColumn+1]; |
409 | | double value = 0.0; |
410 | | for (j = start; j < next; j++) { |
411 | | int jRow = row[j]; |
412 | | value += x[jRow] * elementByColumn[j]; |
413 | | } |
414 | | start = next; |
415 | | y[iColumn] += value * scalar; |
416 | | } |
417 | | } |
418 | | } else { |
419 | | const int * columnLength = matrix_->getVectorLengths(); |
| 405 | } |
| 406 | #endif |
| 407 | } else { |
| 408 | CoinBigIndex start = columnStart[0]; |
| 409 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 410 | CoinBigIndex j; |
| 411 | CoinBigIndex next = columnStart[iColumn + 1]; |
| 412 | double value = 0.0; |
| 413 | for (j = start; j < next; j++) { |
| 414 | int jRow = row[j]; |
| 415 | value += x[jRow] * elementByColumn[j]; |
| 416 | } |
| 417 | start = next; |
| 418 | y[iColumn] += value * scalar; |
| 419 | } |
| 420 | } |
| 421 | } else { |
| 422 | const int *columnLength = matrix_->getVectorLengths(); |
| 423 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 424 | CoinBigIndex j; |
| 425 | double value = 0.0; |
| 426 | CoinBigIndex start = columnStart[iColumn]; |
| 427 | CoinBigIndex end = start + columnLength[iColumn]; |
| 428 | for (j = start; j < end; j++) { |
| 429 | int jRow = row[j]; |
| 430 | value += x[jRow] * elementByColumn[j]; |
| 431 | } |
| 432 | y[iColumn] += value * scalar; |
| 433 | } |
| 434 | } |
| 435 | } |
| 436 | void ClpPackedMatrix::times(double scalar, |
| 437 | const double *COIN_RESTRICT x, double *COIN_RESTRICT y, |
| 438 | const double *COIN_RESTRICT rowScale, |
| 439 | const double *COIN_RESTRICT columnScale) const |
| 440 | { |
| 441 | if (rowScale) { |
| 442 | int iRow, iColumn; |
| 443 | // get matrix data pointers |
| 444 | const int *COIN_RESTRICT row = matrix_->getIndices(); |
| 445 | const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
| 446 | const double *COIN_RESTRICT elementByColumn = matrix_->getElements(); |
| 447 | if (!(flags_ & 2)) { |
| 448 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 449 | double value = x[iColumn]; |
| 450 | if (value) { |
| 451 | // scaled |
| 452 | value *= scalar * columnScale[iColumn]; |
| 453 | CoinBigIndex start = columnStart[iColumn]; |
| 454 | CoinBigIndex end = columnStart[iColumn + 1]; |
| 455 | CoinBigIndex j; |
| 456 | for (j = start; j < end; j++) { |
| 457 | iRow = row[j]; |
| 458 | y[iRow] += value * elementByColumn[j] * rowScale[iRow]; |
| 459 | } |
| 460 | } |
| 461 | } |
| 462 | } else { |
| 463 | const int *COIN_RESTRICT columnLength = matrix_->getVectorLengths(); |
| 464 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 465 | double value = x[iColumn]; |
| 466 | if (value) { |
| 467 | // scaled |
| 468 | value *= scalar * columnScale[iColumn]; |
| 469 | CoinBigIndex start = columnStart[iColumn]; |
| 470 | CoinBigIndex end = start + columnLength[iColumn]; |
| 471 | CoinBigIndex j; |
| 472 | for (j = start; j < end; j++) { |
| 473 | iRow = row[j]; |
| 474 | y[iRow] += value * elementByColumn[j] * rowScale[iRow]; |
| 475 | } |
| 476 | } |
| 477 | } |
| 478 | } |
| 479 | } else { |
| 480 | times(scalar, x, y); |
| 481 | } |
| 482 | } |
| 483 | void ClpPackedMatrix::transposeTimes(double scalar, |
| 484 | const double *COIN_RESTRICT x, double *COIN_RESTRICT y, |
| 485 | const double *COIN_RESTRICT rowScale, |
| 486 | const double *COIN_RESTRICT columnScale, |
| 487 | double *COIN_RESTRICT spare) const |
| 488 | { |
| 489 | if (rowScale) { |
| 490 | int iColumn; |
| 491 | // get matrix data pointers |
| 492 | const int *COIN_RESTRICT row = matrix_->getIndices(); |
| 493 | const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
| 494 | const int *COIN_RESTRICT columnLength = matrix_->getVectorLengths(); |
| 495 | const double *COIN_RESTRICT elementByColumn = matrix_->getElements(); |
| 496 | if (!spare) { |
| 497 | if (!(flags_ & 2)) { |
| 498 | CoinBigIndex start = columnStart[0]; |
| 499 | if (scalar == -1.0) { |
421 | | CoinBigIndex j; |
422 | | double value = 0.0; |
423 | | CoinBigIndex start = columnStart[iColumn]; |
424 | | CoinBigIndex end = start + columnLength[iColumn]; |
425 | | for (j = start; j < end; j++) { |
426 | | int jRow = row[j]; |
427 | | value += x[jRow] * elementByColumn[j]; |
428 | | } |
429 | | y[iColumn] += value * scalar; |
430 | | } |
431 | | } |
432 | | } |
433 | | void |
434 | | ClpPackedMatrix::times(double scalar, |
435 | | const double * COIN_RESTRICT x, double * COIN_RESTRICT y, |
436 | | const double * COIN_RESTRICT rowScale, |
437 | | const double * COIN_RESTRICT columnScale) const |
438 | | { |
439 | | if (rowScale) { |
440 | | int iRow, iColumn; |
441 | | // get matrix data pointers |
442 | | const int * COIN_RESTRICT row = matrix_->getIndices(); |
443 | | const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
444 | | const double * COIN_RESTRICT elementByColumn = matrix_->getElements(); |
445 | | if (!(flags_ & 2)) { |
446 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
447 | | double value = x[iColumn]; |
448 | | if (value) { |
449 | | // scaled |
450 | | value *= scalar * columnScale[iColumn]; |
451 | | CoinBigIndex start = columnStart[iColumn]; |
452 | | CoinBigIndex end = columnStart[iColumn+1]; |
453 | | CoinBigIndex j; |
454 | | for (j = start; j < end; j++) { |
455 | | iRow = row[j]; |
456 | | y[iRow] += value * elementByColumn[j] * rowScale[iRow]; |
457 | | } |
458 | | } |
459 | | } |
460 | | } else { |
461 | | const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths(); |
462 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
463 | | double value = x[iColumn]; |
464 | | if (value) { |
465 | | // scaled |
466 | | value *= scalar * columnScale[iColumn]; |
467 | | CoinBigIndex start = columnStart[iColumn]; |
468 | | CoinBigIndex end = start + columnLength[iColumn]; |
469 | | CoinBigIndex j; |
470 | | for (j = start; j < end; j++) { |
471 | | iRow = row[j]; |
472 | | y[iRow] += value * elementByColumn[j] * rowScale[iRow]; |
473 | | } |
474 | | } |
475 | | } |
476 | | } |
477 | | } else { |
478 | | times(scalar, x, y); |
479 | | } |
480 | | } |
481 | | void |
482 | | ClpPackedMatrix::transposeTimes( double scalar, |
483 | | const double * COIN_RESTRICT x, double * COIN_RESTRICT y, |
484 | | const double * COIN_RESTRICT rowScale, |
485 | | const double * COIN_RESTRICT columnScale, |
486 | | double * COIN_RESTRICT spare) const |
487 | | { |
488 | | if (rowScale) { |
489 | | int iColumn; |
490 | | // get matrix data pointers |
491 | | const int * COIN_RESTRICT row = matrix_->getIndices(); |
492 | | const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
493 | | const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths(); |
494 | | const double * COIN_RESTRICT elementByColumn = matrix_->getElements(); |
495 | | if (!spare) { |
496 | | if (!(flags_ & 2)) { |
497 | | CoinBigIndex start = columnStart[0]; |
498 | | if (scalar == -1.0) { |
499 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
500 | | CoinBigIndex j; |
501 | | CoinBigIndex next = columnStart[iColumn+1]; |
502 | | double value = 0.0; |
503 | | // scaled |
504 | | for (j = start; j < next; j++) { |
505 | | int jRow = row[j]; |
506 | | value += x[jRow] * elementByColumn[j] * rowScale[jRow]; |
507 | | } |
508 | | start = next; |
509 | | y[iColumn] -= value * columnScale[iColumn]; |
510 | | } |
511 | | } else { |
512 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
513 | | CoinBigIndex j; |
514 | | CoinBigIndex next = columnStart[iColumn+1]; |
515 | | double value = 0.0; |
516 | | // scaled |
517 | | for (j = start; j < next; j++) { |
518 | | int jRow = row[j]; |
519 | | value += x[jRow] * elementByColumn[j] * rowScale[jRow]; |
520 | | } |
521 | | start = next; |
522 | | y[iColumn] += value * scalar * columnScale[iColumn]; |
523 | | } |
524 | | } |
525 | | } else { |
526 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
527 | | CoinBigIndex j; |
528 | | double value = 0.0; |
529 | | // scaled |
530 | | for (j = columnStart[iColumn]; |
531 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
532 | | int jRow = row[j]; |
533 | | value += x[jRow] * elementByColumn[j] * rowScale[jRow]; |
534 | | } |
535 | | y[iColumn] += value * scalar * columnScale[iColumn]; |
536 | | } |
537 | | } |
538 | | } else { |
539 | | // can use spare region |
540 | | int iRow; |
541 | | int numberRows = matrix_->getNumRows(); |
542 | | for (iRow = 0; iRow < numberRows; iRow++) { |
543 | | double value = x[iRow]; |
544 | | if (value) |
545 | | spare[iRow] = value * rowScale[iRow]; |
546 | | else |
547 | | spare[iRow] = 0.0; |
548 | | } |
549 | | if (!(flags_ & 2)) { |
550 | | CoinBigIndex start = columnStart[0]; |
551 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
552 | | CoinBigIndex j; |
553 | | CoinBigIndex next = columnStart[iColumn+1]; |
554 | | double value = 0.0; |
555 | | // scaled |
556 | | for (j = start; j < next; j++) { |
557 | | int jRow = row[j]; |
558 | | value += spare[jRow] * elementByColumn[j]; |
559 | | } |
560 | | start = next; |
561 | | y[iColumn] += value * scalar * columnScale[iColumn]; |
562 | | } |
563 | | } else { |
564 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
565 | | CoinBigIndex j; |
566 | | double value = 0.0; |
567 | | // scaled |
568 | | for (j = columnStart[iColumn]; |
569 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
570 | | int jRow = row[j]; |
571 | | value += spare[jRow] * elementByColumn[j]; |
572 | | } |
573 | | y[iColumn] += value * scalar * columnScale[iColumn]; |
574 | | } |
575 | | } |
576 | | // no need to zero out |
577 | | //for (iRow=0;iRow<numberRows;iRow++) |
578 | | //spare[iRow] = 0.0; |
579 | | } |
580 | | } else { |
581 | | transposeTimes(scalar, x, y); |
582 | | } |
| 501 | CoinBigIndex j; |
| 502 | CoinBigIndex next = columnStart[iColumn + 1]; |
| 503 | double value = 0.0; |
| 504 | // scaled |
| 505 | for (j = start; j < next; j++) { |
| 506 | int jRow = row[j]; |
| 507 | value += x[jRow] * elementByColumn[j] * rowScale[jRow]; |
| 508 | } |
| 509 | start = next; |
| 510 | y[iColumn] -= value * columnScale[iColumn]; |
| 511 | } |
| 512 | } else { |
| 513 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 514 | CoinBigIndex j; |
| 515 | CoinBigIndex next = columnStart[iColumn + 1]; |
| 516 | double value = 0.0; |
| 517 | // scaled |
| 518 | for (j = start; j < next; j++) { |
| 519 | int jRow = row[j]; |
| 520 | value += x[jRow] * elementByColumn[j] * rowScale[jRow]; |
| 521 | } |
| 522 | start = next; |
| 523 | y[iColumn] += value * scalar * columnScale[iColumn]; |
| 524 | } |
| 525 | } |
| 526 | } else { |
| 527 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 528 | CoinBigIndex j; |
| 529 | double value = 0.0; |
| 530 | // scaled |
| 531 | for (j = columnStart[iColumn]; |
| 532 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 533 | int jRow = row[j]; |
| 534 | value += x[jRow] * elementByColumn[j] * rowScale[jRow]; |
| 535 | } |
| 536 | y[iColumn] += value * scalar * columnScale[iColumn]; |
| 537 | } |
| 538 | } |
| 539 | } else { |
| 540 | // can use spare region |
| 541 | int iRow; |
| 542 | int numberRows = matrix_->getNumRows(); |
| 543 | for (iRow = 0; iRow < numberRows; iRow++) { |
| 544 | double value = x[iRow]; |
| 545 | if (value) |
| 546 | spare[iRow] = value * rowScale[iRow]; |
| 547 | else |
| 548 | spare[iRow] = 0.0; |
| 549 | } |
| 550 | if (!(flags_ & 2)) { |
| 551 | CoinBigIndex start = columnStart[0]; |
| 552 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 553 | CoinBigIndex j; |
| 554 | CoinBigIndex next = columnStart[iColumn + 1]; |
| 555 | double value = 0.0; |
| 556 | // scaled |
| 557 | for (j = start; j < next; j++) { |
| 558 | int jRow = row[j]; |
| 559 | value += spare[jRow] * elementByColumn[j]; |
| 560 | } |
| 561 | start = next; |
| 562 | y[iColumn] += value * scalar * columnScale[iColumn]; |
| 563 | } |
| 564 | } else { |
| 565 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 566 | CoinBigIndex j; |
| 567 | double value = 0.0; |
| 568 | // scaled |
| 569 | for (j = columnStart[iColumn]; |
| 570 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 571 | int jRow = row[j]; |
| 572 | value += spare[jRow] * elementByColumn[j]; |
| 573 | } |
| 574 | y[iColumn] += value * scalar * columnScale[iColumn]; |
| 575 | } |
| 576 | } |
| 577 | // no need to zero out |
| 578 | //for (iRow=0;iRow<numberRows;iRow++) |
| 579 | //spare[iRow] = 0.0; |
| 580 | } |
| 581 | } else { |
| 582 | transposeTimes(scalar, x, y); |
| 583 | } |
611 | | void |
612 | | ClpPackedMatrix::transposeTimesSubset( int number, |
613 | | const int * which, |
614 | | const double * COIN_RESTRICT x, double * COIN_RESTRICT y, |
615 | | const double * COIN_RESTRICT rowScale, |
616 | | const double * COIN_RESTRICT columnScale, |
617 | | double * COIN_RESTRICT spare) const |
618 | | { |
619 | | // get matrix data pointers |
620 | | const int * COIN_RESTRICT row = matrix_->getIndices(); |
621 | | const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
622 | | const double * COIN_RESTRICT elementByColumn = matrix_->getElements(); |
623 | | if (!spare || !rowScale) { |
624 | | if (rowScale) { |
625 | | for (int jColumn = 0; jColumn < number; jColumn++) { |
626 | | int iColumn = which[jColumn]; |
627 | | CoinBigIndex j; |
628 | | CoinBigIndex start = columnStart[iColumn]; |
629 | | CoinBigIndex next = columnStart[iColumn+1]; |
630 | | double value = 0.0; |
631 | | for (j = start; j < next; j++) { |
632 | | int jRow = row[j]; |
633 | | value += x[jRow] * elementByColumn[j] * rowScale[jRow]; |
634 | | } |
635 | | y[iColumn] -= value * columnScale[iColumn]; |
636 | | } |
637 | | } else { |
| 612 | void ClpPackedMatrix::transposeTimesSubset(int number, |
| 613 | const int *which, |
| 614 | const double *COIN_RESTRICT x, double *COIN_RESTRICT y, |
| 615 | const double *COIN_RESTRICT rowScale, |
| 616 | const double *COIN_RESTRICT columnScale, |
| 617 | double *COIN_RESTRICT spare) const |
| 618 | { |
| 619 | // get matrix data pointers |
| 620 | const int *COIN_RESTRICT row = matrix_->getIndices(); |
| 621 | const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
| 622 | const double *COIN_RESTRICT elementByColumn = matrix_->getElements(); |
| 623 | if (!spare || !rowScale) { |
| 624 | if (rowScale) { |
| 625 | for (int jColumn = 0; jColumn < number; jColumn++) { |
| 626 | int iColumn = which[jColumn]; |
| 627 | CoinBigIndex j; |
| 628 | CoinBigIndex start = columnStart[iColumn]; |
| 629 | CoinBigIndex next = columnStart[iColumn + 1]; |
| 630 | double value = 0.0; |
| 631 | for (j = start; j < next; j++) { |
| 632 | int jRow = row[j]; |
| 633 | value += x[jRow] * elementByColumn[j] * rowScale[jRow]; |
| 634 | } |
| 635 | y[iColumn] -= value * columnScale[iColumn]; |
| 636 | } |
| 637 | } else { |
639 | | int numberThreads=abcState(); |
640 | | if (numberThreads) { |
641 | | clpTempInfo info[ABOCA_LITE]; |
642 | | int chunk = (number+numberThreads-1)/numberThreads; |
643 | | int n=0; |
644 | | for (int i=0;i<numberThreads;i++) { |
645 | | info[i].spare=y; |
646 | | info[i].work=const_cast<double *>(x); |
647 | | info[i].which=const_cast<int *>(which); |
648 | | info[i].startColumn=n; |
649 | | info[i].element=elementByColumn; |
650 | | info[i].start=columnStart; |
651 | | info[i].row=row; |
652 | | info[i].numberToDo=CoinMin(chunk,number-n); |
653 | | n += chunk; |
654 | | } |
655 | | for (int i=0;i<numberThreads;i++) { |
656 | | cilk_spawn transposeTimesSubsetBit(info[i]); |
657 | | } |
658 | | cilk_sync; |
659 | | } else { |
660 | | #endif |
661 | | for (int jColumn = 0; jColumn < number; jColumn++) { |
662 | | int iColumn = which[jColumn]; |
663 | | CoinBigIndex j; |
664 | | CoinBigIndex start = columnStart[iColumn]; |
665 | | CoinBigIndex next = columnStart[iColumn+1]; |
666 | | double value = 0.0; |
667 | | for (j = start; j < next; j++) { |
668 | | int jRow = row[j]; |
669 | | value += x[jRow] * elementByColumn[j]; |
670 | | } |
671 | | y[iColumn] -= value; |
672 | | } |
| 639 | int numberThreads = abcState(); |
| 640 | if (numberThreads) { |
| 641 | clpTempInfo info[ABOCA_LITE]; |
| 642 | int chunk = (number + numberThreads - 1) / numberThreads; |
| 643 | int n = 0; |
| 644 | for (int i = 0; i < numberThreads; i++) { |
| 645 | info[i].spare = y; |
| 646 | info[i].work = const_cast< double * >(x); |
| 647 | info[i].which = const_cast< int * >(which); |
| 648 | info[i].startColumn = n; |
| 649 | info[i].element = elementByColumn; |
| 650 | info[i].start = columnStart; |
| 651 | info[i].row = row; |
| 652 | info[i].numberToDo = CoinMin(chunk, number - n); |
| 653 | n += chunk; |
| 654 | } |
| 655 | for (int i = 0; i < numberThreads; i++) { |
| 656 | cilk_spawn transposeTimesSubsetBit(info[i]); |
| 657 | } |
| 658 | cilk_sync; |
| 659 | } else { |
| 660 | #endif |
| 661 | for (int jColumn = 0; jColumn < number; jColumn++) { |
| 662 | int iColumn = which[jColumn]; |
| 663 | CoinBigIndex j; |
| 664 | CoinBigIndex start = columnStart[iColumn]; |
| 665 | CoinBigIndex next = columnStart[iColumn + 1]; |
| 666 | double value = 0.0; |
| 667 | for (j = start; j < next; j++) { |
| 668 | int jRow = row[j]; |
| 669 | value += x[jRow] * elementByColumn[j]; |
| 670 | } |
| 671 | y[iColumn] -= value; |
| 672 | } |
723 | | int numberRows = model->numberRows(); |
724 | | ClpPackedMatrix* rowCopy = |
725 | | static_cast< ClpPackedMatrix*>(model->rowCopy()); |
726 | | bool packed = rowArray->packedMode(); |
727 | | double factor = (numberRows < 100) ? 0.25 : 0.35; |
728 | | factor = 0.5; |
729 | | // We may not want to do by row if there may be cache problems |
730 | | // It would be nice to find L2 cache size - for moment 512K |
731 | | // Be slightly optimistic |
732 | | if (numberActiveColumns_ * sizeof(double) > 1000000) { |
733 | | if (numberRows * 10 < numberActiveColumns_) |
734 | | factor *= 0.333333333; |
735 | | else if (numberRows * 4 < numberActiveColumns_) |
736 | | factor *= 0.5; |
737 | | else if (numberRows * 2 < numberActiveColumns_) |
738 | | factor *= 0.66666666667; |
739 | | //if (model->numberIterations()%50==0) |
740 | | //printf("%d nonzero\n",numberInRowArray); |
741 | | } |
742 | | // if not packed then bias a bit more towards by column |
743 | | if (!packed) |
744 | | factor *= 0.9; |
745 | | assert (!y->getNumElements()); |
746 | | double multiplierX = 0.8; |
747 | | double factor2 = factor * multiplierX; |
748 | | if (packed && rowCopy_ && numberInRowArray > 2 && numberInRowArray > factor2 * numberRows && |
749 | | numberInRowArray < 0.9 * numberRows && 0) { |
750 | | rowCopy_->transposeTimes(model, rowCopy->matrix_, rowArray, y, columnArray); |
751 | | return; |
752 | | } |
753 | | if (columnCopy_) |
754 | | factor *= 0.7; |
755 | | if (numberInRowArray > factor * numberRows || !rowCopy) { |
756 | | // do by column |
757 | | // If no gaps - can do a bit faster |
758 | | if (!(flags_ & 2) || columnCopy_) { |
759 | | transposeTimesByColumn( model, scalar, |
760 | | rowArray, y, columnArray); |
761 | | return; |
762 | | } |
763 | | int iColumn; |
764 | | // get matrix data pointers |
765 | | const int * COIN_RESTRICT row = matrix_->getIndices(); |
766 | | const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
767 | | const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths(); |
768 | | const double * COIN_RESTRICT elementByColumn = matrix_->getElements(); |
769 | | const double * COIN_RESTRICT rowScale = model->rowScale(); |
| 722 | int numberRows = model->numberRows(); |
| 723 | ClpPackedMatrix *rowCopy = static_cast< ClpPackedMatrix * >(model->rowCopy()); |
| 724 | bool packed = rowArray->packedMode(); |
| 725 | double factor = (numberRows < 100) ? 0.25 : 0.35; |
| 726 | factor = 0.5; |
| 727 | // We may not want to do by row if there may be cache problems |
| 728 | // It would be nice to find L2 cache size - for moment 512K |
| 729 | // Be slightly optimistic |
| 730 | if (numberActiveColumns_ * sizeof(double) > 1000000) { |
| 731 | if (numberRows * 10 < numberActiveColumns_) |
| 732 | factor *= 0.333333333; |
| 733 | else if (numberRows * 4 < numberActiveColumns_) |
| 734 | factor *= 0.5; |
| 735 | else if (numberRows * 2 < numberActiveColumns_) |
| 736 | factor *= 0.66666666667; |
| 737 | //if (model->numberIterations()%50==0) |
| 738 | //printf("%d nonzero\n",numberInRowArray); |
| 739 | } |
| 740 | // if not packed then bias a bit more towards by column |
| 741 | if (!packed) |
| 742 | factor *= 0.9; |
| 743 | assert(!y->getNumElements()); |
| 744 | double multiplierX = 0.8; |
| 745 | double factor2 = factor * multiplierX; |
| 746 | if (packed && rowCopy_ && numberInRowArray > 2 && numberInRowArray > factor2 * numberRows && numberInRowArray < 0.9 * numberRows && 0) { |
| 747 | rowCopy_->transposeTimes(model, rowCopy->matrix_, rowArray, y, columnArray); |
| 748 | return; |
| 749 | } |
| 750 | if (columnCopy_) |
| 751 | factor *= 0.7; |
| 752 | if (numberInRowArray > factor * numberRows || !rowCopy) { |
| 753 | // do by column |
| 754 | // If no gaps - can do a bit faster |
| 755 | if (!(flags_ & 2) || columnCopy_) { |
| 756 | transposeTimesByColumn(model, scalar, |
| 757 | rowArray, y, columnArray); |
| 758 | return; |
| 759 | } |
| 760 | int iColumn; |
| 761 | // get matrix data pointers |
| 762 | const int *COIN_RESTRICT row = matrix_->getIndices(); |
| 763 | const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
| 764 | const int *COIN_RESTRICT columnLength = matrix_->getVectorLengths(); |
| 765 | const double *COIN_RESTRICT elementByColumn = matrix_->getElements(); |
| 766 | const double *COIN_RESTRICT rowScale = model->rowScale(); |
781 | | if (packed) { |
782 | | // need to expand pi into y |
783 | | assert(y->capacity() >= numberRows); |
784 | | double * COIN_RESTRICT piOld = pi; |
785 | | pi = y->denseVector(); |
786 | | const int * COIN_RESTRICT whichRow = rowArray->getIndices(); |
787 | | int i; |
788 | | if (!rowScale) { |
789 | | // modify pi so can collapse to one loop |
790 | | if (scalar == -1.0) { |
791 | | for (i = 0; i < numberInRowArray; i++) { |
792 | | int iRow = whichRow[i]; |
793 | | pi[iRow] = -piOld[i]; |
794 | | } |
795 | | } else { |
796 | | for (i = 0; i < numberInRowArray; i++) { |
797 | | int iRow = whichRow[i]; |
798 | | pi[iRow] = scalar * piOld[i]; |
799 | | } |
800 | | } |
801 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
802 | | double value = 0.0; |
803 | | CoinBigIndex j; |
804 | | for (j = columnStart[iColumn]; |
805 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
806 | | int iRow = row[j]; |
807 | | value += pi[iRow] * elementByColumn[j]; |
808 | | } |
809 | | if (fabs(value) > zeroTolerance) { |
810 | | array[numberNonZero] = value; |
811 | | index[numberNonZero++] = iColumn; |
812 | | } |
813 | | } |
814 | | } else { |
| 778 | if (packed) { |
| 779 | // need to expand pi into y |
| 780 | assert(y->capacity() >= numberRows); |
| 781 | double *COIN_RESTRICT piOld = pi; |
| 782 | pi = y->denseVector(); |
| 783 | const int *COIN_RESTRICT whichRow = rowArray->getIndices(); |
| 784 | int i; |
| 785 | if (!rowScale) { |
| 786 | // modify pi so can collapse to one loop |
| 787 | if (scalar == -1.0) { |
| 788 | for (i = 0; i < numberInRowArray; i++) { |
| 789 | int iRow = whichRow[i]; |
| 790 | pi[iRow] = -piOld[i]; |
| 791 | } |
| 792 | } else { |
| 793 | for (i = 0; i < numberInRowArray; i++) { |
| 794 | int iRow = whichRow[i]; |
| 795 | pi[iRow] = scalar * piOld[i]; |
| 796 | } |
| 797 | } |
| 798 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 799 | double value = 0.0; |
| 800 | CoinBigIndex j; |
| 801 | for (j = columnStart[iColumn]; |
| 802 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 803 | int iRow = row[j]; |
| 804 | value += pi[iRow] * elementByColumn[j]; |
| 805 | } |
| 806 | if (fabs(value) > zeroTolerance) { |
| 807 | array[numberNonZero] = value; |
| 808 | index[numberNonZero++] = iColumn; |
| 809 | } |
| 810 | } |
| 811 | } else { |
816 | | if (model->clpScaledMatrix()) |
817 | | printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__); |
818 | | #endif |
819 | | // scaled |
820 | | // modify pi so can collapse to one loop |
821 | | if (scalar == -1.0) { |
822 | | for (i = 0; i < numberInRowArray; i++) { |
823 | | int iRow = whichRow[i]; |
824 | | pi[iRow] = -piOld[i] * rowScale[iRow]; |
825 | | } |
826 | | } else { |
827 | | for (i = 0; i < numberInRowArray; i++) { |
828 | | int iRow = whichRow[i]; |
829 | | pi[iRow] = scalar * piOld[i] * rowScale[iRow]; |
830 | | } |
831 | | } |
832 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
833 | | double value = 0.0; |
834 | | CoinBigIndex j; |
835 | | const double * columnScale = model->columnScale(); |
836 | | for (j = columnStart[iColumn]; |
837 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
838 | | int iRow = row[j]; |
839 | | value += pi[iRow] * elementByColumn[j]; |
840 | | } |
841 | | value *= columnScale[iColumn]; |
842 | | if (fabs(value) > zeroTolerance) { |
843 | | array[numberNonZero] = value; |
844 | | index[numberNonZero++] = iColumn; |
845 | | } |
846 | | } |
847 | | } |
848 | | // zero out |
849 | | for (i = 0; i < numberInRowArray; i++) { |
850 | | int iRow = whichRow[i]; |
851 | | pi[iRow] = 0.0; |
852 | | } |
853 | | } else { |
854 | | if (!rowScale) { |
855 | | if (scalar == -1.0) { |
856 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
857 | | double value = 0.0; |
858 | | CoinBigIndex j; |
859 | | for (j = columnStart[iColumn]; |
860 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
861 | | int iRow = row[j]; |
862 | | value += pi[iRow] * elementByColumn[j]; |
863 | | } |
864 | | if (fabs(value) > zeroTolerance) { |
865 | | index[numberNonZero++] = iColumn; |
866 | | array[iColumn] = -value; |
867 | | } |
868 | | } |
869 | | } else { |
870 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
871 | | double value = 0.0; |
872 | | CoinBigIndex j; |
873 | | for (j = columnStart[iColumn]; |
874 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
875 | | int iRow = row[j]; |
876 | | value += pi[iRow] * elementByColumn[j]; |
877 | | } |
878 | | value *= scalar; |
879 | | if (fabs(value) > zeroTolerance) { |
880 | | index[numberNonZero++] = iColumn; |
881 | | array[iColumn] = value; |
882 | | } |
883 | | } |
884 | | } |
885 | | } else { |
| 813 | if (model->clpScaledMatrix()) |
| 814 | printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__); |
| 815 | #endif |
| 816 | // scaled |
| 817 | // modify pi so can collapse to one loop |
| 818 | if (scalar == -1.0) { |
| 819 | for (i = 0; i < numberInRowArray; i++) { |
| 820 | int iRow = whichRow[i]; |
| 821 | pi[iRow] = -piOld[i] * rowScale[iRow]; |
| 822 | } |
| 823 | } else { |
| 824 | for (i = 0; i < numberInRowArray; i++) { |
| 825 | int iRow = whichRow[i]; |
| 826 | pi[iRow] = scalar * piOld[i] * rowScale[iRow]; |
| 827 | } |
| 828 | } |
| 829 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 830 | double value = 0.0; |
| 831 | CoinBigIndex j; |
| 832 | const double *columnScale = model->columnScale(); |
| 833 | for (j = columnStart[iColumn]; |
| 834 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 835 | int iRow = row[j]; |
| 836 | value += pi[iRow] * elementByColumn[j]; |
| 837 | } |
| 838 | value *= columnScale[iColumn]; |
| 839 | if (fabs(value) > zeroTolerance) { |
| 840 | array[numberNonZero] = value; |
| 841 | index[numberNonZero++] = iColumn; |
| 842 | } |
| 843 | } |
| 844 | } |
| 845 | // zero out |
| 846 | for (i = 0; i < numberInRowArray; i++) { |
| 847 | int iRow = whichRow[i]; |
| 848 | pi[iRow] = 0.0; |
| 849 | } |
| 850 | } else { |
| 851 | if (!rowScale) { |
| 852 | if (scalar == -1.0) { |
| 853 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 854 | double value = 0.0; |
| 855 | CoinBigIndex j; |
| 856 | for (j = columnStart[iColumn]; |
| 857 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 858 | int iRow = row[j]; |
| 859 | value += pi[iRow] * elementByColumn[j]; |
| 860 | } |
| 861 | if (fabs(value) > zeroTolerance) { |
| 862 | index[numberNonZero++] = iColumn; |
| 863 | array[iColumn] = -value; |
| 864 | } |
| 865 | } |
| 866 | } else { |
| 867 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 868 | double value = 0.0; |
| 869 | CoinBigIndex j; |
| 870 | for (j = columnStart[iColumn]; |
| 871 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 872 | int iRow = row[j]; |
| 873 | value += pi[iRow] * elementByColumn[j]; |
| 874 | } |
| 875 | value *= scalar; |
| 876 | if (fabs(value) > zeroTolerance) { |
| 877 | index[numberNonZero++] = iColumn; |
| 878 | array[iColumn] = value; |
| 879 | } |
| 880 | } |
| 881 | } |
| 882 | } else { |
887 | | if (model->clpScaledMatrix()) |
888 | | printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__); |
889 | | #endif |
890 | | // scaled |
891 | | if (scalar == -1.0) { |
892 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
893 | | double value = 0.0; |
894 | | CoinBigIndex j; |
895 | | const double * columnScale = model->columnScale(); |
896 | | for (j = columnStart[iColumn]; |
897 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
898 | | int iRow = row[j]; |
899 | | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
900 | | } |
901 | | value *= columnScale[iColumn]; |
902 | | if (fabs(value) > zeroTolerance) { |
903 | | index[numberNonZero++] = iColumn; |
904 | | array[iColumn] = -value; |
905 | | } |
906 | | } |
907 | | } else { |
908 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
909 | | double value = 0.0; |
910 | | CoinBigIndex j; |
911 | | const double * columnScale = model->columnScale(); |
912 | | for (j = columnStart[iColumn]; |
913 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
914 | | int iRow = row[j]; |
915 | | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
916 | | } |
917 | | value *= scalar * columnScale[iColumn]; |
918 | | if (fabs(value) > zeroTolerance) { |
919 | | index[numberNonZero++] = iColumn; |
920 | | array[iColumn] = value; |
921 | | } |
922 | | } |
923 | | } |
924 | | } |
925 | | } |
926 | | columnArray->setNumElements(numberNonZero); |
927 | | y->setNumElements(0); |
928 | | } else { |
929 | | // do by row |
930 | | rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray); |
931 | | } |
932 | | if (packed) |
933 | | columnArray->setPackedMode(true); |
934 | | if (0) { |
935 | | columnArray->checkClean(); |
936 | | int numberNonZero = columnArray->getNumElements();; |
937 | | int * index = columnArray->getIndices(); |
938 | | double * array = columnArray->denseVector(); |
939 | | int i; |
940 | | for (i = 0; i < numberNonZero; i++) { |
941 | | int j = index[i]; |
942 | | double value; |
943 | | if (packed) |
944 | | value = array[i]; |
945 | | else |
946 | | value = array[j]; |
947 | | printf("Ti %d %d %g\n", i, j, value); |
948 | | } |
949 | | } |
| 884 | if (model->clpScaledMatrix()) |
| 885 | printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__); |
| 886 | #endif |
| 887 | // scaled |
| 888 | if (scalar == -1.0) { |
| 889 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 890 | double value = 0.0; |
| 891 | CoinBigIndex j; |
| 892 | const double *columnScale = model->columnScale(); |
| 893 | for (j = columnStart[iColumn]; |
| 894 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 895 | int iRow = row[j]; |
| 896 | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
| 897 | } |
| 898 | value *= columnScale[iColumn]; |
| 899 | if (fabs(value) > zeroTolerance) { |
| 900 | index[numberNonZero++] = iColumn; |
| 901 | array[iColumn] = -value; |
| 902 | } |
| 903 | } |
| 904 | } else { |
| 905 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 906 | double value = 0.0; |
| 907 | CoinBigIndex j; |
| 908 | const double *columnScale = model->columnScale(); |
| 909 | for (j = columnStart[iColumn]; |
| 910 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 911 | int iRow = row[j]; |
| 912 | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
| 913 | } |
| 914 | value *= scalar * columnScale[iColumn]; |
| 915 | if (fabs(value) > zeroTolerance) { |
| 916 | index[numberNonZero++] = iColumn; |
| 917 | array[iColumn] = value; |
| 918 | } |
| 919 | } |
| 920 | } |
| 921 | } |
| 922 | } |
| 923 | columnArray->setNumElements(numberNonZero); |
| 924 | y->setNumElements(0); |
| 925 | } else { |
| 926 | // do by row |
| 927 | rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray); |
| 928 | } |
| 929 | if (packed) |
| 930 | columnArray->setPackedMode(true); |
| 931 | if (0) { |
| 932 | columnArray->checkClean(); |
| 933 | int numberNonZero = columnArray->getNumElements(); |
| 934 | ; |
| 935 | int *index = columnArray->getIndices(); |
| 936 | double *array = columnArray->denseVector(); |
| 937 | int i; |
| 938 | for (i = 0; i < numberNonZero; i++) { |
| 939 | int j = index[i]; |
| 940 | double value; |
| 941 | if (packed) |
| 942 | value = array[i]; |
| 943 | else |
| 944 | value = array[j]; |
| 945 | printf("Ti %d %d %g\n", i, j, value); |
| 946 | } |
| 947 | } |
961 | | void |
962 | | ClpPackedMatrix::transposeTimesByColumn(const ClpSimplex * model, double scalar, |
963 | | const CoinIndexedVector * rowArray, |
964 | | CoinIndexedVector * y, |
965 | | CoinIndexedVector * columnArray) const |
966 | | { |
967 | | double * COIN_RESTRICT pi = rowArray->denseVector(); |
968 | | int numberNonZero = 0; |
969 | | int * COIN_RESTRICT index = columnArray->getIndices(); |
970 | | double * COIN_RESTRICT array = columnArray->denseVector(); |
971 | | int numberInRowArray = rowArray->getNumElements(); |
972 | | // maybe I need one in OsiSimplex |
973 | | double zeroTolerance = model->zeroTolerance(); |
974 | | bool packed = rowArray->packedMode(); |
975 | | // do by column |
976 | | int iColumn; |
977 | | // get matrix data pointers |
978 | | const int * COIN_RESTRICT row = matrix_->getIndices(); |
979 | | const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
980 | | const double * COIN_RESTRICT elementByColumn = matrix_->getElements(); |
981 | | const double * COIN_RESTRICT rowScale = model->rowScale(); |
982 | | assert (!y->getNumElements()); |
983 | | assert (numberActiveColumns_ > 0); |
984 | | const ClpPackedMatrix * thisMatrix = this; |
| 959 | void ClpPackedMatrix::transposeTimesByColumn(const ClpSimplex *model, double scalar, |
| 960 | const CoinIndexedVector *rowArray, |
| 961 | CoinIndexedVector *y, |
| 962 | CoinIndexedVector *columnArray) const |
| 963 | { |
| 964 | double *COIN_RESTRICT pi = rowArray->denseVector(); |
| 965 | int numberNonZero = 0; |
| 966 | int *COIN_RESTRICT index = columnArray->getIndices(); |
| 967 | double *COIN_RESTRICT array = columnArray->denseVector(); |
| 968 | int numberInRowArray = rowArray->getNumElements(); |
| 969 | // maybe I need one in OsiSimplex |
| 970 | double zeroTolerance = model->zeroTolerance(); |
| 971 | bool packed = rowArray->packedMode(); |
| 972 | // do by column |
| 973 | int iColumn; |
| 974 | // get matrix data pointers |
| 975 | const int *COIN_RESTRICT row = matrix_->getIndices(); |
| 976 | const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
| 977 | const double *COIN_RESTRICT elementByColumn = matrix_->getElements(); |
| 978 | const double *COIN_RESTRICT rowScale = model->rowScale(); |
| 979 | assert(!y->getNumElements()); |
| 980 | assert(numberActiveColumns_ > 0); |
| 981 | const ClpPackedMatrix *thisMatrix = this; |
1001 | | if (packed) { |
1002 | | // need to expand pi into y |
1003 | | assert(y->capacity() >= model->numberRows()); |
1004 | | double * COIN_RESTRICT piOld = pi; |
1005 | | pi = y->denseVector(); |
1006 | | const int * COIN_RESTRICT whichRow = rowArray->getIndices(); |
1007 | | int i; |
1008 | | if (!rowScale) { |
1009 | | // modify pi so can collapse to one loop |
1010 | | if (scalar == -1.0) { |
1011 | | //yA += numberInRowArray; |
1012 | | for (i = 0; i < numberInRowArray; i++) { |
1013 | | int iRow = whichRow[i]; |
1014 | | pi[iRow] = -piOld[i]; |
| 998 | if (packed) { |
| 999 | // need to expand pi into y |
| 1000 | assert(y->capacity() >= model->numberRows()); |
| 1001 | double *COIN_RESTRICT piOld = pi; |
| 1002 | pi = y->denseVector(); |
| 1003 | const int *COIN_RESTRICT whichRow = rowArray->getIndices(); |
| 1004 | int i; |
| 1005 | if (!rowScale) { |
| 1006 | // modify pi so can collapse to one loop |
| 1007 | if (scalar == -1.0) { |
| 1008 | //yA += numberInRowArray; |
| 1009 | for (i = 0; i < numberInRowArray; i++) { |
| 1010 | int iRow = whichRow[i]; |
| 1011 | pi[iRow] = -piOld[i]; |
| 1012 | } |
| 1013 | } else { |
| 1014 | for (i = 0; i < numberInRowArray; i++) { |
| 1015 | int iRow = whichRow[i]; |
| 1016 | pi[iRow] = scalar * piOld[i]; |
| 1017 | } |
| 1018 | } |
| 1019 | if (!columnCopy_) { |
| 1020 | if ((model->specialOptions(), 131072) != 0) { |
| 1021 | if (model->spareIntArray_[0] > 0) { |
| 1022 | CoinIndexedVector *spareArray = model->rowArray(3); |
| 1023 | // also do dualColumn stuff |
| 1024 | double *COIN_RESTRICT spare = spareArray->denseVector(); |
| 1025 | int *COIN_RESTRICT spareIndex = spareArray->getIndices(); |
| 1026 | const double *COIN_RESTRICT reducedCost = model->djRegion(0); |
| 1027 | double multiplier[] = { -1.0, 1.0 }; |
| 1028 | double dualT = -model->currentDualTolerance(); |
| 1029 | double acceptablePivot = model->spareDoubleArray_[0]; |
| 1030 | // We can also see if infeasible or pivoting on free |
| 1031 | double tentativeTheta = 1.0e15; |
| 1032 | double upperTheta = 1.0e31; |
| 1033 | int addSequence = model->numberColumns(); |
| 1034 | const unsigned char *COIN_RESTRICT statusArray = model->statusArray() + addSequence; |
| 1035 | int numberRemaining = 0; |
| 1036 | assert(scalar == -1.0); |
| 1037 | for (i = 0; i < numberInRowArray; i++) { |
| 1038 | int iSequence = whichRow[i]; |
| 1039 | int iStatus = (statusArray[iSequence] & 3) - 1; |
| 1040 | if (iStatus) { |
| 1041 | double mult = multiplier[iStatus - 1]; |
| 1042 | double alpha = piOld[i] * mult; |
| 1043 | double oldValue; |
| 1044 | double value; |
| 1045 | if (alpha > 0.0) { |
| 1046 | oldValue = reducedCost[iSequence] * mult; |
| 1047 | value = oldValue - tentativeTheta * alpha; |
| 1048 | if (value < dualT) { |
| 1049 | value = oldValue - upperTheta * alpha; |
| 1050 | if (value < dualT && alpha >= acceptablePivot) { |
| 1051 | upperTheta = (oldValue - dualT) / alpha; |
| 1052 | //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta); |
1016 | | } else { |
1017 | | for (i = 0; i < numberInRowArray; i++) { |
1018 | | int iRow = whichRow[i]; |
1019 | | pi[iRow] = scalar * piOld[i]; |
1020 | | } |
1021 | | } |
1022 | | if (!columnCopy_) { |
1023 | | if ((model->specialOptions(), 131072) != 0) { |
1024 | | if(model->spareIntArray_[0] > 0) { |
1025 | | CoinIndexedVector * spareArray = model->rowArray(3); |
1026 | | // also do dualColumn stuff |
1027 | | double * COIN_RESTRICT spare = spareArray->denseVector(); |
1028 | | int * COIN_RESTRICT spareIndex = spareArray->getIndices(); |
1029 | | const double * COIN_RESTRICT reducedCost = model->djRegion(0); |
1030 | | double multiplier[] = { -1.0, 1.0}; |
1031 | | double dualT = - model->currentDualTolerance(); |
1032 | | double acceptablePivot = model->spareDoubleArray_[0]; |
1033 | | // We can also see if infeasible or pivoting on free |
1034 | | double tentativeTheta = 1.0e15; |
1035 | | double upperTheta = 1.0e31; |
1036 | | int addSequence = model->numberColumns(); |
1037 | | const unsigned char * COIN_RESTRICT statusArray = model->statusArray() + addSequence; |
1038 | | int numberRemaining = 0; |
1039 | | assert (scalar == -1.0); |
1040 | | for (i = 0; i < numberInRowArray; i++) { |
1041 | | int iSequence = whichRow[i]; |
1042 | | int iStatus = (statusArray[iSequence] & 3) - 1; |
1043 | | if (iStatus) { |
1044 | | double mult = multiplier[iStatus-1]; |
1045 | | double alpha = piOld[i] * mult; |
1046 | | double oldValue; |
1047 | | double value; |
1048 | | if (alpha > 0.0) { |
1049 | | oldValue = reducedCost[iSequence] * mult; |
1050 | | value = oldValue - tentativeTheta * alpha; |
1051 | | if (value < dualT) { |
1052 | | value = oldValue - upperTheta * alpha; |
1053 | | if (value < dualT && alpha >= acceptablePivot) { |
1054 | | upperTheta = (oldValue - dualT) / alpha; |
1055 | | //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta); |
1056 | | } |
1057 | | // add to list |
1058 | | spare[numberRemaining] = alpha * mult; |
1059 | | spareIndex[numberRemaining++] = iSequence + addSequence; |
1060 | | } |
1061 | | } |
1062 | | } |
1063 | | } |
1064 | | numberNonZero = |
1065 | | thisMatrix->gutsOfTransposeTimesUnscaled(pi, |
1066 | | columnArray->getIndices(), |
1067 | | columnArray->denseVector(), |
1068 | | model->statusArray(), |
1069 | | spareIndex, |
1070 | | spare, |
1071 | | model->djRegion(1), |
1072 | | upperTheta, |
1073 | | acceptablePivot, |
1074 | | model->currentDualTolerance(), |
1075 | | numberRemaining, |
1076 | | zeroTolerance); |
1077 | | model->spareDoubleArray_[0] = upperTheta; |
1078 | | spareArray->setNumElements(numberRemaining); |
| 1054 | // add to list |
| 1055 | spare[numberRemaining] = alpha * mult; |
| 1056 | spareIndex[numberRemaining++] = iSequence + addSequence; |
| 1057 | } |
| 1058 | } |
| 1059 | } |
| 1060 | } |
| 1061 | numberNonZero = thisMatrix->gutsOfTransposeTimesUnscaled(pi, |
| 1062 | columnArray->getIndices(), |
| 1063 | columnArray->denseVector(), |
| 1064 | model->statusArray(), |
| 1065 | spareIndex, |
| 1066 | spare, |
| 1067 | model->djRegion(1), |
| 1068 | upperTheta, |
| 1069 | acceptablePivot, |
| 1070 | model->currentDualTolerance(), |
| 1071 | numberRemaining, |
| 1072 | zeroTolerance); |
| 1073 | model->spareDoubleArray_[0] = upperTheta; |
| 1074 | spareArray->setNumElements(numberRemaining); |
1122 | | if (model->clpScaledMatrix()) |
1123 | | printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__); |
1124 | | #endif |
1125 | | // scaled |
1126 | | // modify pi so can collapse to one loop |
1127 | | if (scalar == -1.0) { |
1128 | | //yC += numberInRowArray; |
1129 | | for (i = 0; i < numberInRowArray; i++) { |
1130 | | int iRow = whichRow[i]; |
1131 | | pi[iRow] = -piOld[i] * rowScale[iRow]; |
1132 | | } |
1133 | | } else { |
1134 | | for (i = 0; i < numberInRowArray; i++) { |
1135 | | int iRow = whichRow[i]; |
1136 | | pi[iRow] = scalar * piOld[i] * rowScale[iRow]; |
1137 | | } |
1138 | | } |
1139 | | const double * columnScale = model->columnScale(); |
1140 | | if (!columnCopy_) { |
1141 | | if ((model->specialOptions(), 131072) != 0) |
1142 | | numberNonZero = |
1143 | | gutsOfTransposeTimesScaled(pi, columnScale, |
1144 | | columnArray->getIndices(), |
1145 | | columnArray->denseVector(), |
1146 | | model->statusArray(), |
1147 | | zeroTolerance); |
1148 | | else |
1149 | | numberNonZero = |
1150 | | gutsOfTransposeTimesScaled(pi, columnScale, |
1151 | | columnArray->getIndices(), |
1152 | | columnArray->denseVector(), |
1153 | | zeroTolerance); |
1154 | | columnArray->setNumElements(numberNonZero); |
1155 | | //xC++; |
1156 | | } else { |
1157 | | if ((model->moreSpecialOptions()&8)!=0 && |
1158 | | model->algorithm()<0) { |
1159 | | columnCopy_->transposeTimes(model, pi, columnArray, |
1160 | | model->rowArray(3),rowArray); |
1161 | | model->spareIntArray_[0]=-2; |
1162 | | } else { |
1163 | | columnCopy_->transposeTimes(model, pi, columnArray); |
1164 | | } |
1165 | | numberNonZero = columnArray->getNumElements(); |
1166 | | //xD++; |
1167 | | } |
1168 | | } |
1169 | | // zero out |
1170 | | int numberRows = model->numberRows(); |
1171 | | if (numberInRowArray * 4 < numberRows) { |
1172 | | for (i = 0; i < numberInRowArray; i++) { |
1173 | | int iRow = whichRow[i]; |
1174 | | pi[iRow] = 0.0; |
1175 | | } |
1176 | | } else { |
1177 | | CoinZeroN(pi, numberRows); |
1178 | | } |
1179 | | //int kA=xA+xB; |
1180 | | //int kC=xC+xD; |
1181 | | //if ((kA+kC)%10000==0) |
1182 | | //printf("AA %d %d %g, CC %d %d %g\n", |
1183 | | // xA,xB,kA ? yA/(double)(kA): 0.0,xC,xD,kC ? yC/(double) (kC) :0.0); |
1184 | | } else { |
1185 | | if (!rowScale) { |
1186 | | if (scalar == -1.0) { |
1187 | | double value = 0.0; |
1188 | | CoinBigIndex j; |
1189 | | CoinBigIndex end = columnStart[1]; |
1190 | | for (j = columnStart[0]; j < end; j++) { |
1191 | | int iRow = row[j]; |
1192 | | value += pi[iRow] * elementByColumn[j]; |
1193 | | } |
1194 | | for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) { |
1195 | | CoinBigIndex start = end; |
1196 | | end = columnStart[iColumn+2]; |
1197 | | if (fabs(value) > zeroTolerance) { |
1198 | | array[iColumn] = -value; |
1199 | | index[numberNonZero++] = iColumn; |
1200 | | } |
1201 | | value = 0.0; |
1202 | | for (j = start; j < end; j++) { |
1203 | | int iRow = row[j]; |
1204 | | value += pi[iRow] * elementByColumn[j]; |
1205 | | } |
1206 | | } |
1207 | | if (fabs(value) > zeroTolerance) { |
1208 | | array[iColumn] = -value; |
1209 | | index[numberNonZero++] = iColumn; |
1210 | | } |
1211 | | } else { |
1212 | | double value = 0.0; |
1213 | | CoinBigIndex j; |
1214 | | CoinBigIndex end = columnStart[1]; |
1215 | | for (j = columnStart[0]; j < end; j++) { |
1216 | | int iRow = row[j]; |
1217 | | value += pi[iRow] * elementByColumn[j]; |
1218 | | } |
1219 | | for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) { |
1220 | | value *= scalar; |
1221 | | CoinBigIndex start = end; |
1222 | | end = columnStart[iColumn+2]; |
1223 | | if (fabs(value) > zeroTolerance) { |
1224 | | array[iColumn] = value; |
1225 | | index[numberNonZero++] = iColumn; |
1226 | | } |
1227 | | value = 0.0; |
1228 | | for (j = start; j < end; j++) { |
1229 | | int iRow = row[j]; |
1230 | | value += pi[iRow] * elementByColumn[j]; |
1231 | | } |
1232 | | } |
1233 | | value *= scalar; |
1234 | | if (fabs(value) > zeroTolerance) { |
1235 | | array[iColumn] = value; |
1236 | | index[numberNonZero++] = iColumn; |
1237 | | } |
1238 | | } |
1239 | | } else { |
| 1115 | if (model->clpScaledMatrix()) |
| 1116 | printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__); |
| 1117 | #endif |
| 1118 | // scaled |
| 1119 | // modify pi so can collapse to one loop |
| 1120 | if (scalar == -1.0) { |
| 1121 | //yC += numberInRowArray; |
| 1122 | for (i = 0; i < numberInRowArray; i++) { |
| 1123 | int iRow = whichRow[i]; |
| 1124 | pi[iRow] = -piOld[i] * rowScale[iRow]; |
| 1125 | } |
| 1126 | } else { |
| 1127 | for (i = 0; i < numberInRowArray; i++) { |
| 1128 | int iRow = whichRow[i]; |
| 1129 | pi[iRow] = scalar * piOld[i] * rowScale[iRow]; |
| 1130 | } |
| 1131 | } |
| 1132 | const double *columnScale = model->columnScale(); |
| 1133 | if (!columnCopy_) { |
| 1134 | if ((model->specialOptions(), 131072) != 0) |
| 1135 | numberNonZero = gutsOfTransposeTimesScaled(pi, columnScale, |
| 1136 | columnArray->getIndices(), |
| 1137 | columnArray->denseVector(), |
| 1138 | model->statusArray(), |
| 1139 | zeroTolerance); |
| 1140 | else |
| 1141 | numberNonZero = gutsOfTransposeTimesScaled(pi, columnScale, |
| 1142 | columnArray->getIndices(), |
| 1143 | columnArray->denseVector(), |
| 1144 | zeroTolerance); |
| 1145 | columnArray->setNumElements(numberNonZero); |
| 1146 | //xC++; |
| 1147 | } else { |
| 1148 | if ((model->moreSpecialOptions() & 8) != 0 && model->algorithm() < 0) { |
| 1149 | columnCopy_->transposeTimes(model, pi, columnArray, |
| 1150 | model->rowArray(3), rowArray); |
| 1151 | model->spareIntArray_[0] = -2; |
| 1152 | } else { |
| 1153 | columnCopy_->transposeTimes(model, pi, columnArray); |
| 1154 | } |
| 1155 | numberNonZero = columnArray->getNumElements(); |
| 1156 | //xD++; |
| 1157 | } |
| 1158 | } |
| 1159 | // zero out |
| 1160 | int numberRows = model->numberRows(); |
| 1161 | if (numberInRowArray * 4 < numberRows) { |
| 1162 | for (i = 0; i < numberInRowArray; i++) { |
| 1163 | int iRow = whichRow[i]; |
| 1164 | pi[iRow] = 0.0; |
| 1165 | } |
| 1166 | } else { |
| 1167 | CoinZeroN(pi, numberRows); |
| 1168 | } |
| 1169 | //int kA=xA+xB; |
| 1170 | //int kC=xC+xD; |
| 1171 | //if ((kA+kC)%10000==0) |
| 1172 | //printf("AA %d %d %g, CC %d %d %g\n", |
| 1173 | // xA,xB,kA ? yA/(double)(kA): 0.0,xC,xD,kC ? yC/(double) (kC) :0.0); |
| 1174 | } else { |
| 1175 | if (!rowScale) { |
| 1176 | if (scalar == -1.0) { |
| 1177 | double value = 0.0; |
| 1178 | CoinBigIndex j; |
| 1179 | CoinBigIndex end = columnStart[1]; |
| 1180 | for (j = columnStart[0]; j < end; j++) { |
| 1181 | int iRow = row[j]; |
| 1182 | value += pi[iRow] * elementByColumn[j]; |
| 1183 | } |
| 1184 | for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) { |
| 1185 | CoinBigIndex start = end; |
| 1186 | end = columnStart[iColumn + 2]; |
| 1187 | if (fabs(value) > zeroTolerance) { |
| 1188 | array[iColumn] = -value; |
| 1189 | index[numberNonZero++] = iColumn; |
| 1190 | } |
| 1191 | value = 0.0; |
| 1192 | for (j = start; j < end; j++) { |
| 1193 | int iRow = row[j]; |
| 1194 | value += pi[iRow] * elementByColumn[j]; |
| 1195 | } |
| 1196 | } |
| 1197 | if (fabs(value) > zeroTolerance) { |
| 1198 | array[iColumn] = -value; |
| 1199 | index[numberNonZero++] = iColumn; |
| 1200 | } |
| 1201 | } else { |
| 1202 | double value = 0.0; |
| 1203 | CoinBigIndex j; |
| 1204 | CoinBigIndex end = columnStart[1]; |
| 1205 | for (j = columnStart[0]; j < end; j++) { |
| 1206 | int iRow = row[j]; |
| 1207 | value += pi[iRow] * elementByColumn[j]; |
| 1208 | } |
| 1209 | for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) { |
| 1210 | value *= scalar; |
| 1211 | CoinBigIndex start = end; |
| 1212 | end = columnStart[iColumn + 2]; |
| 1213 | if (fabs(value) > zeroTolerance) { |
| 1214 | array[iColumn] = value; |
| 1215 | index[numberNonZero++] = iColumn; |
| 1216 | } |
| 1217 | value = 0.0; |
| 1218 | for (j = start; j < end; j++) { |
| 1219 | int iRow = row[j]; |
| 1220 | value += pi[iRow] * elementByColumn[j]; |
| 1221 | } |
| 1222 | } |
| 1223 | value *= scalar; |
| 1224 | if (fabs(value) > zeroTolerance) { |
| 1225 | array[iColumn] = value; |
| 1226 | index[numberNonZero++] = iColumn; |
| 1227 | } |
| 1228 | } |
| 1229 | } else { |
1241 | | if (model->clpScaledMatrix()) |
1242 | | printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__); |
1243 | | #endif |
1244 | | // scaled |
1245 | | if (scalar == -1.0) { |
1246 | | const double * COIN_RESTRICT columnScale = model->columnScale(); |
1247 | | double value = 0.0; |
1248 | | double scale = columnScale[0]; |
1249 | | CoinBigIndex j; |
1250 | | CoinBigIndex end = columnStart[1]; |
1251 | | for (j = columnStart[0]; j < end; j++) { |
1252 | | int iRow = row[j]; |
1253 | | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
1254 | | } |
1255 | | for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) { |
1256 | | value *= scale; |
1257 | | CoinBigIndex start = end; |
1258 | | end = columnStart[iColumn+2]; |
1259 | | scale = columnScale[iColumn+1]; |
1260 | | if (fabs(value) > zeroTolerance) { |
1261 | | array[iColumn] = -value; |
1262 | | index[numberNonZero++] = iColumn; |
1263 | | } |
1264 | | value = 0.0; |
1265 | | for (j = start; j < end; j++) { |
1266 | | int iRow = row[j]; |
1267 | | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
1268 | | } |
1269 | | } |
1270 | | value *= scale; |
1271 | | if (fabs(value) > zeroTolerance) { |
1272 | | array[iColumn] = -value; |
1273 | | index[numberNonZero++] = iColumn; |
1274 | | } |
1275 | | } else { |
1276 | | const double * COIN_RESTRICT columnScale = model->columnScale(); |
1277 | | double value = 0.0; |
1278 | | double scale = columnScale[0] * scalar; |
1279 | | CoinBigIndex j; |
1280 | | CoinBigIndex end = columnStart[1]; |
1281 | | for (j = columnStart[0]; j < end; j++) { |
1282 | | int iRow = row[j]; |
1283 | | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
1284 | | } |
1285 | | for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) { |
1286 | | value *= scale; |
1287 | | CoinBigIndex start = end; |
1288 | | end = columnStart[iColumn+2]; |
1289 | | scale = columnScale[iColumn+1] * scalar; |
1290 | | if (fabs(value) > zeroTolerance) { |
1291 | | array[iColumn] = value; |
1292 | | index[numberNonZero++] = iColumn; |
1293 | | } |
1294 | | value = 0.0; |
1295 | | for (j = start; j < end; j++) { |
1296 | | int iRow = row[j]; |
1297 | | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
1298 | | } |
1299 | | } |
1300 | | value *= scale; |
1301 | | if (fabs(value) > zeroTolerance) { |
1302 | | array[iColumn] = value; |
1303 | | index[numberNonZero++] = iColumn; |
1304 | | } |
1305 | | } |
1306 | | } |
1307 | | } |
1308 | | columnArray->setNumElements(numberNonZero); |
1309 | | y->setNumElements(0); |
1310 | | if (packed) |
1311 | | columnArray->setPackedMode(true); |
| 1231 | if (model->clpScaledMatrix()) |
| 1232 | printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__); |
| 1233 | #endif |
| 1234 | // scaled |
| 1235 | if (scalar == -1.0) { |
| 1236 | const double *COIN_RESTRICT columnScale = model->columnScale(); |
| 1237 | double value = 0.0; |
| 1238 | double scale = columnScale[0]; |
| 1239 | CoinBigIndex j; |
| 1240 | CoinBigIndex end = columnStart[1]; |
| 1241 | for (j = columnStart[0]; j < end; j++) { |
| 1242 | int iRow = row[j]; |
| 1243 | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
| 1244 | } |
| 1245 | for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) { |
| 1246 | value *= scale; |
| 1247 | CoinBigIndex start = end; |
| 1248 | end = columnStart[iColumn + 2]; |
| 1249 | scale = columnScale[iColumn + 1]; |
| 1250 | if (fabs(value) > zeroTolerance) { |
| 1251 | array[iColumn] = -value; |
| 1252 | index[numberNonZero++] = iColumn; |
| 1253 | } |
| 1254 | value = 0.0; |
| 1255 | for (j = start; j < end; j++) { |
| 1256 | int iRow = row[j]; |
| 1257 | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
| 1258 | } |
| 1259 | } |
| 1260 | value *= scale; |
| 1261 | if (fabs(value) > zeroTolerance) { |
| 1262 | array[iColumn] = -value; |
| 1263 | index[numberNonZero++] = iColumn; |
| 1264 | } |
| 1265 | } else { |
| 1266 | const double *COIN_RESTRICT columnScale = model->columnScale(); |
| 1267 | double value = 0.0; |
| 1268 | double scale = columnScale[0] * scalar; |
| 1269 | CoinBigIndex j; |
| 1270 | CoinBigIndex end = columnStart[1]; |
| 1271 | for (j = columnStart[0]; j < end; j++) { |
| 1272 | int iRow = row[j]; |
| 1273 | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
| 1274 | } |
| 1275 | for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) { |
| 1276 | value *= scale; |
| 1277 | CoinBigIndex start = end; |
| 1278 | end = columnStart[iColumn + 2]; |
| 1279 | scale = columnScale[iColumn + 1] * scalar; |
| 1280 | if (fabs(value) > zeroTolerance) { |
| 1281 | array[iColumn] = value; |
| 1282 | index[numberNonZero++] = iColumn; |
| 1283 | } |
| 1284 | value = 0.0; |
| 1285 | for (j = start; j < end; j++) { |
| 1286 | int iRow = row[j]; |
| 1287 | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
| 1288 | } |
| 1289 | } |
| 1290 | value *= scale; |
| 1291 | if (fabs(value) > zeroTolerance) { |
| 1292 | array[iColumn] = value; |
| 1293 | index[numberNonZero++] = iColumn; |
| 1294 | } |
| 1295 | } |
| 1296 | } |
| 1297 | } |
| 1298 | columnArray->setNumElements(numberNonZero); |
| 1299 | y->setNumElements(0); |
| 1300 | if (packed) |
| 1301 | columnArray->setPackedMode(true); |
1315 | | void |
1316 | | ClpPackedMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar, |
1317 | | const CoinIndexedVector * rowArray, |
1318 | | CoinIndexedVector * y, |
1319 | | CoinIndexedVector * columnArray) const |
1320 | | { |
1321 | | columnArray->clear(); |
1322 | | double * COIN_RESTRICT pi = rowArray->denseVector(); |
1323 | | int numberNonZero = 0; |
1324 | | int * COIN_RESTRICT index = columnArray->getIndices(); |
1325 | | double * COIN_RESTRICT array = columnArray->denseVector(); |
1326 | | int numberInRowArray = rowArray->getNumElements(); |
1327 | | // maybe I need one in OsiSimplex |
1328 | | double zeroTolerance = model->zeroTolerance(); |
1329 | | const int * COIN_RESTRICT column = matrix_->getIndices(); |
1330 | | const CoinBigIndex * COIN_RESTRICT rowStart = getVectorStarts(); |
1331 | | const double * COIN_RESTRICT element = getElements(); |
1332 | | const int * COIN_RESTRICT whichRow = rowArray->getIndices(); |
1333 | | bool packed = rowArray->packedMode(); |
1334 | | if (numberInRowArray > 2) { |
1335 | | // do by rows |
1336 | | // ** Row copy is already scaled |
1337 | | int iRow; |
1338 | | int i; |
1339 | | int numberOriginal = 0; |
1340 | | if (packed) { |
1341 | | int * COIN_RESTRICT index = columnArray->getIndices(); |
1342 | | double * COIN_RESTRICT array = columnArray->denseVector(); |
| 1305 | void ClpPackedMatrix::transposeTimesByRow(const ClpSimplex *model, double scalar, |
| 1306 | const CoinIndexedVector *rowArray, |
| 1307 | CoinIndexedVector *y, |
| 1308 | CoinIndexedVector *columnArray) const |
| 1309 | { |
| 1310 | columnArray->clear(); |
| 1311 | double *COIN_RESTRICT pi = rowArray->denseVector(); |
| 1312 | int numberNonZero = 0; |
| 1313 | int *COIN_RESTRICT index = columnArray->getIndices(); |
| 1314 | double *COIN_RESTRICT array = columnArray->denseVector(); |
| 1315 | int numberInRowArray = rowArray->getNumElements(); |
| 1316 | // maybe I need one in OsiSimplex |
| 1317 | double zeroTolerance = model->zeroTolerance(); |
| 1318 | const int *COIN_RESTRICT column = matrix_->getIndices(); |
| 1319 | const CoinBigIndex *COIN_RESTRICT rowStart = getVectorStarts(); |
| 1320 | const double *COIN_RESTRICT element = getElements(); |
| 1321 | const int *COIN_RESTRICT whichRow = rowArray->getIndices(); |
| 1322 | bool packed = rowArray->packedMode(); |
| 1323 | if (numberInRowArray > 2) { |
| 1324 | // do by rows |
| 1325 | // ** Row copy is already scaled |
| 1326 | int iRow; |
| 1327 | int i; |
| 1328 | int numberOriginal = 0; |
| 1329 | if (packed) { |
| 1330 | int *COIN_RESTRICT index = columnArray->getIndices(); |
| 1331 | double *COIN_RESTRICT array = columnArray->denseVector(); |
1380 | | double * COIN_RESTRICT array2 = y->denseVector(); |
1381 | | numberNonZero=gutsOfTransposeTimesByRowGE3(rowArray,index,array, |
1382 | | array2,zeroTolerance,scalar); |
1383 | | #endif |
1384 | | } else { |
1385 | | numberNonZero = gutsOfTransposeTimesByRowGEK(rowArray, index, array, |
1386 | | numberColumns, zeroTolerance, scalar); |
1387 | | } |
1388 | | columnArray->setNumElements(numberNonZero); |
1389 | | } else { |
1390 | | double * COIN_RESTRICT markVector = y->denseVector(); |
1391 | | numberNonZero = 0; |
1392 | | // and set up mark as char array |
1393 | | char * COIN_RESTRICT marked = reinterpret_cast<char *> (markVector); |
1394 | | for (i = 0; i < numberOriginal; i++) { |
1395 | | int iColumn = index[i]; |
1396 | | marked[iColumn] = 0; |
1397 | | } |
| 1369 | double *COIN_RESTRICT array2 = y->denseVector(); |
| 1370 | numberNonZero = gutsOfTransposeTimesByRowGE3(rowArray, index, array, |
| 1371 | array2, zeroTolerance, scalar); |
| 1372 | #endif |
| 1373 | } else { |
| 1374 | numberNonZero = gutsOfTransposeTimesByRowGEK(rowArray, index, array, |
| 1375 | numberColumns, zeroTolerance, scalar); |
| 1376 | } |
| 1377 | columnArray->setNumElements(numberNonZero); |
| 1378 | } else { |
| 1379 | double *COIN_RESTRICT markVector = y->denseVector(); |
| 1380 | numberNonZero = 0; |
| 1381 | // and set up mark as char array |
| 1382 | char *COIN_RESTRICT marked = reinterpret_cast< char * >(markVector); |
| 1383 | for (i = 0; i < numberOriginal; i++) { |
| 1384 | int iColumn = index[i]; |
| 1385 | marked[iColumn] = 0; |
| 1386 | } |
1432 | | double value; |
1433 | | if (packed) { |
1434 | | gutsOfTransposeTimesByRowEQ2(rowArray, columnArray, y, zeroTolerance, scalar); |
1435 | | numberNonZero = columnArray->getNumElements(); |
1436 | | } else { |
1437 | | int iRow = whichRow[0]; |
1438 | | value = pi[iRow] * scalar; |
1439 | | for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) { |
1440 | | int iColumn = column[j]; |
1441 | | double value2 = value * element[j]; |
1442 | | index[numberNonZero++] = iColumn; |
1443 | | array[iColumn] = value2; |
1444 | | } |
1445 | | iRow = whichRow[1]; |
1446 | | value = pi[iRow] * scalar; |
1447 | | for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) { |
1448 | | int iColumn = column[j]; |
1449 | | double value2 = value * element[j]; |
1450 | | // I am assuming no zeros in matrix |
1451 | | if (array[iColumn]) |
1452 | | value2 += array[iColumn]; |
1453 | | else |
1454 | | index[numberNonZero++] = iColumn; |
1455 | | array[iColumn] = value2; |
1456 | | } |
1457 | | // get rid of tiny values and zero out marked |
1458 | | numberOriginal = numberNonZero; |
1459 | | numberNonZero = 0; |
1460 | | for (i = 0; i < numberOriginal; i++) { |
1461 | | int iColumn = index[i]; |
1462 | | if (fabs(array[iColumn]) > zeroTolerance) { |
1463 | | index[numberNonZero++] = iColumn; |
1464 | | } else { |
1465 | | array[iColumn] = 0.0; |
1466 | | } |
1467 | | } |
1468 | | } |
1469 | | } else if (numberInRowArray == 1) { |
1470 | | // Just one row |
1471 | | int iRow = rowArray->getIndices()[0]; |
1472 | | numberNonZero = 0; |
1473 | | CoinBigIndex j; |
1474 | | if (packed) { |
1475 | | gutsOfTransposeTimesByRowEQ1(rowArray, columnArray, zeroTolerance, scalar); |
1476 | | numberNonZero = columnArray->getNumElements(); |
1477 | | } else { |
1478 | | double value = pi[iRow] * scalar; |
1479 | | for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) { |
1480 | | int iColumn = column[j]; |
1481 | | double value2 = value * element[j]; |
1482 | | if (fabs(value2) > zeroTolerance) { |
1483 | | index[numberNonZero++] = iColumn; |
1484 | | array[iColumn] = value2; |
1485 | | } |
1486 | | } |
1487 | | } |
1488 | | } |
1489 | | columnArray->setNumElements(numberNonZero); |
1490 | | y->setNumElements(0); |
| 1421 | double value; |
| 1422 | if (packed) { |
| 1423 | gutsOfTransposeTimesByRowEQ2(rowArray, columnArray, y, zeroTolerance, scalar); |
| 1424 | numberNonZero = columnArray->getNumElements(); |
| 1425 | } else { |
| 1426 | int iRow = whichRow[0]; |
| 1427 | value = pi[iRow] * scalar; |
| 1428 | for (j = rowStart[iRow]; j < rowStart[iRow + 1]; j++) { |
| 1429 | int iColumn = column[j]; |
| 1430 | double value2 = value * element[j]; |
| 1431 | index[numberNonZero++] = iColumn; |
| 1432 | array[iColumn] = value2; |
| 1433 | } |
| 1434 | iRow = whichRow[1]; |
| 1435 | value = pi[iRow] * scalar; |
| 1436 | for (j = rowStart[iRow]; j < rowStart[iRow + 1]; j++) { |
| 1437 | int iColumn = column[j]; |
| 1438 | double value2 = value * element[j]; |
| 1439 | // I am assuming no zeros in matrix |
| 1440 | if (array[iColumn]) |
| 1441 | value2 += array[iColumn]; |
| 1442 | else |
| 1443 | index[numberNonZero++] = iColumn; |
| 1444 | array[iColumn] = value2; |
| 1445 | } |
| 1446 | // get rid of tiny values and zero out marked |
| 1447 | numberOriginal = numberNonZero; |
| 1448 | numberNonZero = 0; |
| 1449 | for (i = 0; i < numberOriginal; i++) { |
| 1450 | int iColumn = index[i]; |
| 1451 | if (fabs(array[iColumn]) > zeroTolerance) { |
| 1452 | index[numberNonZero++] = iColumn; |
| 1453 | } else { |
| 1454 | array[iColumn] = 0.0; |
| 1455 | } |
| 1456 | } |
| 1457 | } |
| 1458 | } else if (numberInRowArray == 1) { |
| 1459 | // Just one row |
| 1460 | int iRow = rowArray->getIndices()[0]; |
| 1461 | numberNonZero = 0; |
| 1462 | CoinBigIndex j; |
| 1463 | if (packed) { |
| 1464 | gutsOfTransposeTimesByRowEQ1(rowArray, columnArray, zeroTolerance, scalar); |
| 1465 | numberNonZero = columnArray->getNumElements(); |
| 1466 | } else { |
| 1467 | double value = pi[iRow] * scalar; |
| 1468 | for (j = rowStart[iRow]; j < rowStart[iRow + 1]; j++) { |
| 1469 | int iColumn = column[j]; |
| 1470 | double value2 = value * element[j]; |
| 1471 | if (fabs(value2) > zeroTolerance) { |
| 1472 | index[numberNonZero++] = iColumn; |
| 1473 | array[iColumn] = value2; |
| 1474 | } |
| 1475 | } |
| 1476 | } |
| 1477 | } |
| 1478 | columnArray->setNumElements(numberNonZero); |
| 1479 | y->setNumElements(0); |
1553 | | int |
1554 | | ClpPackedMatrix::gutsOfTransposeTimesScaled(const double * COIN_RESTRICT pi, |
1555 | | const double * COIN_RESTRICT columnScale, |
1556 | | int * COIN_RESTRICT index, |
1557 | | double * COIN_RESTRICT array, |
1558 | | const double zeroTolerance) const |
1559 | | { |
1560 | | int numberNonZero = 0; |
1561 | | // get matrix data pointers |
1562 | | const int * COIN_RESTRICT row = matrix_->getIndices(); |
1563 | | const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
1564 | | const double * COIN_RESTRICT elementByColumn = matrix_->getElements(); |
1565 | | double value = 0.0; |
1566 | | double scale = columnScale[0]; |
1567 | | CoinBigIndex j; |
1568 | | CoinBigIndex end = columnStart[1]; |
1569 | | for (j = columnStart[0]; j < end; j++) { |
1570 | | int iRow = row[j]; |
1571 | | value += pi[iRow] * elementByColumn[j]; |
1572 | | } |
1573 | | int iColumn; |
1574 | | for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) { |
1575 | | value *= scale; |
1576 | | CoinBigIndex start = end; |
1577 | | scale = columnScale[iColumn+1]; |
1578 | | end = columnStart[iColumn+2]; |
1579 | | if (fabs(value) > zeroTolerance) { |
1580 | | array[numberNonZero] = value; |
1581 | | index[numberNonZero++] = iColumn; |
1582 | | } |
1583 | | value = 0.0; |
1584 | | for (j = start; j < end; j++) { |
1585 | | int iRow = row[j]; |
1586 | | value += pi[iRow] * elementByColumn[j]; |
1587 | | } |
1588 | | } |
1589 | | value *= scale; |
1590 | | if (fabs(value) > zeroTolerance) { |
1591 | | array[numberNonZero] = value; |
1592 | | index[numberNonZero++] = iColumn; |
1593 | | } |
1594 | | return numberNonZero; |
| 1541 | int ClpPackedMatrix::gutsOfTransposeTimesScaled(const double *COIN_RESTRICT pi, |
| 1542 | const double *COIN_RESTRICT columnScale, |
| 1543 | int *COIN_RESTRICT index, |
| 1544 | double *COIN_RESTRICT array, |
| 1545 | const double zeroTolerance) const |
| 1546 | { |
| 1547 | int numberNonZero = 0; |
| 1548 | // get matrix data pointers |
| 1549 | const int *COIN_RESTRICT row = matrix_->getIndices(); |
| 1550 | const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
| 1551 | const double *COIN_RESTRICT elementByColumn = matrix_->getElements(); |
| 1552 | double value = 0.0; |
| 1553 | double scale = columnScale[0]; |
| 1554 | CoinBigIndex j; |
| 1555 | CoinBigIndex end = columnStart[1]; |
| 1556 | for (j = columnStart[0]; j < end; j++) { |
| 1557 | int iRow = row[j]; |
| 1558 | value += pi[iRow] * elementByColumn[j]; |
| 1559 | } |
| 1560 | int iColumn; |
| 1561 | for (iColumn = 0; iColumn < numberActiveColumns_ - 1; iColumn++) { |
| 1562 | value *= scale; |
| 1563 | CoinBigIndex start = end; |
| 1564 | scale = columnScale[iColumn + 1]; |
| 1565 | end = columnStart[iColumn + 2]; |
| 1566 | if (fabs(value) > zeroTolerance) { |
| 1567 | array[numberNonZero] = value; |
| 1568 | index[numberNonZero++] = iColumn; |
| 1569 | } |
| 1570 | value = 0.0; |
| 1571 | for (j = start; j < end; j++) { |
| 1572 | int iRow = row[j]; |
| 1573 | value += pi[iRow] * elementByColumn[j]; |
| 1574 | } |
| 1575 | } |
| 1576 | value *= scale; |
| 1577 | if (fabs(value) > zeroTolerance) { |
| 1578 | array[numberNonZero] = value; |
| 1579 | index[numberNonZero++] = iColumn; |
| 1580 | } |
| 1581 | return numberNonZero; |
1664 | | int numberThreads=abcState(); |
1665 | | if (numberThreads) { |
1666 | | clpTempInfo info[ABOCA_LITE]; |
1667 | | int chunk = (numberActiveColumns_+numberThreads-1)/numberThreads; |
1668 | | int n=0; |
1669 | | for (int i=0;i<numberThreads;i++) { |
1670 | | info[i].which = index+n; |
1671 | | info[i].infeas = array+n; |
1672 | | info[i].work=const_cast<double *>(pi); |
1673 | | info[i].status=const_cast<unsigned char *>(status); |
1674 | | info[i].startColumn=n; |
1675 | | info[i].element=elementByColumn; |
1676 | | info[i].start=columnStart; |
1677 | | info[i].row=row; |
1678 | | info[i].numberToDo=CoinMin(chunk,numberActiveColumns_-n); |
1679 | | info[i].tolerance=zeroTolerance; |
1680 | | n += chunk; |
1681 | | } |
1682 | | for (int i=0;i<numberThreads;i++) { |
1683 | | cilk_spawn transposeTimesUnscaledBit(info[i]); |
1684 | | } |
1685 | | cilk_sync; |
1686 | | for (int i=0;i<numberThreads;i++) |
1687 | | numberNonZero += info[i].numberAdded; |
1688 | | moveAndZero(info,2,NULL); |
1689 | | } else { |
1690 | | #endif |
1691 | | double value = 0.0; |
1692 | | int jColumn = -1; |
1693 | | for (int iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
1694 | | bool wanted = ((status[iColumn] & 3) != 1); |
1695 | | if (fabs(value) > zeroTolerance) { |
1696 | | array[numberNonZero] = value; |
1697 | | index[numberNonZero++] = jColumn; |
1698 | | } |
1699 | | value = 0.0; |
1700 | | if (wanted) { |
1701 | | CoinBigIndex start = columnStart[iColumn]; |
1702 | | CoinBigIndex end = columnStart[iColumn+1]; |
1703 | | jColumn = iColumn; |
1704 | | int n = static_cast<int>(end - start); |
1705 | | bool odd = (n & 1) != 0; |
1706 | | n = n >> 1; |
1707 | | const int * COIN_RESTRICT rowThis = row + start; |
1708 | | const double * COIN_RESTRICT elementThis = elementByColumn + start; |
1709 | | for (; n; n--) { |
1710 | | int iRow0 = *rowThis; |
1711 | | int iRow1 = *(rowThis + 1); |
1712 | | rowThis += 2; |
1713 | | value += pi[iRow0] * (*elementThis); |
1714 | | value += pi[iRow1] * (*(elementThis + 1)); |
1715 | | elementThis += 2; |
1716 | | } |
1717 | | if (odd) { |
1718 | | int iRow = *rowThis; |
1719 | | value += pi[iRow] * (*elementThis); |
1720 | | } |
1721 | | } |
1722 | | } |
1723 | | if (fabs(value) > zeroTolerance) { |
1724 | | array[numberNonZero] = value; |
1725 | | index[numberNonZero++] = jColumn; |
1726 | | } |
| 1650 | int numberThreads = abcState(); |
| 1651 | if (numberThreads) { |
| 1652 | clpTempInfo info[ABOCA_LITE]; |
| 1653 | int chunk = (numberActiveColumns_ + numberThreads - 1) / numberThreads; |
| 1654 | int n = 0; |
| 1655 | for (int i = 0; i < numberThreads; i++) { |
| 1656 | info[i].which = index + n; |
| 1657 | info[i].infeas = array + n; |
| 1658 | info[i].work = const_cast< double * >(pi); |
| 1659 | info[i].status = const_cast< unsigned char * >(status); |
| 1660 | info[i].startColumn = n; |
| 1661 | info[i].element = elementByColumn; |
| 1662 | info[i].start = columnStart; |
| 1663 | info[i].row = row; |
| 1664 | info[i].numberToDo = CoinMin(chunk, numberActiveColumns_ - n); |
| 1665 | info[i].tolerance = zeroTolerance; |
| 1666 | n += chunk; |
| 1667 | } |
| 1668 | for (int i = 0; i < numberThreads; i++) { |
| 1669 | cilk_spawn transposeTimesUnscaledBit(info[i]); |
| 1670 | } |
| 1671 | cilk_sync; |
| 1672 | for (int i = 0; i < numberThreads; i++) |
| 1673 | numberNonZero += info[i].numberAdded; |
| 1674 | moveAndZero(info, 2, NULL); |
| 1675 | } else { |
| 1676 | #endif |
| 1677 | double value = 0.0; |
| 1678 | int jColumn = -1; |
| 1679 | for (int iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 1680 | bool wanted = ((status[iColumn] & 3) != 1); |
| 1681 | if (fabs(value) > zeroTolerance) { |
| 1682 | array[numberNonZero] = value; |
| 1683 | index[numberNonZero++] = jColumn; |
| 1684 | } |
| 1685 | value = 0.0; |
| 1686 | if (wanted) { |
| 1687 | CoinBigIndex start = columnStart[iColumn]; |
| 1688 | CoinBigIndex end = columnStart[iColumn + 1]; |
| 1689 | jColumn = iColumn; |
| 1690 | int n = static_cast< int >(end - start); |
| 1691 | bool odd = (n & 1) != 0; |
| 1692 | n = n >> 1; |
| 1693 | const int *COIN_RESTRICT rowThis = row + start; |
| 1694 | const double *COIN_RESTRICT elementThis = elementByColumn + start; |
| 1695 | for (; n; n--) { |
| 1696 | int iRow0 = *rowThis; |
| 1697 | int iRow1 = *(rowThis + 1); |
| 1698 | rowThis += 2; |
| 1699 | value += pi[iRow0] * (*elementThis); |
| 1700 | value += pi[iRow1] * (*(elementThis + 1)); |
| 1701 | elementThis += 2; |
| 1702 | } |
| 1703 | if (odd) { |
| 1704 | int iRow = *rowThis; |
| 1705 | value += pi[iRow] * (*elementThis); |
| 1706 | } |
| 1707 | } |
| 1708 | } |
| 1709 | if (fabs(value) > zeroTolerance) { |
| 1710 | array[numberNonZero] = value; |
| 1711 | index[numberNonZero++] = jColumn; |
| 1712 | } |
1781 | | double mult = multiplier[wanted-1]; |
1782 | | double alpha = value * mult; |
1783 | | array[numberNonZero] = value; |
1784 | | index[numberNonZero++] = iColumn; |
1785 | | if (alpha > 0.0) { |
1786 | | double oldValue = reducedCost[iColumn] * mult; |
1787 | | double value = oldValue - tentativeTheta * alpha; |
1788 | | if (value < dualT) { |
1789 | | value = oldValue - upperTheta * alpha; |
1790 | | if (value < dualT && alpha >= acceptablePivot) { |
1791 | | upperTheta = (oldValue - dualT) / alpha; |
1792 | | //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta); |
1793 | | } |
1794 | | // add to list |
1795 | | spareArray[numberRemaining] = alpha * mult; |
1796 | | spareIndex[numberRemaining++] = iColumn; |
1797 | | } |
1798 | | } |
1799 | | } |
1800 | | } |
1801 | | } |
1802 | | info.numberAdded=numberNonZero; |
1803 | | info.numberRemaining=numberRemaining; |
1804 | | info.upperTheta=upperTheta; |
| 1767 | double mult = multiplier[wanted - 1]; |
| 1768 | double alpha = value * mult; |
| 1769 | array[numberNonZero] = value; |
| 1770 | index[numberNonZero++] = iColumn; |
| 1771 | if (alpha > 0.0) { |
| 1772 | double oldValue = reducedCost[iColumn] * mult; |
| 1773 | double value = oldValue - tentativeTheta * alpha; |
| 1774 | if (value < dualT) { |
| 1775 | value = oldValue - upperTheta * alpha; |
| 1776 | if (value < dualT && alpha >= acceptablePivot) { |
| 1777 | upperTheta = (oldValue - dualT) / alpha; |
| 1778 | //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta); |
| 1779 | } |
| 1780 | // add to list |
| 1781 | spareArray[numberRemaining] = alpha * mult; |
| 1782 | spareIndex[numberRemaining++] = iColumn; |
| 1783 | } |
| 1784 | } |
| 1785 | } |
| 1786 | } |
| 1787 | } |
| 1788 | info.numberAdded = numberNonZero; |
| 1789 | info.numberRemaining = numberRemaining; |
| 1790 | info.upperTheta = upperTheta; |
1809 | | int |
1810 | | ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi, |
1811 | | int * COIN_RESTRICT index, |
1812 | | double * COIN_RESTRICT array, |
1813 | | const unsigned char * COIN_RESTRICT status, |
1814 | | int * COIN_RESTRICT spareIndex, |
1815 | | double * COIN_RESTRICT spareArray, |
1816 | | const double * COIN_RESTRICT reducedCost, |
1817 | | double & upperThetaP, |
1818 | | double acceptablePivot, |
1819 | | double dualTolerance, |
1820 | | int & numberRemainingP, |
1821 | | const double zeroTolerance) const |
1822 | | { |
1823 | | int numberRemaining = numberRemainingP; |
1824 | | double upperTheta = upperThetaP; |
1825 | | int numberNonZero = 0; |
1826 | | // get matrix data pointers |
1827 | | const int * COIN_RESTRICT row = matrix_->getIndices(); |
1828 | | const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
1829 | | const double * COIN_RESTRICT elementByColumn = matrix_->getElements(); |
| 1795 | int ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double *COIN_RESTRICT pi, |
| 1796 | int *COIN_RESTRICT index, |
| 1797 | double *COIN_RESTRICT array, |
| 1798 | const unsigned char *COIN_RESTRICT status, |
| 1799 | int *COIN_RESTRICT spareIndex, |
| 1800 | double *COIN_RESTRICT spareArray, |
| 1801 | const double *COIN_RESTRICT reducedCost, |
| 1802 | double &upperThetaP, |
| 1803 | double acceptablePivot, |
| 1804 | double dualTolerance, |
| 1805 | int &numberRemainingP, |
| 1806 | const double zeroTolerance) const |
| 1807 | { |
| 1808 | int numberRemaining = numberRemainingP; |
| 1809 | double upperTheta = upperThetaP; |
| 1810 | int numberNonZero = 0; |
| 1811 | // get matrix data pointers |
| 1812 | const int *COIN_RESTRICT row = matrix_->getIndices(); |
| 1813 | const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
| 1814 | const double *COIN_RESTRICT elementByColumn = matrix_->getElements(); |
1831 | | int numberThreads=abcState(); |
1832 | | if (numberThreads) { |
1833 | | clpTempInfo info[ABOCA_LITE]; |
1834 | | int chunk = (numberActiveColumns_+numberThreads-1)/numberThreads; |
1835 | | int n=0; |
1836 | | for (int i=0;i<numberThreads;i++) { |
1837 | | info[i].upperTheta = upperThetaP; |
1838 | | info[i].acceptablePivot = acceptablePivot; |
1839 | | info[i].reducedCost = const_cast<double *>(reducedCost); |
1840 | | info[i].index = spareIndex+n+numberRemaining; |
1841 | | info[i].spare = spareArray+n+numberRemaining; |
1842 | | info[i].which = index+n; |
1843 | | info[i].infeas = array+n; |
1844 | | info[i].work=const_cast<double *>(pi); |
1845 | | info[i].status=const_cast<unsigned char *>(status); |
1846 | | info[i].startColumn=n; |
1847 | | info[i].element=elementByColumn; |
1848 | | info[i].start=columnStart; |
1849 | | info[i].row=row; |
1850 | | info[i].numberToDo=CoinMin(chunk,numberActiveColumns_-n); |
1851 | | info[i].tolerance=zeroTolerance; |
1852 | | info[i].dualTolerance=dualTolerance; |
1853 | | n += chunk; |
1854 | | } |
1855 | | for (int i=0;i<numberThreads;i++) { |
1856 | | cilk_spawn transposeTimesUnscaledBit2(info[i]); |
1857 | | } |
1858 | | cilk_sync; |
1859 | | for (int i=0;i<numberThreads;i++) { |
1860 | | numberNonZero += info[i].numberAdded; |
1861 | | numberRemaining += info[i].numberRemaining; |
1862 | | upperTheta = CoinMin(upperTheta,static_cast<double>(info[i].upperTheta)); |
1863 | | } |
1864 | | moveAndZero(info,1,NULL); |
1865 | | moveAndZero(info,2,NULL); |
1866 | | } else { |
1867 | | #endif |
1868 | | double tentativeTheta = 1.0e15; |
1869 | | double multiplier[] = { -1.0, 1.0}; |
1870 | | double dualT = - dualTolerance; |
1871 | | for (int iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
1872 | | int wanted = (status[iColumn] & 3) - 1; |
1873 | | if (wanted) { |
1874 | | double value = 0.0; |
1875 | | CoinBigIndex start = columnStart[iColumn]; |
1876 | | CoinBigIndex end = columnStart[iColumn+1]; |
1877 | | int n = static_cast<int>(end - start); |
| 1816 | int numberThreads = abcState(); |
| 1817 | if (numberThreads) { |
| 1818 | clpTempInfo info[ABOCA_LITE]; |
| 1819 | int chunk = (numberActiveColumns_ + numberThreads - 1) / numberThreads; |
| 1820 | int n = 0; |
| 1821 | for (int i = 0; i < numberThreads; i++) { |
| 1822 | info[i].upperTheta = upperThetaP; |
| 1823 | info[i].acceptablePivot = acceptablePivot; |
| 1824 | info[i].reducedCost = const_cast< double * >(reducedCost); |
| 1825 | info[i].index = spareIndex + n + numberRemaining; |
| 1826 | info[i].spare = spareArray + n + numberRemaining; |
| 1827 | info[i].which = index + n; |
| 1828 | info[i].infeas = array + n; |
| 1829 | info[i].work = const_cast< double * >(pi); |
| 1830 | info[i].status = const_cast< unsigned char * >(status); |
| 1831 | info[i].startColumn = n; |
| 1832 | info[i].element = elementByColumn; |
| 1833 | info[i].start = columnStart; |
| 1834 | info[i].row = row; |
| 1835 | info[i].numberToDo = CoinMin(chunk, numberActiveColumns_ - n); |
| 1836 | info[i].tolerance = zeroTolerance; |
| 1837 | info[i].dualTolerance = dualTolerance; |
| 1838 | n += chunk; |
| 1839 | } |
| 1840 | for (int i = 0; i < numberThreads; i++) { |
| 1841 | cilk_spawn transposeTimesUnscaledBit2(info[i]); |
| 1842 | } |
| 1843 | cilk_sync; |
| 1844 | for (int i = 0; i < numberThreads; i++) { |
| 1845 | numberNonZero += info[i].numberAdded; |
| 1846 | numberRemaining += info[i].numberRemaining; |
| 1847 | upperTheta = CoinMin(upperTheta, static_cast< double >(info[i].upperTheta)); |
| 1848 | } |
| 1849 | moveAndZero(info, 1, NULL); |
| 1850 | moveAndZero(info, 2, NULL); |
| 1851 | } else { |
| 1852 | #endif |
| 1853 | double tentativeTheta = 1.0e15; |
| 1854 | double multiplier[] = { -1.0, 1.0 }; |
| 1855 | double dualT = -dualTolerance; |
| 1856 | for (int iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 1857 | int wanted = (status[iColumn] & 3) - 1; |
| 1858 | if (wanted) { |
| 1859 | double value = 0.0; |
| 1860 | CoinBigIndex start = columnStart[iColumn]; |
| 1861 | CoinBigIndex end = columnStart[iColumn + 1]; |
| 1862 | int n = static_cast< int >(end - start); |
1896 | | const int * COIN_RESTRICT rowThis = &row[end-16]; |
1897 | | const double * COIN_RESTRICT elementThis = &elementByColumn[end-16]; |
1898 | | bool odd = (n & 1) != 0; |
1899 | | n = n >> 1; |
1900 | | double value2 = 0.0; |
1901 | | if (odd) { |
1902 | | int iRow = row[start]; |
1903 | | value2 = pi[iRow] * elementByColumn[start]; |
1904 | | } |
1905 | | switch (n) { |
1906 | | default: { |
1907 | | if (odd) { |
1908 | | start++; |
1909 | | } |
1910 | | n -= 8; |
1911 | | for (; n; n--) { |
1912 | | int iRow0 = row[start]; |
1913 | | int iRow1 = row[start+1]; |
1914 | | value += pi[iRow0] * elementByColumn[start]; |
1915 | | value2 += pi[iRow1] * elementByColumn[start+1]; |
1916 | | start += 2; |
1917 | | } |
1918 | | case 8: { |
1919 | | int iRow0 = rowThis[16-16]; |
1920 | | int iRow1 = rowThis[16-15]; |
1921 | | value += pi[iRow0] * elementThis[16-16]; |
1922 | | value2 += pi[iRow1] * elementThis[16-15]; |
1923 | | } |
1924 | | case 7: { |
1925 | | int iRow0 = rowThis[16-14]; |
1926 | | int iRow1 = rowThis[16-13]; |
1927 | | value += pi[iRow0] * elementThis[16-14]; |
1928 | | value2 += pi[iRow1] * elementThis[16-13]; |
1929 | | } |
1930 | | case 6: { |
1931 | | int iRow0 = rowThis[16-12]; |
1932 | | int iRow1 = rowThis[16-11]; |
1933 | | value += pi[iRow0] * elementThis[16-12]; |
1934 | | value2 += pi[iRow1] * elementThis[16-11]; |
1935 | | } |
1936 | | case 5: { |
1937 | | int iRow0 = rowThis[16-10]; |
1938 | | int iRow1 = rowThis[16-9]; |
1939 | | value += pi[iRow0] * elementThis[16-10]; |
1940 | | value2 += pi[iRow1] * elementThis[16-9]; |
1941 | | } |
1942 | | case 4: { |
1943 | | int iRow0 = rowThis[16-8]; |
1944 | | int iRow1 = rowThis[16-7]; |
1945 | | value += pi[iRow0] * elementThis[16-8]; |
1946 | | value2 += pi[iRow1] * elementThis[16-7]; |
1947 | | } |
1948 | | case 3: { |
1949 | | int iRow0 = rowThis[16-6]; |
1950 | | int iRow1 = rowThis[16-5]; |
1951 | | value += pi[iRow0] * elementThis[16-6]; |
1952 | | value2 += pi[iRow1] * elementThis[16-5]; |
1953 | | } |
1954 | | case 2: { |
1955 | | int iRow0 = rowThis[16-4]; |
1956 | | int iRow1 = rowThis[16-3]; |
1957 | | value += pi[iRow0] * elementThis[16-4]; |
1958 | | value2 += pi[iRow1] * elementThis[16-3]; |
1959 | | } |
1960 | | case 1: { |
1961 | | int iRow0 = rowThis[16-2]; |
1962 | | int iRow1 = rowThis[16-1]; |
1963 | | value += pi[iRow0] * elementThis[16-2]; |
1964 | | value2 += pi[iRow1] * elementThis[16-1]; |
1965 | | } |
1966 | | case 0: |
1967 | | ; |
1968 | | } |
1969 | | } |
1970 | | value += value2; |
1971 | | #endif |
1972 | | if (fabs(value) > zeroTolerance) { |
1973 | | double mult = multiplier[wanted-1]; |
1974 | | double alpha = value * mult; |
1975 | | array[numberNonZero] = value; |
1976 | | index[numberNonZero++] = iColumn; |
1977 | | if (alpha > 0.0) { |
1978 | | double oldValue = reducedCost[iColumn] * mult; |
1979 | | double value = oldValue - tentativeTheta * alpha; |
1980 | | if (value < dualT) { |
1981 | | value = oldValue - upperTheta * alpha; |
1982 | | if (value < dualT && alpha >= acceptablePivot) { |
1983 | | upperTheta = (oldValue - dualT) / alpha; |
1984 | | //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta); |
1985 | | } |
1986 | | // add to list |
1987 | | spareArray[numberRemaining] = alpha * mult; |
1988 | | spareIndex[numberRemaining++] = iColumn; |
1989 | | } |
1990 | | } |
1991 | | } |
1992 | | } |
1993 | | } |
| 1881 | const int *COIN_RESTRICT rowThis = &row[end - 16]; |
| 1882 | const double *COIN_RESTRICT elementThis = &elementByColumn[end - 16]; |
| 1883 | bool odd = (n & 1) != 0; |
| 1884 | n = n >> 1; |
| 1885 | double value2 = 0.0; |
| 1886 | if (odd) { |
| 1887 | int iRow = row[start]; |
| 1888 | value2 = pi[iRow] * elementByColumn[start]; |
| 1889 | } |
| 1890 | switch (n) { |
| 1891 | default: { |
| 1892 | if (odd) { |
| 1893 | start++; |
| 1894 | } |
| 1895 | n -= 8; |
| 1896 | for (; n; n--) { |
| 1897 | int iRow0 = row[start]; |
| 1898 | int iRow1 = row[start + 1]; |
| 1899 | value += pi[iRow0] * elementByColumn[start]; |
| 1900 | value2 += pi[iRow1] * elementByColumn[start + 1]; |
| 1901 | start += 2; |
| 1902 | } |
| 1903 | case 8: { |
| 1904 | int iRow0 = rowThis[16 - 16]; |
| 1905 | int iRow1 = rowThis[16 - 15]; |
| 1906 | value += pi[iRow0] * elementThis[16 - 16]; |
| 1907 | value2 += pi[iRow1] * elementThis[16 - 15]; |
| 1908 | } |
| 1909 | case 7: { |
| 1910 | int iRow0 = rowThis[16 - 14]; |
| 1911 | int iRow1 = rowThis[16 - 13]; |
| 1912 | value += pi[iRow0] * elementThis[16 - 14]; |
| 1913 | value2 += pi[iRow1] * elementThis[16 - 13]; |
| 1914 | } |
| 1915 | case 6: { |
| 1916 | int iRow0 = rowThis[16 - 12]; |
| 1917 | int iRow1 = rowThis[16 - 11]; |
| 1918 | value += pi[iRow0] * elementThis[16 - 12]; |
| 1919 | value2 += pi[iRow1] * elementThis[16 - 11]; |
| 1920 | } |
| 1921 | case 5: { |
| 1922 | int iRow0 = rowThis[16 - 10]; |
| 1923 | int iRow1 = rowThis[16 - 9]; |
| 1924 | value += pi[iRow0] * elementThis[16 - 10]; |
| 1925 | value2 += pi[iRow1] * elementThis[16 - 9]; |
| 1926 | } |
| 1927 | case 4: { |
| 1928 | int iRow0 = rowThis[16 - 8]; |
| 1929 | int iRow1 = rowThis[16 - 7]; |
| 1930 | value += pi[iRow0] * elementThis[16 - 8]; |
| 1931 | value2 += pi[iRow1] * elementThis[16 - 7]; |
| 1932 | } |
| 1933 | case 3: { |
| 1934 | int iRow0 = rowThis[16 - 6]; |
| 1935 | int iRow1 = rowThis[16 - 5]; |
| 1936 | value += pi[iRow0] * elementThis[16 - 6]; |
| 1937 | value2 += pi[iRow1] * elementThis[16 - 5]; |
| 1938 | } |
| 1939 | case 2: { |
| 1940 | int iRow0 = rowThis[16 - 4]; |
| 1941 | int iRow1 = rowThis[16 - 3]; |
| 1942 | value += pi[iRow0] * elementThis[16 - 4]; |
| 1943 | value2 += pi[iRow1] * elementThis[16 - 3]; |
| 1944 | } |
| 1945 | case 1: { |
| 1946 | int iRow0 = rowThis[16 - 2]; |
| 1947 | int iRow1 = rowThis[16 - 1]; |
| 1948 | value += pi[iRow0] * elementThis[16 - 2]; |
| 1949 | value2 += pi[iRow1] * elementThis[16 - 1]; |
| 1950 | } |
| 1951 | case 0:; |
| 1952 | } |
| 1953 | } |
| 1954 | value += value2; |
| 1955 | #endif |
| 1956 | if (fabs(value) > zeroTolerance) { |
| 1957 | double mult = multiplier[wanted - 1]; |
| 1958 | double alpha = value * mult; |
| 1959 | array[numberNonZero] = value; |
| 1960 | index[numberNonZero++] = iColumn; |
| 1961 | if (alpha > 0.0) { |
| 1962 | double oldValue = reducedCost[iColumn] * mult; |
| 1963 | double value = oldValue - tentativeTheta * alpha; |
| 1964 | if (value < dualT) { |
| 1965 | value = oldValue - upperTheta * alpha; |
| 1966 | if (value < dualT && alpha >= acceptablePivot) { |
| 1967 | upperTheta = (oldValue - dualT) / alpha; |
| 1968 | //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta); |
| 1969 | } |
| 1970 | // add to list |
| 1971 | spareArray[numberRemaining] = alpha * mult; |
| 1972 | spareIndex[numberRemaining++] = iColumn; |
| 1973 | } |
| 1974 | } |
| 1975 | } |
| 1976 | } |
| 1977 | } |
2002 | | int |
2003 | | ClpPackedMatrix::gutsOfTransposeTimesScaled(const double * COIN_RESTRICT pi, |
2004 | | const double * COIN_RESTRICT columnScale, |
2005 | | int * COIN_RESTRICT index, |
2006 | | double * COIN_RESTRICT array, |
2007 | | const unsigned char * COIN_RESTRICT status, const double zeroTolerance) const |
2008 | | { |
2009 | | int numberNonZero = 0; |
2010 | | // get matrix data pointers |
2011 | | const int * COIN_RESTRICT row = matrix_->getIndices(); |
2012 | | const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
2013 | | const double * COIN_RESTRICT elementByColumn = matrix_->getElements(); |
2014 | | double value = 0.0; |
2015 | | int jColumn = -1; |
2016 | | for (int iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
2017 | | bool wanted = ((status[iColumn] & 3) != 1); |
2018 | | if (fabs(value) > zeroTolerance) { |
2019 | | array[numberNonZero] = value; |
2020 | | index[numberNonZero++] = jColumn; |
2021 | | } |
2022 | | value = 0.0; |
2023 | | if (wanted) { |
2024 | | double scale = columnScale[iColumn]; |
2025 | | CoinBigIndex start = columnStart[iColumn]; |
2026 | | CoinBigIndex end = columnStart[iColumn+1]; |
2027 | | jColumn = iColumn; |
2028 | | for (CoinBigIndex j = start; j < end; j++) { |
2029 | | int iRow = row[j]; |
2030 | | value += pi[iRow] * elementByColumn[j]; |
2031 | | } |
2032 | | value *= scale; |
2033 | | } |
2034 | | } |
2035 | | if (fabs(value) > zeroTolerance) { |
2036 | | array[numberNonZero] = value; |
2037 | | index[numberNonZero++] = jColumn; |
2038 | | } |
2039 | | return numberNonZero; |
| 1986 | int ClpPackedMatrix::gutsOfTransposeTimesScaled(const double *COIN_RESTRICT pi, |
| 1987 | const double *COIN_RESTRICT columnScale, |
| 1988 | int *COIN_RESTRICT index, |
| 1989 | double *COIN_RESTRICT array, |
| 1990 | const unsigned char *COIN_RESTRICT status, const double zeroTolerance) const |
| 1991 | { |
| 1992 | int numberNonZero = 0; |
| 1993 | // get matrix data pointers |
| 1994 | const int *COIN_RESTRICT row = matrix_->getIndices(); |
| 1995 | const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
| 1996 | const double *COIN_RESTRICT elementByColumn = matrix_->getElements(); |
| 1997 | double value = 0.0; |
| 1998 | int jColumn = -1; |
| 1999 | for (int iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 2000 | bool wanted = ((status[iColumn] & 3) != 1); |
| 2001 | if (fabs(value) > zeroTolerance) { |
| 2002 | array[numberNonZero] = value; |
| 2003 | index[numberNonZero++] = jColumn; |
| 2004 | } |
| 2005 | value = 0.0; |
| 2006 | if (wanted) { |
| 2007 | double scale = columnScale[iColumn]; |
| 2008 | CoinBigIndex start = columnStart[iColumn]; |
| 2009 | CoinBigIndex end = columnStart[iColumn + 1]; |
| 2010 | jColumn = iColumn; |
| 2011 | for (CoinBigIndex j = start; j < end; j++) { |
| 2012 | int iRow = row[j]; |
| 2013 | value += pi[iRow] * elementByColumn[j]; |
| 2014 | } |
| 2015 | value *= scale; |
| 2016 | } |
| 2017 | } |
| 2018 | if (fabs(value) > zeroTolerance) { |
| 2019 | array[numberNonZero] = value; |
| 2020 | index[numberNonZero++] = jColumn; |
| 2021 | } |
| 2022 | return numberNonZero; |
2065 | | int |
2066 | | ClpPackedMatrix::gutsOfTransposeTimesByRowGEK(const CoinIndexedVector * COIN_RESTRICT piVector, |
2067 | | int * COIN_RESTRICT index, |
2068 | | double * COIN_RESTRICT output, |
2069 | | int numberColumns, |
2070 | | const double tolerance, |
2071 | | const double scalar) const |
2072 | | { |
2073 | | const double * COIN_RESTRICT pi = piVector->denseVector(); |
2074 | | int numberInRowArray = piVector->getNumElements(); |
2075 | | const int * COIN_RESTRICT column = matrix_->getIndices(); |
2076 | | const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts(); |
2077 | | const double * COIN_RESTRICT element = matrix_->getElements(); |
2078 | | const int * COIN_RESTRICT whichRow = piVector->getIndices(); |
2079 | | // ** Row copy is already scaled |
2080 | | for (int i = 0; i < numberInRowArray; i++) { |
2081 | | int iRow = whichRow[i]; |
2082 | | double value = pi[i] * scalar; |
2083 | | CoinBigIndex start = rowStart[iRow]; |
2084 | | CoinBigIndex end = rowStart[iRow+1]; |
2085 | | int n = static_cast<int>(end - start); |
2086 | | const int * COIN_RESTRICT columnThis = column + start; |
2087 | | const double * COIN_RESTRICT elementThis = element + start; |
| 2048 | int ClpPackedMatrix::gutsOfTransposeTimesByRowGEK(const CoinIndexedVector *COIN_RESTRICT piVector, |
| 2049 | int *COIN_RESTRICT index, |
| 2050 | double *COIN_RESTRICT output, |
| 2051 | int numberColumns, |
| 2052 | const double tolerance, |
| 2053 | const double scalar) const |
| 2054 | { |
| 2055 | const double *COIN_RESTRICT pi = piVector->denseVector(); |
| 2056 | int numberInRowArray = piVector->getNumElements(); |
| 2057 | const int *COIN_RESTRICT column = matrix_->getIndices(); |
| 2058 | const CoinBigIndex *COIN_RESTRICT rowStart = matrix_->getVectorStarts(); |
| 2059 | const double *COIN_RESTRICT element = matrix_->getElements(); |
| 2060 | const int *COIN_RESTRICT whichRow = piVector->getIndices(); |
| 2061 | // ** Row copy is already scaled |
| 2062 | for (int i = 0; i < numberInRowArray; i++) { |
| 2063 | int iRow = whichRow[i]; |
| 2064 | double value = pi[i] * scalar; |
| 2065 | CoinBigIndex start = rowStart[iRow]; |
| 2066 | CoinBigIndex end = rowStart[iRow + 1]; |
| 2067 | int n = static_cast< int >(end - start); |
| 2068 | const int *COIN_RESTRICT columnThis = column + start; |
| 2069 | const double *COIN_RESTRICT elementThis = element + start; |
2102 | | int numberThreads=abcState(); |
2103 | | if (numberThreads) { |
2104 | | clpTempInfo info[ABOCA_LITE]; |
2105 | | int chunk = (numberColumns+numberThreads-1)/numberThreads; |
2106 | | int n=0; |
2107 | | for (int i=0;i<numberThreads;i++) { |
2108 | | info[i].which = index+n; |
2109 | | info[i].infeas = output+n; |
2110 | | info[i].startColumn=n; |
2111 | | info[i].numberToDo=CoinMin(chunk,numberColumns-n); |
2112 | | info[i].tolerance=tolerance; |
2113 | | n += chunk; |
2114 | | } |
2115 | | for (int i=0;i<numberThreads;i++) { |
2116 | | cilk_spawn packDownBit(info[i]); |
2117 | | } |
2118 | | cilk_sync; |
2119 | | for (int i=0;i<numberThreads;i++) { |
2120 | | numberNonZero += info[i].numberAdded; |
2121 | | } |
2122 | | moveAndZero(info,2,NULL); |
2123 | | } else { |
2124 | | #endif |
2125 | | for (int i = 0; i < numberColumns; i++) { |
2126 | | double value = output[i]; |
2127 | | if (value) { |
2128 | | output[i] = 0.0; |
2129 | | if (fabs(value) > tolerance) { |
2130 | | output[numberNonZero] = value; |
2131 | | index[numberNonZero++] = i; |
2132 | | } |
2133 | | } |
2134 | | } |
| 2084 | int numberThreads = abcState(); |
| 2085 | if (numberThreads) { |
| 2086 | clpTempInfo info[ABOCA_LITE]; |
| 2087 | int chunk = (numberColumns + numberThreads - 1) / numberThreads; |
| 2088 | int n = 0; |
| 2089 | for (int i = 0; i < numberThreads; i++) { |
| 2090 | info[i].which = index + n; |
| 2091 | info[i].infeas = output + n; |
| 2092 | info[i].startColumn = n; |
| 2093 | info[i].numberToDo = CoinMin(chunk, numberColumns - n); |
| 2094 | info[i].tolerance = tolerance; |
| 2095 | n += chunk; |
| 2096 | } |
| 2097 | for (int i = 0; i < numberThreads; i++) { |
| 2098 | cilk_spawn packDownBit(info[i]); |
| 2099 | } |
| 2100 | cilk_sync; |
| 2101 | for (int i = 0; i < numberThreads; i++) { |
| 2102 | numberNonZero += info[i].numberAdded; |
| 2103 | } |
| 2104 | moveAndZero(info, 2, NULL); |
| 2105 | } else { |
| 2106 | #endif |
| 2107 | for (int i = 0; i < numberColumns; i++) { |
| 2108 | double value = output[i]; |
| 2109 | if (value) { |
| 2110 | output[i] = 0.0; |
| 2111 | if (fabs(value) > tolerance) { |
| 2112 | output[numberNonZero] = value; |
| 2113 | index[numberNonZero++] = i; |
| 2114 | } |
| 2115 | } |
| 2116 | } |
2145 | | void |
2146 | | ClpPackedMatrix::gutsOfTransposeTimesByRowEQ2(const CoinIndexedVector * piVector, CoinIndexedVector * output, |
2147 | | CoinIndexedVector * spareVector, const double tolerance, const double scalar) const |
2148 | | { |
2149 | | double * COIN_RESTRICT pi = piVector->denseVector(); |
2150 | | int numberNonZero = 0; |
2151 | | int * COIN_RESTRICT index = output->getIndices(); |
2152 | | double * COIN_RESTRICT array = output->denseVector(); |
2153 | | const int * COIN_RESTRICT column = matrix_->getIndices(); |
2154 | | const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts(); |
2155 | | const double * COIN_RESTRICT element = matrix_->getElements(); |
2156 | | const int * COIN_RESTRICT whichRow = piVector->getIndices(); |
2157 | | int iRow0 = whichRow[0]; |
2158 | | int iRow1 = whichRow[1]; |
2159 | | double pi0 = pi[0]; |
2160 | | double pi1 = pi[1]; |
2161 | | if (rowStart[iRow0+1] - rowStart[iRow0] > |
2162 | | rowStart[iRow1+1] - rowStart[iRow1]) { |
2163 | | // do one with fewer first |
2164 | | iRow0 = iRow1; |
2165 | | iRow1 = whichRow[0]; |
2166 | | pi0 = pi1; |
2167 | | pi1 = pi[0]; |
2168 | | } |
2169 | | // and set up mark as char array |
2170 | | char * COIN_RESTRICT marked = reinterpret_cast<char *> (index + output->capacity()); |
2171 | | int * COIN_RESTRICT lookup = spareVector->getIndices(); |
2172 | | double value = pi0 * scalar; |
2173 | | CoinBigIndex j; |
2174 | | for (j = rowStart[iRow0]; j < rowStart[iRow0+1]; j++) { |
2175 | | int iColumn = column[j]; |
2176 | | double elValue = element[j]; |
2177 | | double value2 = value * elValue; |
2178 | | array[numberNonZero] = value2; |
2179 | | marked[iColumn] = 1; |
2180 | | lookup[iColumn] = numberNonZero; |
2181 | | index[numberNonZero++] = iColumn; |
2182 | | } |
2183 | | value = pi1 * scalar; |
2184 | | for (j = rowStart[iRow1]; j < rowStart[iRow1+1]; j++) { |
2185 | | int iColumn = column[j]; |
2186 | | double elValue = element[j]; |
2187 | | double value2 = value * elValue; |
2188 | | // I am assuming no zeros in matrix |
2189 | | if (marked[iColumn]) { |
2190 | | int iLookup = lookup[iColumn]; |
2191 | | array[iLookup] += value2; |
2192 | | } else { |
2193 | | if (fabs(value2) > tolerance) { |
2194 | | array[numberNonZero] = value2; |
2195 | | index[numberNonZero++] = iColumn; |
2196 | | } |
2197 | | } |
2198 | | } |
2199 | | // get rid of tiny values and zero out marked |
2200 | | int i; |
2201 | | int n = numberNonZero; |
2202 | | numberNonZero=0; |
2203 | | for (i = 0; i < n; i++) { |
2204 | | int iColumn = index[i]; |
2205 | | marked[iColumn] = 0; |
2206 | | if (fabs(array[i]) > tolerance) { |
2207 | | array[numberNonZero] = array[i]; |
2208 | | index[numberNonZero++] = index[i]; |
2209 | | } |
2210 | | } |
2211 | | memset(array+numberNonZero,0,(n-numberNonZero)*sizeof(double)); |
2212 | | output->setNumElements(numberNonZero); |
2213 | | spareVector->setNumElements(0); |
| 2127 | void ClpPackedMatrix::gutsOfTransposeTimesByRowEQ2(const CoinIndexedVector *piVector, CoinIndexedVector *output, |
| 2128 | CoinIndexedVector *spareVector, const double tolerance, const double scalar) const |
| 2129 | { |
| 2130 | double *COIN_RESTRICT pi = piVector->denseVector(); |
| 2131 | int numberNonZero = 0; |
| 2132 | int *COIN_RESTRICT index = output->getIndices(); |
| 2133 | double *COIN_RESTRICT array = output->denseVector(); |
| 2134 | const int *COIN_RESTRICT column = matrix_->getIndices(); |
| 2135 | const CoinBigIndex *COIN_RESTRICT rowStart = matrix_->getVectorStarts(); |
| 2136 | const double *COIN_RESTRICT element = matrix_->getElements(); |
| 2137 | const int *COIN_RESTRICT whichRow = piVector->getIndices(); |
| 2138 | int iRow0 = whichRow[0]; |
| 2139 | int iRow1 = whichRow[1]; |
| 2140 | double pi0 = pi[0]; |
| 2141 | double pi1 = pi[1]; |
| 2142 | if (rowStart[iRow0 + 1] - rowStart[iRow0] > rowStart[iRow1 + 1] - rowStart[iRow1]) { |
| 2143 | // do one with fewer first |
| 2144 | iRow0 = iRow1; |
| 2145 | iRow1 = whichRow[0]; |
| 2146 | pi0 = pi1; |
| 2147 | pi1 = pi[0]; |
| 2148 | } |
| 2149 | // and set up mark as char array |
| 2150 | char *COIN_RESTRICT marked = reinterpret_cast< char * >(index + output->capacity()); |
| 2151 | int *COIN_RESTRICT lookup = spareVector->getIndices(); |
| 2152 | double value = pi0 * scalar; |
| 2153 | CoinBigIndex j; |
| 2154 | for (j = rowStart[iRow0]; j < rowStart[iRow0 + 1]; j++) { |
| 2155 | int iColumn = column[j]; |
| 2156 | double elValue = element[j]; |
| 2157 | double value2 = value * elValue; |
| 2158 | array[numberNonZero] = value2; |
| 2159 | marked[iColumn] = 1; |
| 2160 | lookup[iColumn] = numberNonZero; |
| 2161 | index[numberNonZero++] = iColumn; |
| 2162 | } |
| 2163 | value = pi1 * scalar; |
| 2164 | for (j = rowStart[iRow1]; j < rowStart[iRow1 + 1]; j++) { |
| 2165 | int iColumn = column[j]; |
| 2166 | double elValue = element[j]; |
| 2167 | double value2 = value * elValue; |
| 2168 | // I am assuming no zeros in matrix |
| 2169 | if (marked[iColumn]) { |
| 2170 | int iLookup = lookup[iColumn]; |
| 2171 | array[iLookup] += value2; |
| 2172 | } else { |
| 2173 | if (fabs(value2) > tolerance) { |
| 2174 | array[numberNonZero] = value2; |
| 2175 | index[numberNonZero++] = iColumn; |
| 2176 | } |
| 2177 | } |
| 2178 | } |
| 2179 | // get rid of tiny values and zero out marked |
| 2180 | int i; |
| 2181 | int n = numberNonZero; |
| 2182 | numberNonZero = 0; |
| 2183 | for (i = 0; i < n; i++) { |
| 2184 | int iColumn = index[i]; |
| 2185 | marked[iColumn] = 0; |
| 2186 | if (fabs(array[i]) > tolerance) { |
| 2187 | array[numberNonZero] = array[i]; |
| 2188 | index[numberNonZero++] = index[i]; |
| 2189 | } |
| 2190 | } |
| 2191 | memset(array + numberNonZero, 0, (n - numberNonZero) * sizeof(double)); |
| 2192 | output->setNumElements(numberNonZero); |
| 2193 | spareVector->setNumElements(0); |
2216 | | void |
2217 | | ClpPackedMatrix::gutsOfTransposeTimesByRowEQ1(const CoinIndexedVector * piVector, CoinIndexedVector * output, |
2218 | | const double tolerance, const double scalar) const |
2219 | | { |
2220 | | double * COIN_RESTRICT pi = piVector->denseVector(); |
2221 | | int numberNonZero = 0; |
2222 | | int * COIN_RESTRICT index = output->getIndices(); |
2223 | | double * COIN_RESTRICT array = output->denseVector(); |
2224 | | const int * COIN_RESTRICT column = matrix_->getIndices(); |
2225 | | const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts(); |
2226 | | const double * COIN_RESTRICT element = matrix_->getElements(); |
2227 | | int iRow = piVector->getIndices()[0]; |
2228 | | numberNonZero = 0; |
2229 | | CoinBigIndex j; |
2230 | | double value = pi[0] * scalar; |
2231 | | for (j = rowStart[iRow]; j < rowStart[iRow+1]; j++) { |
2232 | | int iColumn = column[j]; |
2233 | | double elValue = element[j]; |
2234 | | double value2 = value * elValue; |
2235 | | if (fabs(value2) > tolerance) { |
2236 | | array[numberNonZero] = value2; |
2237 | | index[numberNonZero++] = iColumn; |
2238 | | } |
2239 | | } |
2240 | | output->setNumElements(numberNonZero); |
| 2196 | void ClpPackedMatrix::gutsOfTransposeTimesByRowEQ1(const CoinIndexedVector *piVector, CoinIndexedVector *output, |
| 2197 | const double tolerance, const double scalar) const |
| 2198 | { |
| 2199 | double *COIN_RESTRICT pi = piVector->denseVector(); |
| 2200 | int numberNonZero = 0; |
| 2201 | int *COIN_RESTRICT index = output->getIndices(); |
| 2202 | double *COIN_RESTRICT array = output->denseVector(); |
| 2203 | const int *COIN_RESTRICT column = matrix_->getIndices(); |
| 2204 | const CoinBigIndex *COIN_RESTRICT rowStart = matrix_->getVectorStarts(); |
| 2205 | const double *COIN_RESTRICT element = matrix_->getElements(); |
| 2206 | int iRow = piVector->getIndices()[0]; |
| 2207 | numberNonZero = 0; |
| 2208 | CoinBigIndex j; |
| 2209 | double value = pi[0] * scalar; |
| 2210 | for (j = rowStart[iRow]; j < rowStart[iRow + 1]; j++) { |
| 2211 | int iColumn = column[j]; |
| 2212 | double elValue = element[j]; |
| 2213 | double value2 = value * elValue; |
| 2214 | if (fabs(value2) > tolerance) { |
| 2215 | array[numberNonZero] = value2; |
| 2216 | index[numberNonZero++] = iColumn; |
| 2217 | } |
| 2218 | } |
| 2219 | output->setNumElements(numberNonZero); |
2245 | | void |
2246 | | ClpPackedMatrix::subsetTransposeTimes(const ClpSimplex * model, |
2247 | | const CoinIndexedVector * rowArray, |
2248 | | const CoinIndexedVector * y, |
2249 | | CoinIndexedVector * columnArray) const |
2250 | | { |
2251 | | columnArray->clear(); |
2252 | | double * COIN_RESTRICT pi = rowArray->denseVector(); |
2253 | | double * COIN_RESTRICT array = columnArray->denseVector(); |
2254 | | int jColumn; |
2255 | | // get matrix data pointers |
2256 | | const int * COIN_RESTRICT row = matrix_->getIndices(); |
2257 | | const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
2258 | | const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths(); |
2259 | | const double * COIN_RESTRICT elementByColumn = matrix_->getElements(); |
2260 | | const double * COIN_RESTRICT rowScale = model->rowScale(); |
2261 | | int numberToDo = y->getNumElements(); |
2262 | | const int * COIN_RESTRICT which = y->getIndices(); |
2263 | | assert (!rowArray->packedMode()); |
2264 | | columnArray->setPacked(); |
2265 | | ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix(); |
2266 | | int flags = flags_; |
2267 | | if (rowScale && scaledMatrix && !(scaledMatrix->flags() & 2)) { |
2268 | | flags = 0; |
2269 | | rowScale = NULL; |
2270 | | // get matrix data pointers |
2271 | | row = scaledMatrix->getIndices(); |
2272 | | columnStart = scaledMatrix->getVectorStarts(); |
2273 | | elementByColumn = scaledMatrix->getElements(); |
2274 | | } |
2275 | | if (!(flags & 2) && numberToDo>2) { |
2276 | | // no gaps |
2277 | | if (!rowScale) { |
2278 | | int iColumn = which[0]; |
2279 | | double value = 0.0; |
2280 | | CoinBigIndex j; |
2281 | | int columnNext = which[1]; |
2282 | | CoinBigIndex startNext=columnStart[columnNext]; |
2283 | | //coin_prefetch_const(row+startNext); |
2284 | | //coin_prefetch_const(elementByColumn+startNext); |
2285 | | CoinBigIndex endNext=columnStart[columnNext+1]; |
2286 | | for (j = columnStart[iColumn]; |
2287 | | j < columnStart[iColumn+1]; j++) { |
2288 | | int iRow = row[j]; |
2289 | | value += pi[iRow] * elementByColumn[j]; |
2290 | | } |
2291 | | for (jColumn = 0; jColumn < numberToDo - 2; jColumn++) { |
2292 | | CoinBigIndex start = startNext; |
2293 | | CoinBigIndex end = endNext; |
2294 | | columnNext = which[jColumn+2]; |
2295 | | startNext=columnStart[columnNext]; |
2296 | | //coin_prefetch_const(row+startNext); |
2297 | | //coin_prefetch_const(elementByColumn+startNext); |
2298 | | endNext=columnStart[columnNext+1]; |
2299 | | array[jColumn] = value; |
2300 | | value = 0.0; |
2301 | | for (j = start; j < end; j++) { |
2302 | | int iRow = row[j]; |
2303 | | value += pi[iRow] * elementByColumn[j]; |
2304 | | } |
2305 | | } |
2306 | | array[jColumn++] = value; |
2307 | | value = 0.0; |
2308 | | for (j = startNext; j < endNext; j++) { |
2309 | | int iRow = row[j]; |
2310 | | value += pi[iRow] * elementByColumn[j]; |
2311 | | } |
2312 | | array[jColumn] = value; |
2313 | | } else { |
| 2224 | void ClpPackedMatrix::subsetTransposeTimes(const ClpSimplex *model, |
| 2225 | const CoinIndexedVector *rowArray, |
| 2226 | const CoinIndexedVector *y, |
| 2227 | CoinIndexedVector *columnArray) const |
| 2228 | { |
| 2229 | columnArray->clear(); |
| 2230 | double *COIN_RESTRICT pi = rowArray->denseVector(); |
| 2231 | double *COIN_RESTRICT array = columnArray->denseVector(); |
| 2232 | int jColumn; |
| 2233 | // get matrix data pointers |
| 2234 | const int *COIN_RESTRICT row = matrix_->getIndices(); |
| 2235 | const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
| 2236 | const int *COIN_RESTRICT columnLength = matrix_->getVectorLengths(); |
| 2237 | const double *COIN_RESTRICT elementByColumn = matrix_->getElements(); |
| 2238 | const double *COIN_RESTRICT rowScale = model->rowScale(); |
| 2239 | int numberToDo = y->getNumElements(); |
| 2240 | const int *COIN_RESTRICT which = y->getIndices(); |
| 2241 | assert(!rowArray->packedMode()); |
| 2242 | columnArray->setPacked(); |
| 2243 | ClpPackedMatrix *scaledMatrix = model->clpScaledMatrix(); |
| 2244 | int flags = flags_; |
| 2245 | if (rowScale && scaledMatrix && !(scaledMatrix->flags() & 2)) { |
| 2246 | flags = 0; |
| 2247 | rowScale = NULL; |
| 2248 | // get matrix data pointers |
| 2249 | row = scaledMatrix->getIndices(); |
| 2250 | columnStart = scaledMatrix->getVectorStarts(); |
| 2251 | elementByColumn = scaledMatrix->getElements(); |
| 2252 | } |
| 2253 | if (!(flags & 2) && numberToDo > 2) { |
| 2254 | // no gaps |
| 2255 | if (!rowScale) { |
| 2256 | int iColumn = which[0]; |
| 2257 | double value = 0.0; |
| 2258 | CoinBigIndex j; |
| 2259 | int columnNext = which[1]; |
| 2260 | CoinBigIndex startNext = columnStart[columnNext]; |
| 2261 | //coin_prefetch_const(row+startNext); |
| 2262 | //coin_prefetch_const(elementByColumn+startNext); |
| 2263 | CoinBigIndex endNext = columnStart[columnNext + 1]; |
| 2264 | for (j = columnStart[iColumn]; |
| 2265 | j < columnStart[iColumn + 1]; j++) { |
| 2266 | int iRow = row[j]; |
| 2267 | value += pi[iRow] * elementByColumn[j]; |
| 2268 | } |
| 2269 | for (jColumn = 0; jColumn < numberToDo - 2; jColumn++) { |
| 2270 | CoinBigIndex start = startNext; |
| 2271 | CoinBigIndex end = endNext; |
| 2272 | columnNext = which[jColumn + 2]; |
| 2273 | startNext = columnStart[columnNext]; |
| 2274 | //coin_prefetch_const(row+startNext); |
| 2275 | //coin_prefetch_const(elementByColumn+startNext); |
| 2276 | endNext = columnStart[columnNext + 1]; |
| 2277 | array[jColumn] = value; |
| 2278 | value = 0.0; |
| 2279 | for (j = start; j < end; j++) { |
| 2280 | int iRow = row[j]; |
| 2281 | value += pi[iRow] * elementByColumn[j]; |
| 2282 | } |
| 2283 | } |
| 2284 | array[jColumn++] = value; |
| 2285 | value = 0.0; |
| 2286 | for (j = startNext; j < endNext; j++) { |
| 2287 | int iRow = row[j]; |
| 2288 | value += pi[iRow] * elementByColumn[j]; |
| 2289 | } |
| 2290 | array[jColumn] = value; |
| 2291 | } else { |
2315 | | if (model->clpScaledMatrix()) |
2316 | | printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__); |
2317 | | #endif |
2318 | | // scaled |
2319 | | const double * columnScale = model->columnScale(); |
2320 | | int iColumn = which[0]; |
2321 | | double value = 0.0; |
2322 | | double scale = columnScale[iColumn]; |
2323 | | CoinBigIndex j; |
2324 | | for (j = columnStart[iColumn]; |
2325 | | j < columnStart[iColumn+1]; j++) { |
2326 | | int iRow = row[j]; |
2327 | | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
2328 | | } |
2329 | | for (jColumn = 0; jColumn < numberToDo - 1; jColumn++) { |
2330 | | int iColumn = which[jColumn+1]; |
2331 | | value *= scale; |
2332 | | scale = columnScale[iColumn]; |
2333 | | CoinBigIndex start = columnStart[iColumn]; |
2334 | | CoinBigIndex end = columnStart[iColumn+1]; |
2335 | | array[jColumn] = value; |
2336 | | value = 0.0; |
2337 | | for (j = start; j < end; j++) { |
2338 | | int iRow = row[j]; |
2339 | | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
2340 | | } |
2341 | | } |
2342 | | value *= scale; |
2343 | | array[jColumn] = value; |
2344 | | } |
2345 | | } else if (numberToDo) { |
2346 | | // gaps |
2347 | | if (!rowScale) { |
2348 | | for (jColumn = 0; jColumn < numberToDo; jColumn++) { |
2349 | | int iColumn = which[jColumn]; |
2350 | | double value = 0.0; |
2351 | | CoinBigIndex j; |
2352 | | for (j = columnStart[iColumn]; |
2353 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
2354 | | int iRow = row[j]; |
2355 | | value += pi[iRow] * elementByColumn[j]; |
2356 | | } |
2357 | | array[jColumn] = value; |
2358 | | } |
2359 | | } else { |
| 2293 | if (model->clpScaledMatrix()) |
| 2294 | printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__); |
| 2295 | #endif |
| 2296 | // scaled |
| 2297 | const double *columnScale = model->columnScale(); |
| 2298 | int iColumn = which[0]; |
| 2299 | double value = 0.0; |
| 2300 | double scale = columnScale[iColumn]; |
| 2301 | CoinBigIndex j; |
| 2302 | for (j = columnStart[iColumn]; |
| 2303 | j < columnStart[iColumn + 1]; j++) { |
| 2304 | int iRow = row[j]; |
| 2305 | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
| 2306 | } |
| 2307 | for (jColumn = 0; jColumn < numberToDo - 1; jColumn++) { |
| 2308 | int iColumn = which[jColumn + 1]; |
| 2309 | value *= scale; |
| 2310 | scale = columnScale[iColumn]; |
| 2311 | CoinBigIndex start = columnStart[iColumn]; |
| 2312 | CoinBigIndex end = columnStart[iColumn + 1]; |
| 2313 | array[jColumn] = value; |
| 2314 | value = 0.0; |
| 2315 | for (j = start; j < end; j++) { |
| 2316 | int iRow = row[j]; |
| 2317 | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
| 2318 | } |
| 2319 | } |
| 2320 | value *= scale; |
| 2321 | array[jColumn] = value; |
| 2322 | } |
| 2323 | } else if (numberToDo) { |
| 2324 | // gaps |
| 2325 | if (!rowScale) { |
| 2326 | for (jColumn = 0; jColumn < numberToDo; jColumn++) { |
| 2327 | int iColumn = which[jColumn]; |
| 2328 | double value = 0.0; |
| 2329 | CoinBigIndex j; |
| 2330 | for (j = columnStart[iColumn]; |
| 2331 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 2332 | int iRow = row[j]; |
| 2333 | value += pi[iRow] * elementByColumn[j]; |
| 2334 | } |
| 2335 | array[jColumn] = value; |
| 2336 | } |
| 2337 | } else { |
2384 | | bool |
2385 | | ClpPackedMatrix::canCombine(const ClpSimplex * model, |
2386 | | const CoinIndexedVector * pi) const |
2387 | | { |
2388 | | int numberInRowArray = pi->getNumElements(); |
2389 | | int numberRows = model->numberRows(); |
2390 | | bool packed = pi->packedMode(); |
2391 | | // factor should be smaller if doing both with two pi vectors |
2392 | | double factor = 0.30; |
2393 | | // We may not want to do by row if there may be cache problems |
2394 | | // It would be nice to find L2 cache size - for moment 512K |
2395 | | // Be slightly optimistic |
2396 | | if (numberActiveColumns_ * sizeof(double) > 1000000) { |
2397 | | if (numberRows * 10 < numberActiveColumns_) |
2398 | | factor *= 0.333333333; |
2399 | | else if (numberRows * 4 < numberActiveColumns_) |
2400 | | factor *= 0.5; |
2401 | | else if (numberRows * 2 < numberActiveColumns_) |
2402 | | factor *= 0.66666666667; |
2403 | | //if (model->numberIterations()%50==0) |
2404 | | //printf("%d nonzero\n",numberInRowArray); |
2405 | | } |
2406 | | // if not packed then bias a bit more towards by column |
2407 | | if (!packed) |
2408 | | factor *= 0.9; |
2409 | | // bias if columnCopy |
2410 | | if (columnCopy_) |
2411 | | factor *= 0.5; |
2412 | | return ((numberInRowArray > factor * numberRows || !model->rowCopy()) && !(flags_ & 2)); |
| 2362 | bool ClpPackedMatrix::canCombine(const ClpSimplex *model, |
| 2363 | const CoinIndexedVector *pi) const |
| 2364 | { |
| 2365 | int numberInRowArray = pi->getNumElements(); |
| 2366 | int numberRows = model->numberRows(); |
| 2367 | bool packed = pi->packedMode(); |
| 2368 | // factor should be smaller if doing both with two pi vectors |
| 2369 | double factor = 0.30; |
| 2370 | // We may not want to do by row if there may be cache problems |
| 2371 | // It would be nice to find L2 cache size - for moment 512K |
| 2372 | // Be slightly optimistic |
| 2373 | if (numberActiveColumns_ * sizeof(double) > 1000000) { |
| 2374 | if (numberRows * 10 < numberActiveColumns_) |
| 2375 | factor *= 0.333333333; |
| 2376 | else if (numberRows * 4 < numberActiveColumns_) |
| 2377 | factor *= 0.5; |
| 2378 | else if (numberRows * 2 < numberActiveColumns_) |
| 2379 | factor *= 0.66666666667; |
| 2380 | //if (model->numberIterations()%50==0) |
| 2381 | //printf("%d nonzero\n",numberInRowArray); |
| 2382 | } |
| 2383 | // if not packed then bias a bit more towards by column |
| 2384 | if (!packed) |
| 2385 | factor *= 0.9; |
| 2386 | // bias if columnCopy |
| 2387 | if (columnCopy_) |
| 2388 | factor *= 0.5; |
| 2389 | return ((numberInRowArray > factor * numberRows || !model->rowCopy()) && !(flags_ & 2)); |
2420 | | transposeTimes2UnscaledBit(clpTempInfo & info) |
2421 | | { |
2422 | | const CoinBigIndex * COIN_RESTRICT columnStart=info.start; |
2423 | | const int * COIN_RESTRICT row = info.row; |
2424 | | const double * COIN_RESTRICT elementByColumn = info.element; |
2425 | | double * COIN_RESTRICT array = info.infeas; |
2426 | | double * COIN_RESTRICT weights = const_cast<double *>(info.lower); |
2427 | | double * COIN_RESTRICT reducedCost = info.reducedCost; |
2428 | | double * COIN_RESTRICT infeas = const_cast<double *>(info.upper); |
2429 | | int * COIN_RESTRICT index = info.which; |
2430 | | unsigned int * COIN_RESTRICT reference = reinterpret_cast<unsigned int *>(info.index); |
2431 | | const unsigned char * COIN_RESTRICT status = info.status; |
2432 | | double zeroTolerance=info.tolerance; |
2433 | | double dualTolerance=info.dualTolerance; |
2434 | | const double * COIN_RESTRICT pi = info.work; |
2435 | | const double * COIN_RESTRICT piWeight = info.cost; |
| 2397 | transposeTimes2UnscaledBit(clpTempInfo &info) |
| 2398 | { |
| 2399 | const CoinBigIndex *COIN_RESTRICT columnStart = info.start; |
| 2400 | const int *COIN_RESTRICT row = info.row; |
| 2401 | const double *COIN_RESTRICT elementByColumn = info.element; |
| 2402 | double *COIN_RESTRICT array = info.infeas; |
| 2403 | double *COIN_RESTRICT weights = const_cast< double * >(info.lower); |
| 2404 | double *COIN_RESTRICT reducedCost = info.reducedCost; |
| 2405 | double *COIN_RESTRICT infeas = const_cast< double * >(info.upper); |
| 2406 | int *COIN_RESTRICT index = info.which; |
| 2407 | unsigned int *COIN_RESTRICT reference = reinterpret_cast< unsigned int * >(info.index); |
| 2408 | const unsigned char *COIN_RESTRICT status = info.status; |
| 2409 | double zeroTolerance = info.tolerance; |
| 2410 | double dualTolerance = info.dualTolerance; |
| 2411 | const double *COIN_RESTRICT pi = info.work; |
| 2412 | const double *COIN_RESTRICT piWeight = info.cost; |
2458 | | // and do other array |
2459 | | double modification = 0.0; |
2460 | | for (int i=0; i < n; i ++) { |
2461 | | int iRow = rowThis[i]; |
2462 | | modification += piWeight[iRow] * elementThis[i]; |
2463 | | } |
2464 | | double thisWeight = weights[iColumn]; |
2465 | | double oldWeight=thisWeight; |
2466 | | double pivot = value * scaleFactor; |
2467 | | double pivotSquared = pivot * pivot; |
2468 | | thisWeight += pivotSquared * devex + pivot * modification; |
2469 | | debug3(iColumn,thisWeight,pivotSquared,devex,pivot,modification,oldWeight); |
2470 | | if (thisWeight < DEVEX_TRY_NORM) { |
2471 | | if (referenceIn < 0.0) { |
2472 | | // steepest |
2473 | | thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); |
2474 | | } else { |
2475 | | // exact |
2476 | | thisWeight = referenceIn * pivotSquared; |
2477 | | if (reference(iColumn)) |
2478 | | thisWeight += 1.0; |
2479 | | thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); |
2480 | | } |
2481 | | } |
2482 | | weights[iColumn] = thisWeight; |
2483 | | if (!killDjs) { |
2484 | | value = reducedCost[iColumn]-value; |
2485 | | reducedCost[iColumn] = value; |
2486 | | if (iStat>1&&iStat<4) { |
2487 | | value *= mmult[iStat-2]; |
2488 | | if (value<=dualTolerance) |
2489 | | value=0.0; |
2490 | | } else { |
2491 | | // free or super basic |
2492 | | if (fabs(value) > FREE_ACCEPT * dualTolerance) { |
2493 | | // we are going to bias towards free (but only if reasonable) |
2494 | | value *= FREE_BIAS; |
2495 | | } else { |
2496 | | value=0.0; |
2497 | | } |
2498 | | } |
2499 | | if (value) { |
2500 | | value *= value; |
2501 | | // store square in list |
2502 | | if (infeas[iColumn]) { |
2503 | | infeas[iColumn] = value; // already there |
2504 | | } else { |
2505 | | array[numberNonZero]=value; |
2506 | | index[numberNonZero++]=iColumn; |
2507 | | } |
2508 | | } else { |
2509 | | array[numberNonZero]=0.0; |
2510 | | index[numberNonZero++]=iColumn; |
2511 | | } |
2512 | | } |
2513 | | } |
2514 | | } |
2515 | | } |
2516 | | info.numberAdded=numberNonZero; |
| 2435 | // and do other array |
| 2436 | double modification = 0.0; |
| 2437 | for (int i = 0; i < n; i++) { |
| 2438 | int iRow = rowThis[i]; |
| 2439 | modification += piWeight[iRow] * elementThis[i]; |
| 2440 | } |
| 2441 | double thisWeight = weights[iColumn]; |
| 2442 | double oldWeight = thisWeight; |
| 2443 | double pivot = value * scaleFactor; |
| 2444 | double pivotSquared = pivot * pivot; |
| 2445 | thisWeight += pivotSquared * devex + pivot * modification; |
| 2446 | debug3(iColumn, thisWeight, pivotSquared, devex, pivot, modification, oldWeight); |
| 2447 | if (thisWeight < DEVEX_TRY_NORM) { |
| 2448 | if (referenceIn < 0.0) { |
| 2449 | // steepest |
| 2450 | thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); |
| 2451 | } else { |
| 2452 | // exact |
| 2453 | thisWeight = referenceIn * pivotSquared; |
| 2454 | if (reference(iColumn)) |
| 2455 | thisWeight += 1.0; |
| 2456 | thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); |
| 2457 | } |
| 2458 | } |
| 2459 | weights[iColumn] = thisWeight; |
| 2460 | if (!killDjs) { |
| 2461 | value = reducedCost[iColumn] - value; |
| 2462 | reducedCost[iColumn] = value; |
| 2463 | if (iStat > 1 && iStat < 4) { |
| 2464 | value *= mmult[iStat - 2]; |
| 2465 | if (value <= dualTolerance) |
| 2466 | value = 0.0; |
| 2467 | } else { |
| 2468 | // free or super basic |
| 2469 | if (fabs(value) > FREE_ACCEPT * dualTolerance) { |
| 2470 | // we are going to bias towards free (but only if reasonable) |
| 2471 | value *= FREE_BIAS; |
| 2472 | } else { |
| 2473 | value = 0.0; |
| 2474 | } |
| 2475 | } |
| 2476 | if (value) { |
| 2477 | value *= value; |
| 2478 | // store square in list |
| 2479 | if (infeas[iColumn]) { |
| 2480 | infeas[iColumn] = value; // already there |
| 2481 | } else { |
| 2482 | array[numberNonZero] = value; |
| 2483 | index[numberNonZero++] = iColumn; |
| 2484 | } |
| 2485 | } else { |
| 2486 | array[numberNonZero] = 0.0; |
| 2487 | index[numberNonZero++] = iColumn; |
| 2488 | } |
| 2489 | } |
| 2490 | } |
| 2491 | } |
| 2492 | } |
| 2493 | info.numberAdded = numberNonZero; |
2519 | | transposeTimes2ScaledBit(clpTempInfo & info) |
2520 | | { |
2521 | | const CoinBigIndex * COIN_RESTRICT columnStart=info.start; |
2522 | | const int * COIN_RESTRICT row = info.row; |
2523 | | const double * COIN_RESTRICT elementByColumn = info.element; |
2524 | | double * COIN_RESTRICT array = info.infeas; |
2525 | | double * COIN_RESTRICT weights = const_cast<double *>(info.lower); |
2526 | | int * COIN_RESTRICT index = info.which; |
2527 | | unsigned int * COIN_RESTRICT reference = reinterpret_cast<unsigned int *>(info.index); |
2528 | | const unsigned char * COIN_RESTRICT status = info.status; |
2529 | | double zeroTolerance=info.tolerance; |
2530 | | double dualTolerance=info.dualTolerance; |
2531 | | double * COIN_RESTRICT reducedCost = info.reducedCost; |
2532 | | double * COIN_RESTRICT infeas = const_cast<double *>(info.upper); |
2533 | | const double * COIN_RESTRICT pi = info.work; |
2534 | | const double * COIN_RESTRICT piWeight = info.cost; |
2535 | | const double * COIN_RESTRICT columnScale = info.solution; |
| 2496 | transposeTimes2ScaledBit(clpTempInfo &info) |
| 2497 | { |
| 2498 | const CoinBigIndex *COIN_RESTRICT columnStart = info.start; |
| 2499 | const int *COIN_RESTRICT row = info.row; |
| 2500 | const double *COIN_RESTRICT elementByColumn = info.element; |
| 2501 | double *COIN_RESTRICT array = info.infeas; |
| 2502 | double *COIN_RESTRICT weights = const_cast< double * >(info.lower); |
| 2503 | int *COIN_RESTRICT index = info.which; |
| 2504 | unsigned int *COIN_RESTRICT reference = reinterpret_cast< unsigned int * >(info.index); |
| 2505 | const unsigned char *COIN_RESTRICT status = info.status; |
| 2506 | double zeroTolerance = info.tolerance; |
| 2507 | double dualTolerance = info.dualTolerance; |
| 2508 | double *COIN_RESTRICT reducedCost = info.reducedCost; |
| 2509 | double *COIN_RESTRICT infeas = const_cast< double * >(info.upper); |
| 2510 | const double *COIN_RESTRICT pi = info.work; |
| 2511 | const double *COIN_RESTRICT piWeight = info.cost; |
| 2512 | const double *COIN_RESTRICT columnScale = info.solution; |
2568 | | // and do other array |
2569 | | double modification = 0.0; |
2570 | | for (int i=0; i < n; i += 2) { |
2571 | | int iRow0 = rowThis[i]; |
2572 | | int iRow1 = rowThis[i+1]; |
2573 | | modification += piWeight[iRow0] * elementThis[i]; |
2574 | | modification += piWeight[iRow1] * elementThis[i+1]; |
2575 | | } |
2576 | | if (odd) { |
2577 | | int iRow = rowThis[n]; |
2578 | | modification += piWeight[iRow] * elementThis[n]; |
2579 | | } |
2580 | | modification *= scale; |
2581 | | double thisWeight = weights[iColumn]; |
2582 | | double pivot = value * scaleFactor; |
2583 | | double pivotSquared = pivot * pivot; |
2584 | | thisWeight += pivotSquared * devex + pivot * modification; |
2585 | | if (thisWeight < DEVEX_TRY_NORM) { |
2586 | | if (referenceIn < 0.0) { |
2587 | | // steepest |
2588 | | thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); |
2589 | | } else { |
2590 | | // exact |
2591 | | thisWeight = referenceIn * pivotSquared; |
2592 | | if (reference(iColumn)) |
2593 | | thisWeight += 1.0; |
2594 | | thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); |
2595 | | } |
2596 | | } |
2597 | | weights[iColumn] = thisWeight; |
2598 | | if (!killDjs) { |
2599 | | value = reducedCost[iColumn]-value; |
2600 | | reducedCost[iColumn] = value; |
2601 | | if (iStat>1&&iStat<4) { |
2602 | | value *= mmult[iStat-2]; |
2603 | | if (value<=dualTolerance) |
2604 | | value=0.0; |
2605 | | } else { |
2606 | | // free or super basic |
2607 | | if (fabs(value) > FREE_ACCEPT * dualTolerance) { |
2608 | | // we are going to bias towards free (but only if reasonable) |
2609 | | value *= FREE_BIAS; |
2610 | | } else { |
2611 | | value=0.0; |
2612 | | } |
2613 | | } |
2614 | | if (value) { |
2615 | | value *= value; |
2616 | | // store square in list |
2617 | | if (infeas[iColumn]) { |
2618 | | infeas[iColumn] = value; // already there |
2619 | | } else { |
2620 | | array[numberNonZero]=value; |
2621 | | index[numberNonZero++]=iColumn; |
2622 | | } |
2623 | | } else { |
2624 | | array[numberNonZero]=0.0; |
2625 | | index[numberNonZero++]=iColumn; |
2626 | | } |
2627 | | } |
2628 | | } |
2629 | | } |
2630 | | } |
2631 | | info.numberAdded=numberNonZero; |
| 2545 | // and do other array |
| 2546 | double modification = 0.0; |
| 2547 | for (int i = 0; i < n; i += 2) { |
| 2548 | int iRow0 = rowThis[i]; |
| 2549 | int iRow1 = rowThis[i + 1]; |
| 2550 | modification += piWeight[iRow0] * elementThis[i]; |
| 2551 | modification += piWeight[iRow1] * elementThis[i + 1]; |
| 2552 | } |
| 2553 | if (odd) { |
| 2554 | int iRow = rowThis[n]; |
| 2555 | modification += piWeight[iRow] * elementThis[n]; |
| 2556 | } |
| 2557 | modification *= scale; |
| 2558 | double thisWeight = weights[iColumn]; |
| 2559 | double pivot = value * scaleFactor; |
| 2560 | double pivotSquared = pivot * pivot; |
| 2561 | thisWeight += pivotSquared * devex + pivot * modification; |
| 2562 | if (thisWeight < DEVEX_TRY_NORM) { |
| 2563 | if (referenceIn < 0.0) { |
| 2564 | // steepest |
| 2565 | thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); |
| 2566 | } else { |
| 2567 | // exact |
| 2568 | thisWeight = referenceIn * pivotSquared; |
| 2569 | if (reference(iColumn)) |
| 2570 | thisWeight += 1.0; |
| 2571 | thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); |
| 2572 | } |
| 2573 | } |
| 2574 | weights[iColumn] = thisWeight; |
| 2575 | if (!killDjs) { |
| 2576 | value = reducedCost[iColumn] - value; |
| 2577 | reducedCost[iColumn] = value; |
| 2578 | if (iStat > 1 && iStat < 4) { |
| 2579 | value *= mmult[iStat - 2]; |
| 2580 | if (value <= dualTolerance) |
| 2581 | value = 0.0; |
| 2582 | } else { |
| 2583 | // free or super basic |
| 2584 | if (fabs(value) > FREE_ACCEPT * dualTolerance) { |
| 2585 | // we are going to bias towards free (but only if reasonable) |
| 2586 | value *= FREE_BIAS; |
| 2587 | } else { |
| 2588 | value = 0.0; |
| 2589 | } |
| 2590 | } |
| 2591 | if (value) { |
| 2592 | value *= value; |
| 2593 | // store square in list |
| 2594 | if (infeas[iColumn]) { |
| 2595 | infeas[iColumn] = value; // already there |
| 2596 | } else { |
| 2597 | array[numberNonZero] = value; |
| 2598 | index[numberNonZero++] = iColumn; |
| 2599 | } |
| 2600 | } else { |
| 2601 | array[numberNonZero] = 0.0; |
| 2602 | index[numberNonZero++] = iColumn; |
| 2603 | } |
| 2604 | } |
| 2605 | } |
| 2606 | } |
| 2607 | } |
| 2608 | info.numberAdded = numberNonZero; |
2638 | | int |
2639 | | ClpPackedMatrix::transposeTimes2(const ClpSimplex * model, |
2640 | | const CoinIndexedVector * pi1, CoinIndexedVector * dj1, |
2641 | | const CoinIndexedVector * pi2, |
2642 | | CoinIndexedVector * spare, |
2643 | | double * COIN_RESTRICT infeas, |
2644 | | double * COIN_RESTRICT reducedCost, |
2645 | | double referenceIn, double devex, |
2646 | | // Array for exact devex to say what is in reference framework |
2647 | | unsigned int * COIN_RESTRICT reference, |
2648 | | double * COIN_RESTRICT weights, double scaleFactor) |
2649 | | { |
2650 | | int returnCode=0; |
2651 | | // put row of tableau in dj1 |
2652 | | double * COIN_RESTRICT pi = pi1->denseVector(); |
2653 | | int numberNonZero = 0; |
2654 | | int * COIN_RESTRICT index = dj1->getIndices(); |
2655 | | double * COIN_RESTRICT array = dj1->denseVector(); |
2656 | | int numberInRowArray = pi1->getNumElements(); |
2657 | | double zeroTolerance = model->zeroTolerance(); |
2658 | | double dualTolerance = model->currentDualTolerance(); |
2659 | | // we can't really trust infeasibilities if there is dual error |
2660 | | // this coding has to mimic coding in checkDualSolution |
2661 | | double error = CoinMin(1.0e-2, model->largestDualError()); |
2662 | | // allow tolerance at least slightly bigger than standard |
2663 | | dualTolerance = dualTolerance + error; |
2664 | | bool packed = pi1->packedMode(); |
2665 | | // do by column |
2666 | | int iColumn; |
2667 | | // get matrix data pointers |
2668 | | const int * COIN_RESTRICT row = matrix_->getIndices(); |
2669 | | const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
2670 | | const double * COIN_RESTRICT elementByColumn = matrix_->getElements(); |
2671 | | const double * COIN_RESTRICT rowScale = model->rowScale(); |
2672 | | assert (!spare->getNumElements()); |
2673 | | assert (numberActiveColumns_ > 0); |
2674 | | double * COIN_RESTRICT piWeight = pi2->denseVector(); |
2675 | | assert (!pi2->packedMode()); |
2676 | | bool killDjs = (scaleFactor == 0.0); |
2677 | | if (!scaleFactor) |
2678 | | scaleFactor = 1.0; |
2679 | | if (packed) { |
2680 | | // need to expand pi into y |
2681 | | assert(spare->capacity() >= model->numberRows()); |
2682 | | double * COIN_RESTRICT piOld = pi; |
2683 | | pi = spare->denseVector(); |
2684 | | const int * COIN_RESTRICT whichRow = pi1->getIndices(); |
2685 | | int i; |
2686 | | ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix(); |
2687 | | if (rowScale && scaledMatrix) { |
2688 | | rowScale = NULL; |
2689 | | // get matrix data pointers |
2690 | | row = scaledMatrix->getIndices(); |
2691 | | columnStart = scaledMatrix->getVectorStarts(); |
2692 | | elementByColumn = scaledMatrix->getElements(); |
2693 | | } |
2694 | | if (!rowScale) { |
2695 | | // modify pi so can collapse to one loop |
2696 | | for (i = 0; i < numberInRowArray; i++) { |
2697 | | int iRow = whichRow[i]; |
2698 | | pi[iRow] = piOld[i]; |
2699 | | } |
2700 | | if (!columnCopy_ || killDjs) { |
2701 | | if (infeas) |
2702 | | returnCode=1; |
| 2615 | int ClpPackedMatrix::transposeTimes2(const ClpSimplex *model, |
| 2616 | const CoinIndexedVector *pi1, CoinIndexedVector *dj1, |
| 2617 | const CoinIndexedVector *pi2, |
| 2618 | CoinIndexedVector *spare, |
| 2619 | double *COIN_RESTRICT infeas, |
| 2620 | double *COIN_RESTRICT reducedCost, |
| 2621 | double referenceIn, double devex, |
| 2622 | // Array for exact devex to say what is in reference framework |
| 2623 | unsigned int *COIN_RESTRICT reference, |
| 2624 | double *COIN_RESTRICT weights, double scaleFactor) |
| 2625 | { |
| 2626 | int returnCode = 0; |
| 2627 | // put row of tableau in dj1 |
| 2628 | double *COIN_RESTRICT pi = pi1->denseVector(); |
| 2629 | int numberNonZero = 0; |
| 2630 | int *COIN_RESTRICT index = dj1->getIndices(); |
| 2631 | double *COIN_RESTRICT array = dj1->denseVector(); |
| 2632 | int numberInRowArray = pi1->getNumElements(); |
| 2633 | double zeroTolerance = model->zeroTolerance(); |
| 2634 | double dualTolerance = model->currentDualTolerance(); |
| 2635 | // we can't really trust infeasibilities if there is dual error |
| 2636 | // this coding has to mimic coding in checkDualSolution |
| 2637 | double error = CoinMin(1.0e-2, model->largestDualError()); |
| 2638 | // allow tolerance at least slightly bigger than standard |
| 2639 | dualTolerance = dualTolerance + error; |
| 2640 | bool packed = pi1->packedMode(); |
| 2641 | // do by column |
| 2642 | int iColumn; |
| 2643 | // get matrix data pointers |
| 2644 | const int *COIN_RESTRICT row = matrix_->getIndices(); |
| 2645 | const CoinBigIndex *COIN_RESTRICT columnStart = matrix_->getVectorStarts(); |
| 2646 | const double *COIN_RESTRICT elementByColumn = matrix_->getElements(); |
| 2647 | const double *COIN_RESTRICT rowScale = model->rowScale(); |
| 2648 | assert(!spare->getNumElements()); |
| 2649 | assert(numberActiveColumns_ > 0); |
| 2650 | double *COIN_RESTRICT piWeight = pi2->denseVector(); |
| 2651 | assert(!pi2->packedMode()); |
| 2652 | bool killDjs = (scaleFactor == 0.0); |
| 2653 | if (!scaleFactor) |
| 2654 | scaleFactor = 1.0; |
| 2655 | if (packed) { |
| 2656 | // need to expand pi into y |
| 2657 | assert(spare->capacity() >= model->numberRows()); |
| 2658 | double *COIN_RESTRICT piOld = pi; |
| 2659 | pi = spare->denseVector(); |
| 2660 | const int *COIN_RESTRICT whichRow = pi1->getIndices(); |
| 2661 | int i; |
| 2662 | ClpPackedMatrix *scaledMatrix = model->clpScaledMatrix(); |
| 2663 | if (rowScale && scaledMatrix) { |
| 2664 | rowScale = NULL; |
| 2665 | // get matrix data pointers |
| 2666 | row = scaledMatrix->getIndices(); |
| 2667 | columnStart = scaledMatrix->getVectorStarts(); |
| 2668 | elementByColumn = scaledMatrix->getElements(); |
| 2669 | } |
| 2670 | if (!rowScale) { |
| 2671 | // modify pi so can collapse to one loop |
| 2672 | for (i = 0; i < numberInRowArray; i++) { |
| 2673 | int iRow = whichRow[i]; |
| 2674 | pi[iRow] = piOld[i]; |
| 2675 | } |
| 2676 | if (!columnCopy_ || killDjs) { |
| 2677 | if (infeas) |
| 2678 | returnCode = 1; |
2704 | | int numberThreads=abcState(); |
2705 | | if (numberThreads) { |
2706 | | clpTempInfo info[ABOCA_LITE]; |
2707 | | int chunk = (numberActiveColumns_+numberThreads-1)/numberThreads; |
2708 | | int n=0; |
2709 | | for (int i=0;i<numberThreads;i++) { |
2710 | | info[i].theta = referenceIn; |
2711 | | info[i].changeObj = devex; |
2712 | | info[i].primalRatio = scaleFactor; |
2713 | | info[i].lower = weights; |
2714 | | info[i].reducedCost = reducedCost; |
2715 | | info[i].upper = infeas; |
2716 | | info[i].which = index+n; |
2717 | | info[i].infeas = array+n; |
2718 | | info[i].index = reinterpret_cast<int *>(reference); |
2719 | | info[i].work=const_cast<double *>(pi); |
2720 | | info[i].cost=const_cast<double *>(piWeight); |
2721 | | info[i].status=const_cast<unsigned char *>(model->statusArray()); |
2722 | | info[i].startColumn=n; |
2723 | | info[i].element=elementByColumn; |
2724 | | info[i].start=columnStart; |
2725 | | info[i].row=row; |
2726 | | info[i].numberToDo=CoinMin(chunk,numberActiveColumns_-n); |
2727 | | info[i].tolerance=zeroTolerance; |
2728 | | info[i].dualTolerance=dualTolerance; |
2729 | | info[i].numberInfeasibilities=killDjs ? 1 : 0; |
2730 | | n += chunk; |
2731 | | } |
2732 | | for (int i=0;i<numberThreads;i++) { |
2733 | | cilk_spawn transposeTimes2UnscaledBit(info[i]); |
2734 | | } |
2735 | | cilk_sync; |
2736 | | for (int i=0;i<numberThreads;i++) { |
2737 | | numberNonZero += info[i].numberAdded; |
2738 | | } |
2739 | | moveAndZero(info,2,NULL); |
2740 | | } else { |
2741 | | #endif |
2742 | | CoinBigIndex j; |
2743 | | CoinBigIndex end = columnStart[0]; |
2744 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
2745 | | CoinBigIndex start = end; |
2746 | | end = columnStart[iColumn+1]; |
2747 | | ClpSimplex::Status status = model->getStatus(iColumn); |
2748 | | if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue; |
2749 | | double value = 0.0; |
2750 | | for (j = start; j < end; j++) { |
2751 | | int iRow = row[j]; |
2752 | | value -= pi[iRow] * elementByColumn[j]; |
2753 | | } |
2754 | | if (fabs(value) > zeroTolerance) { |
2755 | | // and do other array |
2756 | | double modification = 0.0; |
2757 | | for (j = start; j < end; j++) { |
2758 | | int iRow = row[j]; |
2759 | | modification += piWeight[iRow] * elementByColumn[j]; |
2760 | | } |
2761 | | double thisWeight = weights[iColumn]; |
2762 | | double oldWeight=thisWeight; |
2763 | | double pivot = value * scaleFactor; |
2764 | | double pivotSquared = pivot * pivot; |
2765 | | thisWeight += pivotSquared * devex + pivot * modification; |
2766 | | debug3(iColumn,thisWeight,pivotSquared,devex,pivot,modification,oldWeight); |
2767 | | if (thisWeight < DEVEX_TRY_NORM) { |
2768 | | if (referenceIn < 0.0) { |
2769 | | // steepest |
2770 | | thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); |
2771 | | } else { |
2772 | | // exact |
2773 | | thisWeight = referenceIn * pivotSquared; |
2774 | | if (reference(iColumn)) |
2775 | | thisWeight += 1.0; |
2776 | | thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); |
2777 | | } |
2778 | | } |
2779 | | weights[iColumn] = thisWeight; |
2780 | | if (!killDjs) { |
2781 | | value = reducedCost[iColumn]-value; |
2782 | | reducedCost[iColumn] = value; |
2783 | | // simplify status |
2784 | | ClpSimplex::Status status = model->getStatus(iColumn); |
2785 | | |
2786 | | switch(status) { |
2787 | | |
2788 | | case ClpSimplex::basic: |
2789 | | case ClpSimplex::isFixed: |
2790 | | break; |
2791 | | case ClpSimplex::isFree: |
2792 | | case ClpSimplex::superBasic: |
2793 | | if (fabs(value) > FREE_ACCEPT * dualTolerance) { |
2794 | | // we are going to bias towards free (but only if reasonable) |
2795 | | value *= FREE_BIAS; |
2796 | | value *= value; |
2797 | | // store square in list |
2798 | | if (infeas[iColumn]) { |
2799 | | infeas[iColumn] = value; // already there |
2800 | | } else { |
2801 | | array[numberNonZero]=value; |
2802 | | index[numberNonZero++]=iColumn; |
2803 | | } |
2804 | | } else { |
2805 | | array[numberNonZero]=0.0; |
2806 | | index[numberNonZero++]=iColumn; |
2807 | | } |
2808 | | break; |
2809 | | case ClpSimplex::atUpperBound: |
2810 | | if (value > dualTolerance) { |
2811 | | value *= value; |
2812 | | // store square in list |
2813 | | if (infeas[iColumn]) { |
2814 | | infeas[iColumn] = value; // already there |
2815 | | } else { |
2816 | | array[numberNonZero]=value; |
2817 | | index[numberNonZero++]=iColumn; |
2818 | | } |
2819 | | } else { |
2820 | | array[numberNonZero]=0.0; |
2821 | | index[numberNonZero++]=iColumn; |
2822 | | } |
2823 | | break; |
2824 | | case ClpSimplex::atLowerBound: |
2825 | | if (value < -dualTolerance) { |
2826 | | value *= value; |
2827 | | // store square in list |
2828 | | if (infeas[iColumn]) { |
2829 | | infeas[iColumn] = value; // already there |
2830 | | } else { |
2831 | | array[numberNonZero]=value; |
2832 | | index[numberNonZero++]=iColumn; |
2833 | | } |
2834 | | } else { |
2835 | | array[numberNonZero]=0.0; |
2836 | | index[numberNonZero++]=iColumn; |
2837 | | } |
2838 | | } |
2839 | | } |
2840 | | } |
| 2680 | int numberThreads = abcState(); |
| 2681 | if (numberThreads) { |
| 2682 | clpTempInfo info[ABOCA_LITE]; |
| 2683 | int chunk = (numberActiveColumns_ + numberThreads - 1) / numberThreads; |
| 2684 | int n = 0; |
| 2685 | for (int i = 0; i < numberThreads; i++) { |
| 2686 | info[i].theta = referenceIn; |
| 2687 | info[i].changeObj = devex; |
| 2688 | info[i].primalRatio = scaleFactor; |
| 2689 | info[i].lower = weights; |
| 2690 | info[i].reducedCost = reducedCost; |
| 2691 | info[i].upper = infeas; |
| 2692 | info[i].which = index + n; |
| 2693 | info[i].infeas = array + n; |
| 2694 | info[i].index = reinterpret_cast< int * >(reference); |
| 2695 | info[i].work = const_cast< double * >(pi); |
| 2696 | info[i].cost = const_cast< double * >(piWeight); |
| 2697 | info[i].status = const_cast< unsigned char * >(model->statusArray()); |
| 2698 | info[i].startColumn = n; |
| 2699 | info[i].element = elementByColumn; |
| 2700 | info[i].start = columnStart; |
| 2701 | info[i].row = row; |
| 2702 | info[i].numberToDo = CoinMin(chunk, numberActiveColumns_ - n); |
| 2703 | info[i].tolerance = zeroTolerance; |
| 2704 | info[i].dualTolerance = dualTolerance; |
| 2705 | info[i].numberInfeasibilities = killDjs ? 1 : 0; |
| 2706 | n += chunk; |
| 2707 | } |
| 2708 | for (int i = 0; i < numberThreads; i++) { |
| 2709 | cilk_spawn transposeTimes2UnscaledBit(info[i]); |
| 2710 | } |
| 2711 | cilk_sync; |
| 2712 | for (int i = 0; i < numberThreads; i++) { |
| 2713 | numberNonZero += info[i].numberAdded; |
| 2714 | } |
| 2715 | moveAndZero(info, 2, NULL); |
| 2716 | } else { |
| 2717 | #endif |
| 2718 | CoinBigIndex j; |
| 2719 | CoinBigIndex end = columnStart[0]; |
| 2720 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 2721 | CoinBigIndex start = end; |
| 2722 | end = columnStart[iColumn + 1]; |
| 2723 | ClpSimplex::Status status = model->getStatus(iColumn); |
| 2724 | if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) |
| 2725 | continue; |
| 2726 | double value = 0.0; |
| 2727 | for (j = start; j < end; j++) { |
| 2728 | int iRow = row[j]; |
| 2729 | value -= pi[iRow] * elementByColumn[j]; |
| 2730 | } |
| 2731 | if (fabs(value) > zeroTolerance) { |
| 2732 | // and do other array |
| 2733 | double modification = 0.0; |
| 2734 | for (j = start; j < end; j++) { |
| 2735 | int iRow = row[j]; |
| 2736 | modification += piWeight[iRow] * elementByColumn[j]; |
| 2737 | } |
| 2738 | double thisWeight = weights[iColumn]; |
| 2739 | double oldWeight = thisWeight; |
| 2740 | double pivot = value * scaleFactor; |
| 2741 | double pivotSquared = pivot * pivot; |
| 2742 | thisWeight += pivotSquared * devex + pivot * modification; |
| 2743 | debug3(iColumn, thisWeight, pivotSquared, devex, pivot, modification, oldWeight); |
| 2744 | if (thisWeight < DEVEX_TRY_NORM) { |
| 2745 | if (referenceIn < 0.0) { |
| 2746 | // steepest |
| 2747 | thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); |
| 2748 | } else { |
| 2749 | // exact |
| 2750 | thisWeight = referenceIn * pivotSquared; |
| 2751 | if (reference(iColumn)) |
| 2752 | thisWeight += 1.0; |
| 2753 | thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); |
| 2754 | } |
| 2755 | } |
| 2756 | weights[iColumn] = thisWeight; |
| 2757 | if (!killDjs) { |
| 2758 | value = reducedCost[iColumn] - value; |
| 2759 | reducedCost[iColumn] = value; |
| 2760 | // simplify status |
| 2761 | ClpSimplex::Status status = model->getStatus(iColumn); |
| 2762 | |
| 2763 | switch (status) { |
| 2764 | |
| 2765 | case ClpSimplex::basic: |
| 2766 | case ClpSimplex::isFixed: |
| 2767 | break; |
| 2768 | case ClpSimplex::isFree: |
| 2769 | case ClpSimplex::superBasic: |
| 2770 | if (fabs(value) > FREE_ACCEPT * dualTolerance) { |
| 2771 | // we are going to bias towards free (but only if reasonable) |
| 2772 | value *= FREE_BIAS; |
| 2773 | value *= value; |
| 2774 | // store square in list |
| 2775 | if (infeas[iColumn]) { |
| 2776 | infeas[iColumn] = value; // already there |
| 2777 | } else { |
| 2778 | array[numberNonZero] = value; |
| 2779 | index[numberNonZero++] = iColumn; |
2843 | | } |
2844 | | #endif |
2845 | | } else { |
2846 | | // use special column copy |
2847 | | // reset back |
2848 | | assert (!killDjs); |
2849 | | //if (killDjs) |
2850 | | // scaleFactor = 0.0; |
2851 | | if (infeas) |
2852 | | returnCode=1; |
2853 | | columnCopy_->transposeTimes2(model, pi, dj1, piWeight, |
2854 | | infeas, reducedCost, |
2855 | | referenceIn, devex, |
2856 | | reference, weights, scaleFactor); |
2857 | | numberNonZero = dj1->getNumElements(); |
2858 | | } |
2859 | | } else { |
2860 | | // scaled |
2861 | | // modify pi so can collapse to one loop |
2862 | | for (i = 0; i < numberInRowArray; i++) { |
2863 | | int iRow = whichRow[i]; |
2864 | | pi[iRow] = piOld[i] * rowScale[iRow]; |
2865 | | } |
2866 | | // can also scale piWeight as not used again |
2867 | | int numberWeight = pi2->getNumElements(); |
2868 | | const int * indexWeight = pi2->getIndices(); |
2869 | | for (i = 0; i < numberWeight; i++) { |
2870 | | int iRow = indexWeight[i]; |
2871 | | piWeight[iRow] *= rowScale[iRow]; |
2872 | | } |
2873 | | if (!columnCopy_ || killDjs) { |
2874 | | if (infeas) |
2875 | | returnCode=1; |
2876 | | const double * COIN_RESTRICT columnScale = model->columnScale(); |
| 2820 | } |
| 2821 | #endif |
| 2822 | } else { |
| 2823 | // use special column copy |
| 2824 | // reset back |
| 2825 | assert(!killDjs); |
| 2826 | //if (killDjs) |
| 2827 | // scaleFactor = 0.0; |
| 2828 | if (infeas) |
| 2829 | returnCode = 1; |
| 2830 | columnCopy_->transposeTimes2(model, pi, dj1, piWeight, |
| 2831 | infeas, reducedCost, |
| 2832 | referenceIn, devex, |
| 2833 | reference, weights, scaleFactor); |
| 2834 | numberNonZero = dj1->getNumElements(); |
| 2835 | } |
| 2836 | } else { |
| 2837 | // scaled |
| 2838 | // modify pi so can collapse to one loop |
| 2839 | for (i = 0; i < numberInRowArray; i++) { |
| 2840 | int iRow = whichRow[i]; |
| 2841 | pi[iRow] = piOld[i] * rowScale[iRow]; |
| 2842 | } |
| 2843 | // can also scale piWeight as not used again |
| 2844 | int numberWeight = pi2->getNumElements(); |
| 2845 | const int *indexWeight = pi2->getIndices(); |
| 2846 | for (i = 0; i < numberWeight; i++) { |
| 2847 | int iRow = indexWeight[i]; |
| 2848 | piWeight[iRow] *= rowScale[iRow]; |
| 2849 | } |
| 2850 | if (!columnCopy_ || killDjs) { |
| 2851 | if (infeas) |
| 2852 | returnCode = 1; |
| 2853 | const double *COIN_RESTRICT columnScale = model->columnScale(); |
2878 | | int numberThreads=abcState(); |
2879 | | if (numberThreads) { |
2880 | | clpTempInfo info[ABOCA_LITE]; |
2881 | | int chunk = (numberActiveColumns_+numberThreads-1)/numberThreads; |
2882 | | int n=0; |
2883 | | for (int i=0;i<numberThreads;i++) { |
2884 | | info[i].theta = referenceIn; |
2885 | | info[i].changeObj = devex; |
2886 | | info[i].primalRatio = scaleFactor; |
2887 | | info[i].lower = weights; |
2888 | | info[i].reducedCost = reducedCost; |
2889 | | info[i].upper = infeas; |
2890 | | info[i].solution = const_cast<double *>(columnScale); |
2891 | | info[i].which = index+n; |
2892 | | info[i].infeas = array+n; |
2893 | | info[i].index = reinterpret_cast<int *>(reference); |
2894 | | info[i].work=const_cast<double *>(pi); |
2895 | | info[i].cost=const_cast<double *>(piWeight); |
2896 | | info[i].status=const_cast<unsigned char *>(model->statusArray()); |
2897 | | info[i].startColumn=n; |
2898 | | info[i].element=elementByColumn; |
2899 | | info[i].start=columnStart; |
2900 | | info[i].row=row; |
2901 | | info[i].numberToDo=CoinMin(chunk,numberActiveColumns_-n); |
2902 | | info[i].tolerance=zeroTolerance; |
2903 | | info[i].dualTolerance=dualTolerance; |
2904 | | info[i].numberInfeasibilities=killDjs ? 1 : 0; |
2905 | | n += chunk; |
2906 | | } |
2907 | | for (int i=0;i<numberThreads;i++) { |
2908 | | cilk_spawn transposeTimes2ScaledBit(info[i]); |
2909 | | } |
2910 | | cilk_sync; |
2911 | | for (int i=0;i<numberThreads;i++) { |
2912 | | numberNonZero += info[i].numberAdded; |
2913 | | } |
2914 | | moveAndZero(info,2,NULL); |
2915 | | if (infeas) { |
2916 | | returnCode=1; |
2917 | | dj1->setNumElements(numberNonZero); |
2918 | | } |
2919 | | } else { |
2920 | | #endif |
2921 | | CoinBigIndex j; |
2922 | | CoinBigIndex end = columnStart[0]; |
2923 | | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
2924 | | CoinBigIndex start = end; |
2925 | | end = columnStart[iColumn+1]; |
2926 | | ClpSimplex::Status status = model->getStatus(iColumn); |
2927 | | if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue; |
2928 | | double scale = columnScale[iColumn]; |
2929 | | double value = 0.0; |
2930 | | for (j = start; j < end; j++) { |
2931 | | int iRow = row[j]; |
2932 | | value -= pi[iRow] * elementByColumn[j]; |
2933 | | } |
2934 | | value *= scale; |
2935 | | if (fabs(value) > zeroTolerance) { |
2936 | | double modification = 0.0; |
2937 | | for (j = start; j < end; j++) { |
2938 | | int iRow = row[j]; |
2939 | | modification += piWeight[iRow] * elementByColumn[j]; |
2940 | | } |
2941 | | modification *= scale; |
2942 | | double thisWeight = weights[iColumn]; |
2943 | | double pivot = value * scaleFactor; |
2944 | | double pivotSquared = pivot * pivot; |
2945 | | thisWeight += pivotSquared * devex + pivot * modification; |
2946 | | if (thisWeight < DEVEX_TRY_NORM) { |
2947 | | if (referenceIn < 0.0) { |
2948 | | // steepest |
2949 | | thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); |
2950 | | } else { |
2951 | | // exact |
2952 | | thisWeight = referenceIn * pivotSquared; |
2953 | | if (reference(iColumn)) |
2954 | | thisWeight += 1.0; |
2955 | | thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); |
2956 | | } |
2957 | | } |
2958 | | weights[iColumn] = thisWeight; |
2959 | | if (!killDjs) { |
2960 | | value = reducedCost[iColumn]-value; |
2961 | | reducedCost[iColumn] = value; |
2962 | | // simplify status |
2963 | | ClpSimplex::Status status = model->getStatus(iColumn); |
2964 | | |
2965 | | switch(status) { |
2966 | | |
2967 | | case ClpSimplex::basic: |
2968 | | case ClpSimplex::isFixed: |
2969 | | break; |
2970 | | case ClpSimplex::isFree: |
2971 | | case ClpSimplex::superBasic: |
2972 | | if (fabs(value) > FREE_ACCEPT * dualTolerance) { |
2973 | | // we are going to bias towards free (but only if reasonable) |
2974 | | value *= FREE_BIAS; |
2975 | | value *= value; |
2976 | | // store square in list |
2977 | | if (infeas[iColumn]) { |
2978 | | infeas[iColumn] = value; // already there |
2979 | | } else { |
2980 | | array[numberNonZero]=value; |
2981 | | index[numberNonZero++]=iColumn; |
2982 | | } |
2983 | | } else { |
2984 | | array[numberNonZero]=0.0; |
2985 | | index[numberNonZero++]=iColumn; |
2986 | | } |
2987 | | break; |
2988 | | case ClpSimplex::atUpperBound: |
2989 | | if (value > dualTolerance) { |
2990 | | value *= value; |
2991 | | // store square in list |
2992 | | if (infeas[iColumn]) { |
2993 | | infeas[iColumn] = value; // already there |
2994 | | } else { |
2995 | | array[numberNonZero]=value; |
2996 | | index[numberNonZero++]=iColumn; |
2997 | | } |
2998 | | } else { |
2999 | | array[numberNonZero]=0.0; |
3000 | | index[numberNonZero++]=iColumn; |
3001 | | } |
3002 | | break; |
3003 | | case ClpSimplex::atLowerBound: |
3004 | | if (value < -dualTolerance) { |
3005 | | value *= value; |
3006 | | // store square in list |
3007 | | if (infeas[iColumn]) { |
3008 | | infeas[iColumn] = value; // already there |
3009 | | } else { |
3010 | | array[numberNonZero]=value; |
3011 | | index[numberNonZero++]=iColumn; |
3012 | | } |
3013 | | } else { |
3014 | | array[numberNonZero]=0.0; |
3015 | | index[numberNonZero++]=iColumn; |
3016 | | } |
3017 | | } |
3018 | | } |
3019 | | } |
| 2855 | int numberThreads = abcState(); |
| 2856 | if (numberThreads) { |
| 2857 | clpTempInfo info[ABOCA_LITE]; |
| 2858 | int chunk = (numberActiveColumns_ + numberThreads - 1) / numberThreads; |
| 2859 | int n = 0; |
| 2860 | for (int i = 0; i < numberThreads; i++) { |
| 2861 | info[i].theta = referenceIn; |
| 2862 | info[i].changeObj = devex; |
| 2863 | info[i].primalRatio = scaleFactor; |
| 2864 | info[i].lower = weights; |
| 2865 | info[i].reducedCost = reducedCost; |
| 2866 | info[i].upper = infeas; |
| 2867 | info[i].solution = const_cast< double * >(columnScale); |
| 2868 | info[i].which = index + n; |
| 2869 | info[i].infeas = array + n; |
| 2870 | info[i].index = reinterpret_cast< int * >(reference); |
| 2871 | info[i].work = const_cast< double * >(pi); |
| 2872 | info[i].cost = const_cast< double * >(piWeight); |
| 2873 | info[i].status = const_cast< unsigned char * >(model->statusArray()); |
| 2874 | info[i].startColumn = n; |
| 2875 | info[i].element = elementByColumn; |
| 2876 | info[i].start = columnStart; |
| 2877 | info[i].row = row; |
| 2878 | info[i].numberToDo = CoinMin(chunk, numberActiveColumns_ - n); |
| 2879 | info[i].tolerance = zeroTolerance; |
| 2880 | info[i].dualTolerance = dualTolerance; |
| 2881 | info[i].numberInfeasibilities = killDjs ? 1 : 0; |
| 2882 | n += chunk; |
| 2883 | } |
| 2884 | for (int i = 0; i < numberThreads; i++) { |
| 2885 | cilk_spawn transposeTimes2ScaledBit(info[i]); |
| 2886 | } |
| 2887 | cilk_sync; |
| 2888 | for (int i = 0; i < numberThreads; i++) { |
| 2889 | numberNonZero += info[i].numberAdded; |
| 2890 | } |
| 2891 | moveAndZero(info, 2, NULL); |
| 2892 | if (infeas) { |
| 2893 | returnCode = 1; |
| 2894 | dj1->setNumElements(numberNonZero); |
| 2895 | } |
| 2896 | } else { |
| 2897 | #endif |
| 2898 | CoinBigIndex j; |
| 2899 | CoinBigIndex end = columnStart[0]; |
| 2900 | for (iColumn = 0; iColumn < numberActiveColumns_; iColumn++) { |
| 2901 | CoinBigIndex start = end; |
| 2902 | end = columnStart[iColumn + 1]; |
| 2903 | ClpSimplex::Status status = model->getStatus(iColumn); |
| 2904 | if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) |
| 2905 | continue; |
| 2906 | double scale = columnScale[iColumn]; |
| 2907 | double value = 0.0; |
| 2908 | for (j = start; j < end; j++) { |
| 2909 | int iRow = row[j]; |
| 2910 | value -= pi[iRow] * elementByColumn[j]; |
| 2911 | } |
| 2912 | value *= scale; |
| 2913 | if (fabs(value) > zeroTolerance) { |
| 2914 | double modification = 0.0; |
| 2915 | for (j = start; j < end; j++) { |
| 2916 | int iRow = row[j]; |
| 2917 | modification += piWeight[iRow] * elementByColumn[j]; |
| 2918 | } |
| 2919 | modification *= scale; |
| 2920 | double thisWeight = weights[iColumn]; |
| 2921 | double pivot = value * scaleFactor; |
| 2922 | double pivotSquared = pivot * pivot; |
| 2923 | thisWeight += pivotSquared * devex + pivot * modification; |
| 2924 | if (thisWeight < DEVEX_TRY_NORM) { |
| 2925 | if (referenceIn < 0.0) { |
| 2926 | // steepest |
| 2927 | thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); |
| 2928 | } else { |
| 2929 | // exact |
| 2930 | thisWeight = referenceIn * pivotSquared; |
| 2931 | if (reference(iColumn)) |
| 2932 | thisWeight += 1.0; |
| 2933 | thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); |
| 2934 | } |
| 2935 | } |
| 2936 | weights[iColumn] = thisWeight; |
| 2937 | if (!killDjs) { |
| 2938 | value = reducedCost[iColumn] - value; |
| 2939 | reducedCost[iColumn] = value; |