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