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: }