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) {
 19:     delete TM;
 20:     TM=0;
 21:   }
 22:   TaoFunctionReturn(0);
 23: }

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

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

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

 36:    Level: developer
 37: @*/
 38: int TaoMat::Compatible(TaoVec *xx, TaoVec *yy, TaoTruth *flag){
 39:   TaoFunctionBegin;
 40:   *flag=TAO_FALSE;
 41:   TaoFunctionReturn(0);
 42: }


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

 50:    Output Parameter:
 51: +  m -  the number of rows
 52: -  n - the number of columns

 54:    Level: intermediate

 56: @*/
 57: int TaoMat::GetDimensions( TaoInt* m, TaoInt* n ){
 58:   TaoFunctionBegin;
 59:   SETERRQ(56,"Operation not defined");
 60:   /* TaoFunctionReturn(1); */
 61: }

 65: /*@C
 66:    Multiply - Computes  ty = this * tx.

 68:    Input Parameter:
 69: .  tx -  the vector to be multiplied by this matrix.

 71:    Output Parameter:
 72: .  ty -  the destination vector


 75:    Level: intermediate

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

 79: @*/
 80: int TaoMat::Multiply(TaoVec* tx,TaoVec* ty){
 81:   TaoFunctionBegin;
 82:   SETERRQ(56,"Operation not defined");
 83:   /* TaoFunctionReturn(1); */
 84: }


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

 92:    Input Parameter:
 93: .  tx -  the vector to be multiplied by this matrix.

 95:    Output Parameter:
 96: .  ty -  the destination vector


 99:    Level: intermediate

101: .seealso TaoMat::Multiply(), TaoMat::MultiplyTransposeAdd()
102: @*/
103: int TaoMat::MultiplyTranspose(TaoVec* tx,TaoVec* ty){
104:   TaoFunctionBegin;
105:   SETERRQ(56,"Operation not defined");
106:   /* TaoFunctionReturn(1); */
107: }


112: /*@C
113:    SetDiagonal - Sets the diagonal elements of this matrix with the elements
114:    of the vector.

116:    Input Parameter:
117: .  tv -  the vector containing the diagonal elements

119:    Level: intermediate

121: .seealso TaoMat::AddDiagonal(),TaoMat::ShiftDiagonal()
122: @*/
123: int TaoMat::SetDiagonal(TaoVec* tv){
124:   TaoFunctionBegin;
125:   SETERRQ(56,"Operation not defined");
126:   /* TaoFunctionReturn(1); */
127: }

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

134:    Input Parameter:
135: .  tv -  the vector containing the diagonal elements


138:    Level: intermediate

140: .seealso TaoMat::SetDiagonal(),TaoMat::ShiftDiagonal()
141: @*/
142: int TaoMat::AddDiagonal(TaoVec* tv){
143:   TaoFunctionBegin;
144:   SETERRQ(56,"Operation not defined");
145:   /* TaoFunctionReturn(1); */
146: }

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

153:    Output Parameter:
154: .  tv -  the vector containing the diagonal elements

156:    Level: intermediate

158: .seealso TaoMat::SetDiagonal()
159: @*/
160: int TaoMat::GetDiagonal(TaoVec* tv){
161:   TaoFunctionBegin;
162:   SETERRQ(56,"Operation not defined");
163:   /* TaoFunctionReturn(1); */
164: }

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

171:    Input Parameter:
172: .  c -  the constant

174:    Level: intermediate

176: .seealso TaoMat::SetDiagonal(), TaoMat::AddDiagonal()
177: @*/
178: int TaoMat::ShiftDiagonal(double c){
179:   TaoFunctionBegin;
180:   SETERRQ(56,"Operation not defined");
181:   /* TaoFunctionReturn(1); */
182: }

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

191:    Input:

193:    Note:
194:    This routine is optional.  A linear solver object can also be used
195:    to implement this operation.

197:    Level: advanced

199: .seealso TaoMat::Solve()
200: @*/
201: int TaoMat::Presolve(){
202:   TaoFunctionBegin;
203:   TaoFunctionReturn(0);
204: }

208: /*@C
209:    Solve - Solves the linear system $this tx = tb$.

211:    Input Parameter:
212: .  tb -  right hand side

214:    Output Parameter:
215: +  tx -  solution vector
216: -  tt -  TAO_TRUE if solved to prescribed accuracy and TAO_FALSE otherwise

218:    Level: advanced

220:    Note:
221:    This routine is optional.  A linear solver object can also be used
222:    to implement this operation.

224: .seealso TaoApplication::GetLinearSolver(),TaoMat::Presolve()
225: @*/
226: int TaoMat::Solve(TaoVec* tb, TaoVec *tx, TaoTruth *tt){
227:   TaoFunctionBegin;
228:   SETERRQ(56,"No TAO linear solver has been set.");
229:   /* TaoFunctionReturn(1); */
230: }


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

238:    Input Parameter:
239: +  row -  the rows of this matrix to be extracted
240: -  column -  the columns of this matrix to be extracted

242:    Output Parameter:
243: -  M - the new matrix


246: .seealso TaoMat::SetReducedMatrix()

248:    Level: intermediate
249: @*/
250: int TaoMat::CreateReducedMatrix(TaoIndexSet* row,TaoIndexSet* col, TaoMat **M){
251:   TaoFunctionBegin;
252:   SETERRQ(56,"Operation not defined");
253:   /* TaoFunctionReturn(1); */
254: };

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

261:    Input Parameter:
262: +  row -  the rows of this matrix to be extracted
263: -  col -  the columns of this matrix to be extracted

265:    Output Parameter:
266: .  M - the full new matrix


269: .seealso TaoMat::CreateReducedMatrix()

271:    Level: intermediate
272: @*/
273: int TaoMat::SetReducedMatrix(TaoMat *M, TaoIndexSet* row,TaoIndexSet* col){
274:   TaoFunctionBegin;
275:   SETERRQ(56,"Operation not defined");
276:   /* TaoFunctionReturn(1); */
277: };


282: /*@C
283:   RowScale - Scales the rows of this matrix.

285:    Input Parameter:
286: .  scale - the scaling parameters

288:    Note:
289:    This method can also be called by using a pointer to TaoVec as
290:    input parameters.

292:    Level: intermediate
293: @*/
294: int TaoMat::RowScale(TaoVec* scale){
295:   TaoFunctionBegin;
296:   SETERRQ(56,"Operation not defined");
297:   /* TaoFunctionReturn(1); */
298: };

302: /*@C
303:   ColScale - Scales the columns of this matrix.

305:    Input Parameter:
306: .  scale - the scaling parameters

308:    Level: intermediate
309: @*/
310: int TaoMat::ColScale(TaoVec* scale){
311:   TaoFunctionBegin;
312:   SETERRQ(56,"Operation not defined");
313:   /* TaoFunctionReturn(1); */
314: };

318: inline static double fischer(double a, double b)
319: {

321: #ifdef TODD

323:   if (TaoAbsScalar(a) > TaoAbsScalar(b)) {
324:     return sqrt(a*a + b*b) - a - b;
325:   }
326:   return sqrt(a*a + b*b) - b - a;

328: #else

330:    // Method suggested by Bob Vanderbei

332:    if (a + b <= 0) {
333:      return sqrt(a*a + b*b) - (a + b);
334:    }
335:    return -2.0*a*b / (sqrt(a*a + b*b) + (a + b));

337: #endif

339: }

343: inline static double norm(double a, double b)
344: {
345:   return sqrt(a*a + b*b);
346: }

350: /*@C
351:    D_Fischer - Calculates an element of the B-subdifferential of the 
352:    Fischer-Burmeister function for complementarity problems.
353:   
354:    Input Parameters:   
355: +  this - the jacobian of tf at tx
356: .  tx - current point
357: .  tf - function evaluated at tx
358: .  tl - lower bounds
359: .  tu - upper bounds
360: .  tt1 - work vector
361: -  tt2 - work vector

363:    Output Parameters:
364: +  tda - diagonal perturbation component of the result
365: -  tdb - row scaling component of the result

367:    Level: intermediate

369: .seealso TaoVec::Fischer()
370: @*/
371: int TaoMat::D_Fischer(TaoVec *tx, TaoVec *tf, TaoVec *tl, TaoVec *tu, 
372:                       TaoVec *tt1, TaoVec *tt2, TaoVec *tda, TaoVec *tdb )
373: {
374:   int i, info;
375:   TaoInt nn1,nn2,nn3,nn4,nn5,nn6,nn7,nn8;
376:   TaoTruth flag;
377:   TaoScalar *x,*f,*l,*u,*da,*db,*t1,*t2;
378:   TaoScalar ai,bi,ci,di,ei;

380:   TaoFunctionBegin;

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

385:   info = tx->GetArray(&x,&nn1);CHKERRQ(info);
386:   info = tf->GetArray(&f,&nn2);CHKERRQ(info);
387:   info = tl->GetArray(&l,&nn3);CHKERRQ(info);
388:   info = tu->GetArray(&u,&nn4);CHKERRQ(info);
389:   info = tda->GetArray(&da,&nn5);CHKERRQ(info);
390:   info = tdb->GetArray(&db,&nn6);CHKERRQ(info);
391:   info = tt1->GetArray(&t1,&nn7);CHKERRQ(info);
392:   info = tt2->GetArray(&t2,&nn8);CHKERRQ(info);

394:   if (nn1!=nn2 || nn2!=nn3 || nn3!=nn4 || nn4!=nn5 || 
395:       nn5!=nn6 || nn6!=nn7 || nn7!=nn8) {
396:     TaoFunctionReturn(1);
397:   }

399:   for (i = 0; i < nn1; i++) {
400:     da[i] = 0;
401:     db[i] = 0;
402:     t1[i] = 0;

404:     if (TaoAbsScalar(f[i]) <= TAO_EPSILON) {
405:       if (l[i] > -TAO_INFINITY && TaoAbsScalar(x[i] - l[i]) <= TAO_EPSILON) {
406:         t1[i] = 1;
407:         da[i] = 1;
408:       }

410:       if (u[i] <  TAO_INFINITY && TaoAbsScalar(u[i] - x[i]) <= TAO_EPSILON) {
411:         t1[i] = 1;
412:         db[i] = 1;
413:       }
414:     }
415:   }

417:   info = tt1->RestoreArray(&t1,&nn7);CHKERRQ(info);
418:   info = tt2->RestoreArray(&t2,&nn8);CHKERRQ(info);
419:   info = this->Multiply(tt1, tt2); CHKERRQ(info);
420:   info = tt2->GetArray(&t2,&nn8); CHKERRQ(info);

422:   for (i = 0; i < nn1; i++) {
423:     if ((l[i] <= -TAO_INFINITY) && (u[i] >= TAO_INFINITY)) {
424:       da[i] = 0;
425:       db[i] = -1;
426:     } 
427:     else if (l[i] <= -TAO_INFINITY) {
428:       if (db[i] >= 1) {
429:         ai = norm(1, t2[i]);

431:         da[i] = -1/ai - 1;
432:         db[i] = -t2[i]/ai - 1;
433:       } 
434:       else {
435:         bi = u[i] - x[i];
436:         ai = norm(bi, f[i]);
437:         ai = TaoMax(TAO_EPSILON, ai);

439:         da[i] = bi / ai - 1;
440:         db[i] = -f[i] / ai - 1;
441:       }
442:     } 
443:     else if (u[i] >=  TAO_INFINITY) {
444:       if (da[i] >= 1) {
445:         ai = norm(1, t2[i]);

447:         da[i] = 1 / ai - 1;
448:         db[i] = t2[i] / ai - 1;
449:       } 
450:       else {
451:         bi = x[i] - l[i];
452:         ai = norm(bi, f[i]);
453:         ai = TaoMax(TAO_EPSILON, ai);

455:         da[i] = bi / ai - 1;
456:         db[i] = f[i] / ai - 1;
457:       }
458:     } 
459:     else if (l[i] == u[i]) {
460:       da[i] = -1;
461:       db[i] = 0;
462:     } 
463:     else {
464:       if (db[i] >= 1) {
465:         ai = norm(1, t2[i]);

467:         ci = 1 / ai + 1;
468:         di = t2[i] / ai + 1;
469:       } 
470:       else {
471:         bi = x[i] - u[i];
472:         ai = norm(bi, f[i]);
473:         ai = TaoMax(TAO_EPSILON, ai);

475:         ci = bi / ai + 1;
476:         di = f[i] / ai + 1;
477:       }

479:       if (da[i] >= 1) {
480:         bi = ci + di*t2[i];
481:         ai = norm(1, bi);

483:         bi = bi / ai - 1;
484:         ai = 1 / ai - 1;
485:       } 
486:       else {
487:         ei = fischer(u[i] - x[i], -f[i]);
488:         ai = norm(x[i] - l[i], ei);
489:         ai = TaoMax(TAO_EPSILON, ai);

491:         bi = ei / ai - 1;
492:         ai = (x[i] - l[i]) / ai - 1;
493:       }

495:       da[i] = ai + bi*ci;
496:       db[i] = bi*di;
497:     }
498:   }

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

503:   info = tx->RestoreArray(&x,&nn1);CHKERRQ(info);
504:   info = tf->RestoreArray(&f,&nn2);CHKERRQ(info);
505:   info = tl->RestoreArray(&l,&nn3);CHKERRQ(info);
506:   info = tu->RestoreArray(&u,&nn4);CHKERRQ(info);
507:   info = tt2->RestoreArray(&t2,&nn8);CHKERRQ(info);
508:   TaoFunctionReturn(0);
509: };

513: inline static double sfischer(double a, double b, double c)
514: {

516: #ifdef TODD

518:   if (TaoAbsScalar(a) > TaoAbsScalar(b)) {
519:     return sqrt(a*a + b*b + 2.0*c*c) - a - b;
520:   }
521:   return sqrt(a*a + b*b + 2.0*c*c) - b - a;

523: #else

525:    // Method suggested by Bob Vanderbei

527:    if (a + b <= 0) {
528:      return sqrt(a*a + b*b + 2.0*c*c) - (a + b);
529:    }
530:    return 2.0*(c*c - a*b) / (sqrt(a*a + b*b + 2.0*c*c) + (a + b));

532: #endif

534: }

538: inline static double snorm(double a, double b, double c)
539: {
540:   return sqrt(a*a + b*b + 2.0*c*c);
541: }

545: /*@C
546:    D_SFischer - Calculates an element of the B-subdifferential of the
547:    smoothed Fischer-Burmeister function for complementarity problems.
548:  
549:    Input Parameters: 
550: +  this - the jacobian of tf at tx
551: .  tx - current point
552: .  tf - function evaluated at tx
553: .  tl - lower bounds
554: .  tu - upper bounds
555: .  mu - smoothing parameter
556: .  tt1 - work vector
557: -  tt2 - work vector

559:    Output Parameter: 
560: +  tda - diagonal perturbation component of the result
561: .  tdb - row scaling component of the result
562: -  tdm - derivative with respect to scaling parameter

564:    Level: intermediate

566: .seealso TaoVec::SFischer()
567: @*/
568: int TaoMat::D_SFischer(TaoVec *tx, TaoVec *tf, 
569:                        TaoVec *tl, TaoVec *tu, double mu, 
570:                        TaoVec *tt1, TaoVec *tt2, 
571:                        TaoVec *tda, TaoVec *tdb, TaoVec *tdm)
572: {
573:   int i, info;
574:   TaoInt nn1, nn2, nn3, nn4, nn5, nn6, nn7;
575:   TaoTruth flag;
576:   TaoScalar *x, *f, *l, *u, *da, *db, *dm;
577:   TaoScalar ai, bi, ci, di, ei, fi;

579:   TaoFunctionBegin;

581:   if ((mu >= -TAO_EPSILON) && (mu <= TAO_EPSILON)) {
582:     tdm->SetToZero();
583:     D_Fischer(tx, tf, tl, tu, tt1, tt2, tda, tdb);
584:   } 
585:   else {
586:     info = this->Compatible(tx, tx, &flag); CHKERRQ(info);

588:     info = tx->GetArray(&x, &nn1); CHKERRQ(info);
589:     info = tf->GetArray(&f, &nn2); CHKERRQ(info);
590:     info = tl->GetArray(&l, &nn3); CHKERRQ(info);
591:     info = tu->GetArray(&u, &nn4); CHKERRQ(info);
592:     info = tda->GetArray(&da, &nn5); CHKERRQ(info);
593:     info = tdb->GetArray(&db, &nn6); CHKERRQ(info);
594:     info = tdm->GetArray(&dm, &nn7); CHKERRQ(info);

596:     if (nn1!=nn2 || nn2!=nn3 || nn3!=nn4 || nn4!=nn5 || nn5!=nn6 || nn6!=nn7) {
597:       TaoFunctionReturn(1);
598:     }

600:     for (i = 0; i < nn1; ++i) {
601:       if ((l[i] <= -TAO_INFINITY) && (u[i] >= TAO_INFINITY)) {
602:         da[i] = -mu;
603:         db[i] = -1;
604:         dm[i] = -x[i];
605:       } 
606:       else if (l[i] <= -TAO_INFINITY) {
607:         bi = u[i] - x[i];
608:         ai = snorm(bi, f[i], mu);
609:         ai = TaoMax(TAO_EPSILON, ai);

611:         da[i] = bi / ai - 1;
612:         db[i] = -f[i] / ai - 1;
613:         dm[i] = 2.0 * mu / ai;
614:       } 
615:       else if (u[i] >=  TAO_INFINITY) {
616:         bi = x[i] - l[i];
617:         ai = snorm(bi, f[i], mu);
618:         ai = TaoMax(TAO_EPSILON, ai);

620:         da[i] = bi / ai - 1;
621:         db[i] = f[i] / ai - 1;
622:         dm[i] = 2.0 * mu / ai;
623:       } 
624:       else if (l[i] == u[i]) {
625:         da[i] = -1;
626:         db[i] = 0;
627:         dm[i] = 0;
628:       } 
629:       else {
630:         bi = x[i] - u[i];
631:         ai = snorm(bi, f[i], mu);
632:         ai = TaoMax(TAO_EPSILON, ai);
633:   
634:         ci = bi / ai + 1;
635:         di = f[i] / ai + 1;
636:         fi = 2.0 * mu / ai;

638:         ei = sfischer(u[i] - x[i], -f[i], mu);
639:         ai = snorm(x[i] - l[i], ei, mu);
640:         ai = TaoMax(TAO_EPSILON, ai);
641:   
642:         bi = ei / ai - 1;
643:         ei = 2.0 * mu / ei;
644:         ai = (x[i] - l[i]) / ai - 1;
645:   
646:         da[i] = ai + bi*ci;
647:         db[i] = bi*di;
648:         dm[i] = ei + bi*fi;
649:       }
650:     }

652:     info = tx->RestoreArray(&x, &nn1); CHKERRQ(info);
653:     info = tf->RestoreArray(&f, &nn2); CHKERRQ(info);
654:     info = tl->RestoreArray(&l, &nn3); CHKERRQ(info);
655:     info = tu->RestoreArray(&u, &nn4); CHKERRQ(info);
656:     info = tda->RestoreArray(&da, &nn5); CHKERRQ(info);
657:     info = tdb->RestoreArray(&db, &nn6); CHKERRQ(info);
658:     info = tdm->RestoreArray(&dm, &nn7); CHKERRQ(info);
659:   }
660:   TaoFunctionReturn(0);
661: }

665: /*@C
666:   View - Views the contents of this matrix.

668:    Input Parameter:

670:    Level: intermediate
671: @*/
672: int TaoMat::View(){
673:   TaoFunctionBegin;
674:   SETERRQ(56,"Operation not defined");
675:   /* TaoFunctionReturn(1); */
676: };

680: /*@C
681:    Norm1 - Computes the 1-norm of the matrix.

683:    Output Parameter:
684: .  nm -  matrix 1-norm value

686:    Level: intermediate
687: @*/
688: int TaoMat::Norm1(double *nm){
689:   TaoFunctionBegin;
690:   SETERRQ(56,"Operation not defined");
691:   /* TaoFunctionReturn(1); */
692: }