Actual source code: taomat.c

      1: #include "tao_general.h"   /*I "tao_solver.h"  I*/
     2:  #include taomat.h
     3:  #include taovec.h

      8: /*@C
      9:    TaoMatDestroy - Destroys the TaoMat object.

     11:    Input Parameter:
     12: .  TM - the matrix

     14:    Level: beginner
     15: @*/
     16: int TaoMatDestroy( TaoMat* TM){
     17:   TaoFunctionBegin;
     18:   if (TM) delete TM;
     19:   TaoFunctionReturn(0);
     20: }

     24: /*
     25:    SetStructure - Identifies some stucture in the matrix.

     27:    Input Parameter:
     28: .  tmstructure -  TAOMAT_SYMMETRIC_PSD, TAOMAT_SYMMETRIC,
     29: TAOMAT_UNSYMMETRIC;

     31:    Level: intermediate
     32: */
     33: int TaoMat::SetStructure(TaoMatStructure tmstucture){
     34:   TaoFunctionBegin;
     35:   //  tmatstructure=tmstucture;
     36:   TaoFunctionReturn(0);
     37: }

     41: /*@C
     42:    Compatible - Confirms that the operation yy = this * xx is well defined.

     44:    Input Parameters:
     45: .  xx,yy -  vectors like those that will be used for matrix-vector multiplication

     47:    Output Parameters:
     48: .  flag - whether or not a matrix vector product can be performed.

     50:    Level: developer
     51: @*/
     52: int TaoMat::Compatible(TaoVec *xx, TaoVec *yy, TaoTruth *flag){
     53:   TaoFunctionBegin;
     54:   *flag=TAO_FALSE;
     55:   TaoFunctionReturn(0);
     56: }

     61: /*@C
     62:    GetDimensions - Gets the dimensions of the rowspace and columnspace of this matrix.

     64:    Output Parameter:
     65: +  m -  the number of rows
     66: -  n - the number of columns

     68:    Level: intermediate

     70: @*/
     71: int TaoMat::GetDimensions( int* m, int* n ){
     72:   TaoFunctionBegin;
     73:   SETERRQ(56,"Operation not defined");
     74:   /* TaoFunctionReturn(1); */
     75: }

     79: /*@C
     80:    Multiply - Computes  ty = this * tx.

     82:    Input Parameter:
     83: .  tx -  the vector to be multiplied by this matrix.

     85:    Output Parameter:
     86: .  ty -  the destination vector

     88:    Note:
     89:    This method can also be called by using pointers to TaoVec as
     90:    input parameters.

     92:    Level: intermediate

     94: .seealso TaoMat::MultiplyAdd(), TaoMat::MultiplyTranspose()

     96: @*/
     97: int TaoMat::Multiply(TaoVec* tx,TaoVec* ty){
     98:   TaoFunctionBegin;
     99:   SETERRQ(56,"Operation not defined");
    100:   /* TaoFunctionReturn(1); */
    101: }

    106: /*@C
    107:    MultiplyTranspose - Multiplies the transpose of this matrix by a TaoVec.

    109:    Input Parameter:
    110: .  tx -  the vector to be multiplied by this matrix.

    112:    Output Parameter:
    113: .  ty -  the destination vector

    115:    Note:
    116:    This method can also be called by using pointers to TaoVec as
    117:    input parameters.

    119:    Level: intermediate

    121: .seealso TaoMat::Multiply(), TaoMat::MultiplyTransposeAdd()
    122: @*/
    123: int TaoMat::MultiplyTranspose(TaoVec* tx,TaoVec* ty){
    124:   TaoFunctionBegin;
    125:   SETERRQ(56,"Operation not defined");
    126:   /* TaoFunctionReturn(1); */
    127: }

    132: /*@C
    133:    SetDiagonal - Sets the diagonal elements of this matrix with the elements
    134:    of the vector.

    136:    Input Parameter:
    137: .  tv -  the vector containing the diagonal elements

    139:    Note:
    140:    This method can also be called by using a pointer to TaoVec as an
    141:    input parameter.

    143:    Level: intermediate

    145: .seealso TaoMat::AddDiagonal(),TaoMat::ShiftDiagonal()
    146: @*/
    147: int TaoMat::SetDiagonal(TaoVec* tv){
    148:   TaoFunctionBegin;
    149:   SETERRQ(56,"Operation not defined");
    150:   /* TaoFunctionReturn(1); */
    151: }

    155: /*@C
    156:    AddDiagonal - Adds the elements of the vector to the diagonal of this matrix.

    158:    Input Parameter:
    159: .  tv -  the vector containing the diagonal elements

    161:    Note:
    162:    This method can also be called by using a pointer to TaoVec as an
    163:    input parameter.

    165:    Level: intermediate

    167: .seealso TaoMat::SetDiagonal(),TaoMat::ShiftDiagonal()
    168: @*/
    169: int TaoMat::AddDiagonal(TaoVec* tv){
    170:   TaoFunctionBegin;
    171:   SETERRQ(56,"Operation not defined");
    172:   /* TaoFunctionReturn(1); */
    173: }

    177: /*@C
    178:    GetDiagonal - Inserts the diagonal elements of this matrix into the vector.

    180:    Output Parameter:
    181: .  tv -  the vector containing the diagonal elements

    183:    Note:
    184:    This method can also be called by using a pointer to TaoVec as an
    185:    input parameter.

    187:    Level: intermediate

    189: .seealso TaoMat::SetDiagonal()
    190: @*/
    191: int TaoMat::GetDiagonal(TaoVec* tv){
    192:   TaoFunctionBegin;
    193:   SETERRQ(56,"Operation not defined");
    194:   /* TaoFunctionReturn(1); */
    195: }

    199: /*@C
    200:    ShiftDiagonal - Adds this constant to the diagonal elements of this matrix.

    202:    Input Parameter:
    203: .  c -  the constant

    205:    Level: intermediate

    207: .seealso TaoMat::SetDiagonal(), TaoMat::AddDiagonal()
    208: @*/
    209: int TaoMat::ShiftDiagonal(double c){
    210:   TaoFunctionBegin;
    211:   SETERRQ(56,"Operation not defined");
    212:   /* TaoFunctionReturn(1); */
    213: }

    217: /*@C
    218:    Presolve - Prepares to solve a linear system with this matrix by
    219:    doing some initial setup (e.g., computing parts of a preconditioner,
    220:    such as matrix factorization).

    222:    Input:

    224:    Note:
    225:    This routine is optional.  A linear solver object can also be used
    226:    to implement this operation.

    228:    Level: advanced

    230: .seealso TaoMat::Solve()
    231: @*/
    232: int TaoMat::Presolve(){
    233:   TaoFunctionBegin;
    234:   TaoFunctionReturn(0);
    235: }
    236: /*
    239: int TaoMat::minQuadraticTrustRegion(TaoVec*tv,TaoVec*tw,double rr, TaoTruth*tt){
    240:   TaoFunctionBegin;
    241:   SETERRQ(56,"No Linear SOlver has been set.");
    242:   / * TaoFunctionReturn(1); * /
    243: }
    244: */
    247: /*@C
    248:    Solve - Solves the linear system $this tx = tb$.

    250:    Input Parameter:
    251: .  tb -  right hand side

    253:    Output Parameter:
    254: +  tx -  solution vector
    255: -  tt -  TAO_TRUE if solved to prescribed accuracy and TAO_FALSE otherwise

    257:    Level: advanced

    259:    Note:
    260:    This routine is optional.  A linear solver object can also be used
    261:    to implement this operation.

    263: .seealso TaoApplication::GetLinearSolver(),TaoMat::Presolve()
    264: @*/
    265: int TaoMat::Solve(TaoVec* tb, TaoVec *tx, TaoTruth *tt){
    266:   TaoFunctionBegin;
    267:   SETERRQ(56,"No TAO linear solver has been set.");
    268:   /* TaoFunctionReturn(1); */
    269: }

    274: /*@C
    275:   CreateReducedMatrix - Creates a new matrix using the specified rows and columns of this matrix.

    277:    Input Parameter:
    278: +  row -  the rows of this matrix to be extracted
    279: -  column -  the columns of this matrix to be extracted

    281:    Output Parameter:
    282: -  M - the new matrix

    284:    Note:
    285:    This method can also be called by using a pointer to TaoIndexSet as
    286:    input parameters.

    288: .seealso TaoMat::SetReducedMatrix()

    290:    Level: intermediate
    291: @*/
    292: int TaoMat::CreateReducedMatrix(TaoIndexSet* row,TaoIndexSet* col, TaoMat **M){
    293:   TaoFunctionBegin;
    294:   SETERRQ(56,"Operation not defined");
    295:   /* TaoFunctionReturn(1); */
    296: };

    300: /*@C
    301:   SetReducedMatrix - Creates a new matrix using the specified rows and columns of this matrix.

    303:    Input Parameter:
    304: +  row -  the rows of this matrix to be extracted
    305: -  col -  the columns of this matrix to be extracted

    307:    Output Parameter:
    308: .  M - the full new matrix

    310:    Note:
    311:    This method can also be called by using a pointer to TaoIndexSet as
    312:    input parameters.

    314: .seealso TaoMat::CreateReducedMatrix()

    316:    Level: intermediate
    317: @*/
    318: int TaoMat::SetReducedMatrix(TaoMat *M, TaoIndexSet* row,TaoIndexSet* col){
    319:   TaoFunctionBegin;
    320:   SETERRQ(56,"Operation not defined");
    321:   /* TaoFunctionReturn(1); */
    322: };

    327: /*@C
    328:   RowScale - Scales the rows of this matrix.

    330:    Input Parameter:
    331: .  scale - the scaling parameters

    333:    Note:
    334:    This method can also be called by using a pointer to TaoVec as
    335:    input parameters.

    337:    Level: intermediate
    338: @*/
    339: int TaoMat::RowScale(TaoVec* scale){
    340:   TaoFunctionBegin;
    341:   SETERRQ(56,"Operation not defined");
    342:   /* TaoFunctionReturn(1); */
    343: };

    347: /*@C
    348:   ColScale - Scales the columns of this matrix.

    350:    Input Parameter:
    351: .  scale - the scaling parameters

    353:    Note:
    354:    This method can also be called by using a pointer to TaoVec as
    355:    input parameters.

    357:    Level: intermediate
    358: @*/
    359: int TaoMat::ColScale(TaoVec* scale){
    360:   TaoFunctionBegin;
    361:   SETERRQ(56,"Operation not defined");
    362:   /* TaoFunctionReturn(1); */
    363: };

    367: inline static double fischer(double a, double b)
    368: {

    370: #ifdef TODD

    372:   if (TaoAbsScalar(a) > TaoAbsScalar(b)) {
    373:     return sqrt(a*a + b*b) - a - b;
    374:   }
    375:   return sqrt(a*a + b*b) - b - a;

    377: #else

    379:    // Method suggested by Bob Vanderbei

    381:    if (a + b <= 0) {
    382:      return sqrt(a*a + b*b) - (a + b);
    383:    }
    384:    return -2.0*a*b / (sqrt(a*a + b*b) + (a + b));

    386: #endif

    388: }

    392: inline static double norm(double a, double b)
    393: {
    394:   return sqrt(a*a + b*b);
    395: }

    399: /*@C
    400:    D_Fischer - Calculates an element of the B-subdifferential of the
    401:    Fischer-Burmeister function for complementarity problems.
    402:
    403:    Input Parameters:
    404: +  this - the jacobian of tf at tx
    405: .  tx - current point
    406: .  tf - function evaluated at tx
    407: .  tl - lower bounds
    408: .  tu - upper bounds
    409: .  tt1 - work vector
    410: -  tt2 - work vector

    412:    Output Parameters:
    413: +  tda - diagonal perturbation component of the result
    414: -  tdb - row scaling component of the result

    416:    Level: intermediate

    418: .seealso TaoVec::Fischer()
    419: @*/
    420: int TaoMat::D_Fischer(TaoVec *tx, TaoVec *tf, TaoVec *tl, TaoVec *tu,
    421:                       TaoVec *tt1, TaoVec *tt2, TaoVec *tda, TaoVec *tdb )
    422: {
    423:   int i,nn1,nn2,nn3,nn4,nn5,nn6,nn7,nn8,info;
    424:   TaoTruth flag;
    425:   TaoScalar *x,*f,*l,*u,*da,*db,*t1,*t2;
    426:   TaoScalar ai,bi,ci,di,ei;

    428:   TaoFunctionBegin;

    430:   info = this->Compatible(tx, tx, &flag); CHKERRQ(info);
    431:   if (flag==TAO_FALSE) {TaoFunctionReturn(1);}

    433:   info = tx->GetArray(&x,&nn1);CHKERRQ(info);
    434:   info = tf->GetArray(&f,&nn2);CHKERRQ(info);
    435:   info = tl->GetArray(&l,&nn3);CHKERRQ(info);
    436:   info = tu->GetArray(&u,&nn4);CHKERRQ(info);
    437:   info = tda->GetArray(&da,&nn5);CHKERRQ(info);
    438:   info = tdb->GetArray(&db,&nn6);CHKERRQ(info);
    439:   info = tt1->GetArray(&t1,&nn7);CHKERRQ(info);
    440:   info = tt2->GetArray(&t2,&nn8);CHKERRQ(info);

    442:   if (nn1!=nn2 || nn2!=nn3 || nn3!=nn4 || nn4!=nn5 ||
    443:       nn5!=nn6 || nn6!=nn7 || nn7!=nn8) {
    444:     TaoFunctionReturn(1);
    445:   }

    447:   for (i = 0; i < nn1; i++) {
    448:     da[i] = 0;
    449:     db[i] = 0;
    450:     t1[i] = 0;

    452:     if (TaoAbsScalar(f[i]) <= TAO_EPSILON) {
    453:       if (l[i] > -TAO_INFINITY && TaoAbsScalar(x[i] - l[i]) <= TAO_EPSILON) {
    454:         t1[i] = 1;
    455:         da[i] = 1;
    456:       }

    458:       if (u[i] <  TAO_INFINITY && TaoAbsScalar(u[i] - x[i]) <= TAO_EPSILON) {
    459:         t1[i] = 1;
    460:         db[i] = 1;
    461:       }
    462:     }
    463:   }

    465:   info = tt1->RestoreArray(&t1,&nn7);CHKERRQ(info);
    466:   info = tt2->RestoreArray(&t2,&nn8);CHKERRQ(info);
    467:   info = this->Multiply(tt1, tt2); CHKERRQ(info);
    468:   info = tt2->GetArray(&t2,&nn8); CHKERRQ(info);

    470:   for (i = 0; i < nn1; i++) {
    471:     if ((l[i] <= -TAO_INFINITY) && (u[i] >= TAO_INFINITY)) {
    472:       da[i] = 0;
    473:       db[i] = -1;
    474:     }
    475:     else if (l[i] <= -TAO_INFINITY) {
    476:       if (db[i] >= 1) {
    477:         ai = norm(1, t2[i]);

    479:         da[i] = -1/ai - 1;
    480:         db[i] = -t2[i]/ai - 1;
    481:       }
    482:       else {
    483:         bi = u[i] - x[i];
    484:         ai = norm(bi, f[i]);
    485:         ai = TaoMax(TAO_EPSILON, ai);

    487:         da[i] = bi / ai - 1;
    488:         db[i] = -f[i] / ai - 1;
    489:       }
    490:     }
    491:     else if (u[i] >=  TAO_INFINITY) {
    492:       if (da[i] >= 1) {
    493:         ai = norm(1, t2[i]);

    495:         da[i] = 1 / ai - 1;
    496:         db[i] = t2[i] / ai - 1;
    497:       }
    498:       else {
    499:         bi = x[i] - l[i];
    500:         ai = norm(bi, f[i]);
    501:         ai = TaoMax(TAO_EPSILON, ai);

    503:         da[i] = bi / ai - 1;
    504:         db[i] = f[i] / ai - 1;
    505:       }
    506:     }
    507:     else if (l[i] == u[i]) {
    508:       da[i] = -1;
    509:       db[i] = 0;
    510:     }
    511:     else {
    512:       if (db[i] >= 1) {
    513:         ai = norm(1, t2[i]);

    515:         ci = 1 / ai + 1;
    516:         di = t2[i] / ai + 1;
    517:       }
    518:       else {
    519:         bi = x[i] - u[i];
    520:         ai = norm(bi, f[i]);
    521:         ai = TaoMax(TAO_EPSILON, ai);

    523:         ci = bi / ai + 1;
    524:         di = f[i] / ai + 1;
    525:       }

    527:       if (da[i] >= 1) {
    528:         bi = ci + di*t2[i];
    529:         ai = norm(1, bi);

    531:         bi = bi / ai - 1;
    532:         ai = 1 / ai - 1;
    533:       }
    534:       else {
    535:         ei = fischer(u[i] - x[i], -f[i]);
    536:         ai = norm(x[i] - l[i], ei);
    537:         ai = TaoMax(TAO_EPSILON, ai);

    539:         bi = ei / ai - 1;
    540:         ai = (x[i] - l[i]) / ai - 1;
    541:       }

    543:       da[i] = ai + bi*ci;
    544:       db[i] = bi*di;
    545:     }
    546:   }

    548:   info = tda->RestoreArray(&da,&nn5);CHKERRQ(info);
    549:   info = tdb->RestoreArray(&db,&nn6);CHKERRQ(info);

    551:   info = tx->RestoreArray(&x,&nn1);CHKERRQ(info);
    552:   info = tf->RestoreArray(&f,&nn2);CHKERRQ(info);
    553:   info = tl->RestoreArray(&l,&nn3);CHKERRQ(info);
    554:   info = tu->RestoreArray(&u,&nn4);CHKERRQ(info);
    555:   info = tt2->RestoreArray(&t2,&nn8);CHKERRQ(info);
    556:   TaoFunctionReturn(0);
    557: };

    561: inline static double sfischer(double a, double b, double c)
    562: {

    564: #ifdef TODD

    566:   if (TaoAbsScalar(a) > TaoAbsScalar(b)) {
    567:     return sqrt(a*a + b*b + 2.0*c*c) - a - b;
    568:   }
    569:   return sqrt(a*a + b*b + 2.0*c*c) - b - a;

    571: #else

    573:    // Method suggested by Bob Vanderbei

    575:    if (a + b <= 0) {
    576:      return sqrt(a*a + b*b + 2.0*c*c) - (a + b);
    577:    }
    578:    return 2.0*(c*c - a*b) / (sqrt(a*a + b*b + 2.0*c*c) + (a + b));

    580: #endif

    582: }

    586: inline static double snorm(double a, double b, double c)
    587: {
    588:   return sqrt(a*a + b*b + 2.0*c*c);
    589: }

    593: /*@C
    594:    SD_Fischer - Calculates an element of the B-subdifferential of the
    595:    smoothed Fischer-Burmeister function for complementarity problems.
    596:
    597:    Input Parameters:
    598: +  this - the jacobian of tf at tx
    599: .  tx - current point
    600: .  tf - function evaluated at tx
    601: .  tl - lower bounds
    602: .  tu - upper bounds
    603: .  mu - smoothing parameter
    604: .  tt1 - work vector
    605: -  tt2 - work vector

    607:    Output Parameter:
    608: +  tda - diagonal perturbation component of the result
    609: .  tdb - row scaling component of the result
    610: -  tdm - derivative with respect to scaling parameter

    612:    Level: intermediate

    614: .seealso TaoVec::SFischer()
    615: @*/
    616: int TaoMat::D_SFischer(TaoVec *tx, TaoVec *tf,
    617:                        TaoVec *tl, TaoVec *tu, double mu,
    618:                        TaoVec *tt1, TaoVec *tt2,
    619:                        TaoVec *tda, TaoVec *tdb, TaoVec *tdm)
    620: {
    621:   int i, nn1, nn2, nn3, nn4, nn5, nn6, nn7, info;
    622:   TaoTruth flag;
    623:   TaoScalar *x, *f, *l, *u, *da, *db, *dm;
    624:   TaoScalar ai, bi, ci, di, ei, fi;

    626:   TaoFunctionBegin;

    628:   if ((mu >= -TAO_EPSILON) && (mu <= TAO_EPSILON)) {
    629:     tdm->SetToZero();
    630:     D_Fischer(tx, tf, tl, tu, tt1, tt2, tda, tdb);
    631:   }
    632:   else {
    633:     info = this->Compatible(tx, tx, &flag); CHKERRQ(info);

    635:     info = tx->GetArray(&x, &nn1); CHKERRQ(info);
    636:     info = tf->GetArray(&f, &nn2); CHKERRQ(info);
    637:     info = tl->GetArray(&l, &nn3); CHKERRQ(info);
    638:     info = tu->GetArray(&u, &nn4); CHKERRQ(info);
    639:     info = tda->GetArray(&da, &nn5); CHKERRQ(info);
    640:     info = tdb->GetArray(&db, &nn6); CHKERRQ(info);
    641:     info = tdm->GetArray(&dm, &nn7); CHKERRQ(info);

    643:     if (nn1!=nn2 || nn2!=nn3 || nn3!=nn4 || nn4!=nn5 || nn5!=nn6 || nn6!=nn7) {
    644:       TaoFunctionReturn(1);
    645:     }

    647:     for (i = 0; i < nn1; ++i) {
    648:       if ((l[i] <= -TAO_INFINITY) && (u[i] >= TAO_INFINITY)) {
    649:         da[i] = -mu;
    650:         db[i] = -1;
    651:         dm[i] = -x[i];
    652:       }
    653:       else if (l[i] <= -TAO_INFINITY) {
    654:         bi = u[i] - x[i];
    655:         ai = snorm(bi, f[i], mu);
    656:         ai = TaoMax(TAO_EPSILON, ai);

    658:         da[i] = bi / ai - 1;
    659:         db[i] = -f[i] / ai - 1;
    660:         dm[i] = 2.0 * mu / ai;
    661:       }
    662:       else if (u[i] >=  TAO_INFINITY) {
    663:         bi = x[i] - l[i];
    664:         ai = snorm(bi, f[i], mu);
    665:         ai = TaoMax(TAO_EPSILON, ai);

    667:         da[i] = bi / ai - 1;
    668:         db[i] = f[i] / ai - 1;
    669:         dm[i] = 2.0 * mu / ai;
    670:       }
    671:       else if (l[i] == u[i]) {
    672:         da[i] = -1;
    673:         db[i] = 0;
    674:         dm[i] = 0;
    675:       }
    676:       else {
    677:         bi = x[i] - u[i];
    678:         ai = snorm(bi, f[i], mu);
    679:         ai = TaoMax(TAO_EPSILON, ai);
    680:
    681:         ci = bi / ai + 1;
    682:         di = f[i] / ai + 1;
    683:         fi = 2.0 * mu / ai;

    685:         ei = sfischer(u[i] - x[i], -f[i], mu);
    686:         ai = snorm(x[i] - l[i], ei, mu);
    687:         ai = TaoMax(TAO_EPSILON, ai);
    688:
    689:         bi = ei / ai - 1;
    690:         ei = 2.0 * mu / ei;
    691:         ai = (x[i] - l[i]) / ai - 1;
    692:
    693:         da[i] = ai + bi*ci;
    694:         db[i] = bi*di;
    695:         dm[i] = ei + bi*fi;
    696:       }
    697:     }

    699:     info = tx->RestoreArray(&x, &nn1); CHKERRQ(info);
    700:     info = tf->RestoreArray(&f, &nn2); CHKERRQ(info);
    701:     info = tl->RestoreArray(&l, &nn3); CHKERRQ(info);
    702:     info = tu->RestoreArray(&u, &nn4); CHKERRQ(info);
    703:     info = tda->RestoreArray(&da, &nn5); CHKERRQ(info);
    704:     info = tdb->RestoreArray(&db, &nn6); CHKERRQ(info);
    705:     info = tdm->RestoreArray(&dm, &nn7); CHKERRQ(info);
    706:   }
    707:   TaoFunctionReturn(0);
    708: }

    712: /*@C
    713:   View - Views the contents of this matrix.

    715:    Input Parameter:

    717:    Level: intermediate
    718: @*/
    719: int TaoMat::View(){
    720:   TaoFunctionBegin;
    721:   SETERRQ(56,"Operation not defined");
    722:   /* TaoFunctionReturn(1); */
    723: };

    727: /*@C
    728:    Norm1 - Computes the 1-norm of the matrix.

    730:    Output Parameter:
    731: .  nm -  matrix 1-norm value

    733:    Level: intermediate
    734: @*/
    735: int TaoMat::Norm1(double *nm){
    736:   TaoFunctionBegin;
    737:   SETERRQ(56,"Operation not defined");
    738:   /* TaoFunctionReturn(1); */
    739: }