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