Actual source code: taovec.c
1: #include "tao_general.h" /*I "tao_solver.h" I*/
2: #include "taovec.h"
6: /*@C
7: TaoVecDestroy - Destroys the TaoVec object.
9: Input Parameter:
10: . vv - the vector
12: Level: beginner
13: @*/
14: int TaoVecDestroy( TaoVec* vv){
15: TaoFunctionBegin;
16: if (vv!=TAO_NULL && vv!=0){
17: delete vv;
18: }
19: vv=TAO_NULL;
20: TaoFunctionReturn(0);
21: }
25: /*@C
26: CloneVecs - Creates an array of pointers to new TaoVec objects. The new
27: objects have the same structure as this one.
29: Input Parameter:
30: . nn - number of new vectors
32: Output Parameter:
33: . tvs - pointer to array TaoVec pointers
35: .seealso TaoVec::Clone()
37: Level: intermediate
38: @*/
39: int TaoVec::CloneVecs(TaoInt nn, TaoVec***tvs){
40: int info;
41: TaoInt i;
42: TaoVec ** ntv;
43: TaoFunctionBegin;
44: ntv = new TaoVec*[nn];
45: for (i=0;i<nn;i++){
46: info = this->Clone(&ntv[i]);CHKERRQ(info);
47: }
48: *tvs=ntv;
49: TaoFunctionReturn(0);
50: }
54: /*@C
55: DestroyVecs - Destroys an array TaoVec objects.
57: Input Parameter:
58: + nn - number of new vectors
59: - tvs - pointer to array TaoVec pointers
61: Level: advanced
63: .seealso TaoVec::CloneVecs()
64: @*/
65: int TaoVec::DestroyVecs(TaoInt nn, TaoVec**ntv){
66: int info;
67: TaoInt i;
68: TaoFunctionBegin;
69: for (i=0;i<nn;i++){
70: info = TaoVecDestroy( ntv[i] );CHKERRQ(info);
71: }
72: delete[] ntv;
73: TaoFunctionReturn(0);
74: }
78: /*@C
79: Clone - Creates a new TaoVec object with the same structure as this one. It does not
80: copy the value to the new vector.
82: Input:
83: . vv - address of a pointer to a TaoVec
85: Output Parameter:
86: . vv - address of a pointer to new TaoVec object
88: .seealso TaoVec::CloneVecs(), TaoVec::CopyFrom()
90: Level: intermediate
91: @*/
92: int TaoVec::Clone( TaoVec* *vv ){
93: TaoFunctionBegin;
94: SETERRQ(56,"Operation not defined");
95: /* TaoFunctionReturn(1); */
96: }
100: /*@C
101: Compatible - Determines whether this vector belongs to the same space as another,
102: and operations such as inner product and sum are well defined.
104: Input Parameter:
105: . vv - TAO vector to which to the comparison is made
107: Output Value:
108: . flag - TAO_TRUE if the two vectors are Compatible and TAO_FALSE otherwise.
110: Level: advanced
112: .seealso TaoVec::GetDimension()
113: @*/
114: int TaoVec::Compatible(TaoVec* vv, TaoTruth *flag){
115: TaoFunctionBegin;
116: if (!flag){
117: TaoFunctionReturn(1);
118: }
119: *flag=TAO_FALSE;
120: TaoFunctionReturn(0);
121: }
125: /*@C
126: SetToConstant - Sets each element of this vector equal to the specified constant.
128: Input Parameter:
129: . c - a constant
131: Level: intermediate
133: .seealso TaoVec::Scale()
134: @*/
135: int TaoVec::SetToConstant( double c ){
136: int info;
137: TaoInt nn,i;
138: TaoScalar *tptr;
140: TaoFunctionBegin;
141: info = this->GetArray(&tptr,&nn);CHKERRQ(info);
142: for (i=0;i<nn;i++){ tptr[i]=c; }
143: info = this->RestoreArray(&tptr,&nn);CHKERRQ(info);
144: TaoFunctionReturn(0);
145: }
149: /*@C
150: SetToZero - Sets each element of this vector equal to zero.
152: Input Parameters: none
154: Level: intermediate
156: .seealso TaoVec::SetToConstant()
157: @*/
158: int TaoVec::SetToZero(){
159: int info;
160: TaoInt nn,i;
161: TaoScalar *tptr;
163: TaoFunctionBegin;
164: info = this->GetArray(&tptr,&nn);CHKERRQ(info);
165: for (i=0;i<nn;i++){ tptr[i]=0; }
166: info = this->RestoreArray(&tptr,&nn);CHKERRQ(info);
167: TaoFunctionReturn(0);
168: }
172: /*@C
173: CopyFrom - Copies the contents of one vector into this vector.
175: Input Parameter:
176: . vv - A TaoVec from which the contents will be copied.
178: Level: intermediate
180: .seealso TaoVec::Axpy(), TaoVec::ScaleCopyFrom()
181: @*/
182: int TaoVec::CopyFrom( TaoVec* vv ){
183: int info;
184: TaoInt nn1,nn2,i;
185: TaoScalar *tptr1,*tptr2;
187: TaoFunctionBegin;
188: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
189: if (vv!=this){
190: info = vv->GetArray(&tptr2,&nn2);CHKERRQ(info);
191: if (nn1!=nn2) {TaoFunctionReturn(1);}
192: for (i=0;i<nn1;i++){ tptr1[i]=tptr2[i]; }
193: info = vv->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
194: }
195: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
196: TaoFunctionReturn(0);
197: }
201: /*@C
202: ScaleCopyFrom - Copies the contents of one vector into this vector and scales it.
204: Input Parameter:
205: + a - the scalar
206: - vv - A TaoVec from which the contents will be copied.
208: Level: intermediate
210: .seealso TaoVec::Axpy(), TaoVec::Aypx()
211: @*/
212: int TaoVec::ScaleCopyFrom( double a, TaoVec* vv ){
213: int info;
214: TaoInt nn1,nn2,i;
215: TaoScalar *tptr1,*tptr2;
217: TaoFunctionBegin;
218: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
219: if (vv!=this){
220: info = vv->GetArray(&tptr2,&nn2);CHKERRQ(info);
221: if (nn1!=nn2) {TaoFunctionReturn(1);}
222: for (i=0;i<nn1;i++){ tptr1[i]=a*tptr2[i]; }
223: info = vv->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
224: } else {
225: this->Scale(a);
226: }
227: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
228: TaoFunctionReturn(0);
229: }
233: /*@C
234: NormInfinity - Computes the infinity norm of this vector.
236: Ouput Parameter:
237: . vnorm - the infinity norm of this vector
239: Level: intermediate
241: .seealso TaoVec::Norm1(), TaoVec::Norm2()
242: @*/
243: int TaoVec::NormInfinity(double *vnorm){
244: int info;
245: TaoInt nn,i;
246: TaoScalar dd=0, *v;
248: TaoFunctionBegin;
249: info = this->GetArray(&v,&nn);CHKERRQ(info);
250: for (i=0;i<nn;i++){
251: if (v[i]<0 && v[i]<-dd) dd=-v[i];
252: else if (v[i]>0 && v[i]>dd) dd=v[i];
253: }
254: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
255: *vnorm=dd;
256: TaoFunctionReturn(0);
257: }
261: /*@C
262: Norm1 - Computes the one-norm of this vector.
264: Ouput Parameter:
265: . vnorm - the one-norm of this vector
267: Level: intermediate
269: .seealso TaoVec::NormInfinity(), TaoVec::Norm2()
270: @*/
271: int TaoVec::Norm1(double *vnorm){
272: int info;
273: TaoInt nn,i;
274: TaoScalar dd=0, *v;
276: TaoFunctionBegin;
277: info = this->GetArray(&v,&nn);CHKERRQ(info);
278: for (i=0;i<nn;i++){
279: if (v[i]<0) dd-=v[i];
280: else dd+=v[i];
281: }
282: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
283: *vnorm = dd;
284: TaoFunctionReturn(0);
285: }
289: /*@C
290: Norm2 - Compute the two-norm of this vector.
292: Ouput Parameter:
293: . vnorm - the two-norm of this vector
295: Level: intermediate
297: .seealso TaoVec::Norm1(), TaoVec::NormInfinity(), TaoVec::Norm2squared()
298: @*/
299: int TaoVec::Norm2(double *vnorm){
300: int info;
301: TaoInt nn,i;
302: TaoScalar dd=0, *v;
304: TaoFunctionBegin;
305: info = this->GetArray(&v,&nn);CHKERRQ(info);
306: for (i=0;i<nn;i++) dd+=v[i]*v[i];
307: *vnorm=sqrt(dd);
308: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
309: TaoFunctionReturn(0);
310: }
314: /*@C
315: Norm2squared - Computes the square of the two norm of this vector.
317: Ouput Parameter:
318: . vnorm2 - the square of the two norm of this vector
320: Level: intermediate
322: .seealso TaoVec::Norm2()
323: @*/
324: int TaoVec::Norm2squared(double *vnorm2){
325: int info;
326: TaoInt nn,i;
327: TaoScalar dd=0, *v;
329: TaoFunctionBegin;
330: info = this->GetArray(&v,&nn);CHKERRQ(info);
331: for (i=0;i<nn;i++) dd+=v[i]*v[i];
332: *vnorm2=dd;
333: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
334: TaoFunctionReturn(0);
335: }
339: /*@C
340: Scale - Multiplies this vector by a scalar.
342: Input Parameter:
343: . alpha - the scalar
345: Level: intermediate
347: .seealso TaoVec::SetToConstant(), TaoVec::Aypx()
348: @*/
349: int TaoVec::Scale( double alpha ){
350: int info;
351: TaoInt nn,i;
352: TaoScalar *v;
354: TaoFunctionBegin;
355: info = this->GetArray(&v,&nn);CHKERRQ(info);
356: for (i=0;i<nn;i++) v[i]*=alpha;
357: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
358: TaoFunctionReturn(0);
359: }
363: /*@C
364: Axpy - Adds a scalar multiple of a vector to this vector. (this += alpha * xx)
366: Input Parameter:
367: + alpha - the scalar
368: - xx - the vector
370: Level: intermediate
372: .seealso TaoVec::CopyFrom(), TaoVec::Aypx()
373: @*/
374: int TaoVec::Axpy( double alpha, TaoVec* xx ){
375: int info;
376: TaoInt nn1,nn2,i;
377: TaoScalar *tptr1,*tptr2;
379: TaoFunctionBegin;
380: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
381: info = xx->GetArray(&tptr2,&nn2);CHKERRQ(info);
382: if (nn1!=nn2) {TaoFunctionReturn(1);}
383: for (i=0;i<nn1;i++){ tptr1[i]+= alpha * tptr2[i]; }
384: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
385: info = xx->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
386: TaoFunctionReturn(0);
387: }
391: /*@C
392: Axpby - Adds a scalar multiple of a vector to a multiple of this vector. (this=alpha*xx + beta*this)
394: Input Parameter:
395: + alpha - the scalar of tx
396: . xx - the vector
397: - beta - the scalar multiple of this vector
399: Level: intermediate
401: .seealso TaoVec::Axpy(), TaoVec::Aypx()
402: @*/
403: int TaoVec::Axpby( double alpha, TaoVec* xx, double beta ){
404: int info;
405: TaoInt nn1,nn2,i;
406: TaoScalar *tptr1,*tptr2;
408: TaoFunctionBegin;
409: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
410: info = xx->GetArray(&tptr2,&nn2);CHKERRQ(info);
411: if (nn1!=nn2) {TaoFunctionReturn(1);}
412: for (i=0;i<nn1;i++){ tptr1[i] = beta * tptr1[i] + alpha * tptr2[i]; }
413: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
414: info = xx->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
415: TaoFunctionReturn(0);
416: }
420: /*@C
421: Aypx - Adds a vector to a scalar multiple of this vector. (this=alpha*this+xx)
423: Input Parameter:
424: + alpha - the scalar
425: - xx - the vector
427: Level: intermediate
429: .seealso TaoVec::Scale(), TaoVec::Axpy()
430: @*/
431: int TaoVec::Aypx( double alpha, TaoVec* xx ){
432: int info;
433: TaoInt nn1,nn2,i;
434: TaoScalar *tptr1,*tptr2;
436: TaoFunctionBegin;
437: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
438: info = xx->GetArray(&tptr2,&nn2);CHKERRQ(info);
439: if (nn1!=nn2) {TaoFunctionReturn(1);}
440: for (i=0;i<nn1;i++){ tptr1[i] = alpha * tptr1[i] + tptr2[i]; }
441: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
442: info = xx->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
443: TaoFunctionReturn(0);
444: }
445:
448: /*@C
449: AddConstant - Adds a constant to each element of this vector.
451: Input Parameter:
452: . alpha - the scalar
454: Level: intermediate
456: .seealso TaoVec::SetToConstant(), TaoVec::Axpy()
457: @*/
458: int TaoVec::AddConstant( double alpha ){
459: int info;
460: TaoInt nn,i;
461: TaoScalar *v;
463: TaoFunctionBegin;
464: info = this->GetArray(&v,&nn);CHKERRQ(info);
465: for (i=0;i<nn;i++) v[i]+=alpha;
466: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
467: TaoFunctionReturn(0);
468: }
472: /*@C
473: Dot - Computes the inner product of this vector with another vector.
475: Input Parameter:
476: . vv - another TaoVec object
478: Output Parameter:
479: . vDotv - the inner product of the two vectors
481: Level: intermediate
483: .seealso TaoVec::Norm()
484: @*/
485: int TaoVec::Dot( TaoVec* vv, double *vDotv ){
486: int info;
487: TaoInt nn1,nn2,i;
488: TaoScalar dd=0,*tptr1,*tptr2;
490: TaoFunctionBegin;
491: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
492: info = vv->GetArray(&tptr2,&nn2);CHKERRQ(info);
493: if (nn1!=nn2) {TaoFunctionReturn(1);}
494: for (i=0;i<nn1;i++) dd+=tptr1[i]*tptr2[i];
495: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
496: info = vv->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
497: *vDotv=dd;
498: TaoFunctionReturn(0);
499: }
503: /*@C
504: Negate - Multiplies the elements of this vector by negative one.
506: Input Parameters: none
508: Level: intermediate
510: .seealso TaoVec::Scale()
511: @*/
512: int TaoVec::Negate(){
513: int info;
514: TaoInt nn,i;
515: TaoScalar *v;
517: TaoFunctionBegin;
518: info = this->GetArray(&v,&nn);CHKERRQ(info);
519: for (i=0;i<nn;i++) v[i]=-v[i];
520: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
521: TaoFunctionReturn(0);
522: }
526: /*@C
527: Reciprocal - Sets each element of this vector to its Reciprocal.
529: Input Parameters: none
531: Level: intermediate
533: .seealso TaoVec::PointwiseDivide()
534: @*/
535: int TaoVec::Reciprocal(){
536: int info;
537: TaoInt nn,i;
538: TaoScalar *v;
540: TaoFunctionBegin;
541: info = this->GetArray(&v,&nn);CHKERRQ(info);
542: for (i=0;i<nn;i++){
543: if (v[i]!=0) v[i]= 1.0/v[i];
544: else v[i]=TAO_INFINITY;
545: }
546: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
547: TaoFunctionReturn(0);
548: }
552: /*@C
553: Sqrt - Sets each element of this vector to its square root.
555: Input Parameters: none
557: Level: intermediate
559: @*/
560: int TaoVec::Sqrt(){
561: int info;
562: TaoInt i;
563: TaoInt nn;
564: TaoScalar *v;
566: TaoFunctionBegin;
567: info = this->GetArray(&v,&nn);CHKERRQ(info);
568: for (i=0;i<nn;i++){
569: if (v[i] >= 0) {
570: v[i] = sqrt(v[i]);
571: }
572: else {
573: v[i] = TAO_INFINITY;
574: }
575: }
576: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
577: TaoFunctionReturn(0);
578: }
582: /*@C
583: Pow - Raises each element of this vector to a power.
585: Input Parameters: none
587: Level: intermediate
589: @*/
590: int TaoVec::Pow(double p){
591: int info;
592: TaoInt nn,i;
593: TaoScalar *v;
595: TaoFunctionBegin;
596: info = this->GetArray(&v,&nn);CHKERRQ(info);
597: for (i=0;i<nn;i++){
598: if (v[i] >= 0) {
599: v[i] = pow(v[i], p);
600: }
601: else {
602: v[i] = TAO_INFINITY;
603: }
604: }
605: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
606: TaoFunctionReturn(0);
607: }
611: /*@C
612: GetDimension - Gets the dimension of the vector space where this vector belongs.
614: Output Parameter:
615: . n - the dimension of the vector space
617: Level: intermediate
619: .seealso TaoVec::Compatible()
620: @*/
621: int TaoVec::GetDimension(TaoInt *n){
622: TaoFunctionBegin;
623: SETERRQ(56,"Operation not defined");
624: /* TaoFunctionReturn(1); */
625: }
629: /*@C
630: PointwiseMultiply - Computes the componentwise multiplication of two vectors
631: and stores the result in this vector.
633: Input Parameters:
634: . vv, ww - the two vectors
636: Level: intermediate
638: .seealso TaoVec::PointwiseDivide()
639: @*/
640: int TaoVec::PointwiseMultiply( TaoVec* vv, TaoVec* ww ){
641: int info;
642: TaoInt nn1,nn2,nn3,i;
643: TaoScalar *tptr1,*tptr2,*tptr3;
645: TaoFunctionBegin;
646: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
647: if (this!=vv){
648: info = vv->GetArray(&tptr2,&nn2);CHKERRQ(info);
649: } else {
650: tptr2=tptr1; nn2=nn1;
651: }
652: if (this!=ww && vv!=ww){
653: info = ww->GetArray(&tptr3,&nn3);CHKERRQ(info);
654: } else if (vv==ww){
655: tptr3=tptr2; nn3=nn2;
656: } else {
657: tptr3=tptr1; nn3=nn1;
658: }
659: if (nn1!=nn2 || nn2!=nn3) {TaoFunctionReturn(1);}
660: for (i=0;i<nn1;i++) tptr1[i]=tptr2[i] * tptr3[i];
661: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
662: if (this!=vv){
663: info = vv->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
664: }
665: if (this!=ww && vv!=ww){
666: info = ww->RestoreArray(&tptr3,&nn3);CHKERRQ(info);
667: }
669: TaoFunctionReturn(0);
670: }
674: /*@C
675: PointwiseDivide - Computes the componentwise division of two vectors
676: and stores the result in this vector.
678: Input Parameters:
679: + vv - the vector of numerators
680: - ww - the vector of denominators
682: Level: intermediate
684: .seealso TaoVec::PointwiseMultiply()
685: @*/
686: int TaoVec::PointwiseDivide( TaoVec* vv , TaoVec* ww){
687: int info;
688: TaoInt nn1,nn2,nn3,i;
689: TaoScalar *tptr1,*tptr2,*tptr3;
691: TaoFunctionBegin;
692: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
693: if (this!=vv){
694: info = vv->GetArray(&tptr2,&nn2);CHKERRQ(info);
695: } else {
696: tptr2=tptr1; nn2=nn1;
697: }
698: if (this!=ww && vv!=ww){
699: info = ww->GetArray(&tptr3,&nn3);CHKERRQ(info);
700: } else if (vv==ww){
701: tptr3=tptr2; nn3=nn2;
702: } else {
703: tptr3=tptr1; nn3=nn1;
704: }
705: if (nn1!=nn2 || nn2!=nn3) {TaoFunctionReturn(1);}
706: for (i=0;i<nn1;i++) tptr1[i]=tptr2[i] / tptr3[i];
707: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
708: if (this!=vv){
709: info = vv->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
710: }
711: if (this!=ww && vv!=ww){
712: info = ww->RestoreArray(&tptr3,&nn3);CHKERRQ(info);
713: }
714: TaoFunctionReturn(0);
715: }
719: /*@C
720: Median - Computes the componentwise median of three vectors
721: and stores the result in this vector.
723: Input Parameters:
724: . vv, ww, xx - the three vectors
726: Level: intermediate
728: @*/
729: int TaoVec::Median( TaoVec* vv, TaoVec* ww, TaoVec* xx){
730: int info;
731: TaoInt nn1,nn2,nn3,nn4,i;
732: TaoScalar *tptr1,*tptr2,*tptr3,*tptr4;
734: TaoFunctionBegin;
735: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
736: if (this!=vv){
737: info = vv->GetArray(&tptr2,&nn2);CHKERRQ(info);
738: } else {
739: tptr2=tptr1; nn2=nn1;
740: }
741: if (this!=ww){
742: info = ww->GetArray(&tptr3,&nn3);CHKERRQ(info);
743: } else {
744: tptr3=tptr1; nn3=nn1;
745: }
746: if (this!=xx){
747: info = xx->GetArray(&tptr4,&nn4);CHKERRQ(info);
748: } else {
749: tptr4=tptr1; nn4=nn1;
750: }
752: if (nn1!=nn2 || nn2!=nn3 || nn3!=nn4) {TaoFunctionReturn(1);}
753: for (i=0;i<nn1;i++){
754: if (tptr2[i]<=tptr3[i] && tptr3[i] <= tptr4[i]){
755: tptr1[i]=tptr3[i];
756: } else if (tptr4[i]<=tptr3[i] && tptr3[i] <= tptr2[i]){
757: tptr1[i]=tptr3[i];
758: } else if (tptr3[i]<=tptr2[i] && tptr2[i] <= tptr4[i]){
759: tptr1[i]=tptr2[i];
760: } else if (tptr4[i]<=tptr2[i] && tptr2[i] <= tptr3[i]){
761: tptr1[i]=tptr2[i];
762: } else {
763: tptr1[i]=tptr4[i];
764: }
765: }
766: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
767: if (this!=vv){
768: info = vv->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
769: }
770: if (this!=ww){
771: info = ww->RestoreArray(&tptr3,&nn3);CHKERRQ(info);
772: }
773: if (this!=xx){
774: info = xx->RestoreArray(&tptr4,&nn4);CHKERRQ(info);
775: }
776: TaoFunctionReturn(0);
777: }
781: /*@C
782: PointwiseMinimum - Computes the componentwise minimum of two vectors
783: and stores the result in this vector.
785: Input Parameters:
786: . vv, ww - the two vectors
788: Level: intermediate
790: .seealso TaoVec::PointwiseMaximum()
791: @*/
792: int TaoVec::PointwiseMinimum( TaoVec* vv, TaoVec* ww){
793: int info;
794: TaoInt nn1,nn2,nn3,i;
795: TaoScalar *tptr1,*tptr2,*tptr3;
797: TaoFunctionBegin;
798: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
799: if (vv!=this){
800: info = vv->GetArray(&tptr2,&nn2);CHKERRQ(info);
801: } else {
802: tptr2=tptr1; nn2=nn1;
803: }
804: if (ww!=this){
805: info = ww->GetArray(&tptr3,&nn3);CHKERRQ(info);
806: } else {
807: tptr3=tptr1; nn3=nn1;
808: }
810: if (nn1!=nn2 || nn2!=nn3) {TaoFunctionReturn(1);}
811: for (i=0;i<nn1;i++) tptr1[i] = TaoMin( tptr2[i] , tptr3[i]);
812: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
813: if (vv!=this){
814: info = vv->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
815: }
816: if (ww!=this){
817: info = ww->RestoreArray(&tptr3,&nn3);CHKERRQ(info);
818: }
819: TaoFunctionReturn(0);
820: }
824: /*@C
825: PointwiseMaximum - Computes the componentwise minimum of two vectors
826: and stores the result in this vector.
828: Input Parameters:
829: . vv, ww - the two vectors
831: Level: intermediate
833: .seealso TaoVec::PointwiseMinimum()
834: @*/
835: int TaoVec::PointwiseMaximum( TaoVec* vv, TaoVec* ww){
836: int info;
837: TaoInt nn1,nn2,nn3,i;
838: TaoScalar *tptr1,*tptr2,*tptr3;
840: TaoFunctionBegin;
841: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
842: if (vv!=this){
843: info = vv->GetArray(&tptr2,&nn2);CHKERRQ(info);
844: } else {
845: tptr2=tptr1; nn2=nn1;
846: }
847: if (ww!=this){
848: info = ww->GetArray(&tptr3,&nn3);CHKERRQ(info);
849: } else {
850: tptr3=tptr1; nn3=nn1;
851: }
852: if (nn1!=nn2 || nn2!=nn3) {TaoFunctionReturn(1);}
853: for (i=0;i<nn1;i++) tptr1[i] = TaoMax( tptr2[i] , tptr3[i]);
854: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
855: if (vv!=this){
856: info = vv->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
857: }
858: if (ww!=this){
859: info = ww->RestoreArray(&tptr3,&nn3);CHKERRQ(info);
860: }
861: TaoFunctionReturn(0);
862: }
866: /*@C
867: Waxpby - Sums two scaled vectors and stores the result in this vector. (this=alpha*xx+beta*yy)
869: Input Parameters:
870: + a - the multiple of the first vector
871: . xx - the first vector
872: . b - the multiple of the second vector
873: - yy - the second vector
875: Level: intermediate
877: .seealso TaoVec::Axpy(), TaoVec::Axpby(), TaoVec::Aypx();
878: @*/
879: int TaoVec::Waxpby ( double a, TaoVec* xx, double b, TaoVec* yy){
880: int info;
881: TaoInt nn1,nn2,nn3,i;
882: TaoScalar *tptr1,*tptr2,*tptr3;
884: TaoFunctionBegin;
885: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
886: info = xx->GetArray(&tptr2,&nn2);CHKERRQ(info);
887: info = yy->GetArray(&tptr3,&nn3);CHKERRQ(info);
888: if (nn1!=nn2 || nn2!=nn3) {TaoFunctionReturn(1);}
889: for (i=0;i<nn1;i++){ tptr1[i] = a * tptr2[i] + b * tptr3[i]; }
890: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
891: info = xx->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
892: info = yy->RestoreArray(&tptr3,&nn3);CHKERRQ(info);
893: TaoFunctionReturn(0);
894: }
898: /*@C
899: AbsoluteValue - Sets each element of this vector equal to its magnitude.
901: Input Parameters: none
903: Level: intermediate
904: @*/
905: int TaoVec::AbsoluteValue(){
906: int info;
907: TaoInt nn,i;
908: TaoScalar *v;
910: TaoFunctionBegin;
911: info = this->GetArray(&v,&nn);CHKERRQ(info);
912: for (i=0;i<nn;i++){
913: v[i]= TaoAbsScalar(v[i]);
914: }
915: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
916: TaoFunctionReturn(0);
917: }
921: /*@C
922: MinElement - Finds the smallest element of this vector.
924: Output Parameter:
925: . val - the smallest value in the vector
927: Level: intermediate
928: @*/
929: int TaoVec::MinElement(double *val){
930: int info;
931: TaoInt nn,i;
932: TaoScalar dd=TAO_INFINITY,*v;
934: TaoFunctionBegin;
935: info = this->GetArray(&v,&nn);CHKERRQ(info);
936: for (i=0;i<nn;i++){
937: if (v[i]<dd) dd=v[i];
938: }
939: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
940: *val = dd;
941: TaoFunctionReturn(0);
942: }
946: /*@C
947: Sign - Replace each element of this vector with -1, 0, or 1, depending on its sign.
949: Level: intermediate
950: @*/
951: int TaoVec::Sign(){
952: int info;
953: TaoInt nn,i;
954: TaoScalar *v;
956: TaoFunctionBegin;
957: info = this->GetArray(&v,&nn);CHKERRQ(info);
958: for (i=0;i<nn;i++){
959: if (v[i]<0){
960: v[i]=-1;
961: } else if (v[i]>0){
962: v[i]=1.0;
963: } else {
964: v[i]=0.0;
965: }
966: }
967: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
968: TaoFunctionReturn(0);
969: }
973: /*@C
974: SetReducedVec - Sets the reduced space of the vector that this
975: vector should represent. The index set describes which
976: elements of the vector should be used. This routine also
977: copies the appropriate elements from the full space vector
978: into the reduced space vector.
980: Input Parameters:
981: + vv - a vector
982: - ss - an index set
984: Level: advanced
986: .seealso TaoVec::CreateReducedVec(), TaoVec::ReducedCopyFromFull(), TaoVec::ReducedXPY(), TaoSelectSubset()
987: @*/
988: int TaoVec::SetReducedVec(TaoVec* vv,TaoIndexSet* ss){
989: TaoFunctionBegin;
990: SETERRQ(56,"Operation not defined");
991: /* TaoFunctionReturn(1); */
992: }
996: /*@C
997: ReducedCopyFromFull - Copies the appropriate elements of the vector into
998: this reduced vector.
999:
1000: Input Parameters:
1001: + ss - an index set
1002: - vv - a full vector
1004: Level: advanced
1006: .seealso TaoVec::CreateReducedVec(), TaoVec::ReducedXPY()
1007: @*/
1008: int TaoVec::ReducedCopyFromFull(TaoVec* vv, TaoIndexSet* ss){
1009: TaoFunctionBegin;
1010: SETERRQ(56,"Operation not defined");
1011: /* TaoFunctionReturn(1); */
1012: }
1016: /*@C
1017: ReducedXPY - Adds a reduced vector to the appropriate elements of this vector.
1018:
1019: Input Parameters:
1020: + vv - the reduced vector
1021: - ss - the index set identifying which elements of this vector should be supplemented
1023: Level: advanced
1025: .seealso TaoVec::CreateReducedVec(), TaoVec::ReducedCopyFromFull()
1026: @*/
1027: int TaoVec::ReducedXPY(TaoVec* vv, TaoIndexSet* ss){
1028: TaoFunctionBegin;
1029: SETERRQ(56,"Operation not defined");
1030: /* TaoFunctionReturn(1); */
1031: }
1035: /*@C
1036: StepMax - Calculates the largest multiple of a vector that can be added to
1037: this vector while keeping each element of this vector nonnegative.
1038:
1039: Input Parameters:
1040: . vv - the step direction
1042: Input Parameters:
1043: . step1 - the maximum stepsize
1045: Note:
1046: This vector should contain all positive elements.
1048: Note:
1049: If there is no maximum steplength, the output scalar may be set
1050: to TAO_INFINITY.
1052: Level: advanced
1053: @*/
1054: int TaoVec::StepMax(TaoVec* vv,double *step1){
1055: int info;
1056: TaoInt nn1,nn2,i;
1057: TaoScalar *xx,*dx;
1058: double stepmax1=TAO_INFINITY;
1060: TaoFunctionBegin;
1061: info = this->GetArray(&xx,&nn1);CHKERRQ(info);
1062: if (this!=vv){
1063: info = vv->GetArray(&dx,&nn2);CHKERRQ(info);
1064: if (nn1!=nn2) {TaoFunctionReturn(1);}
1065: for (i=0;i<nn1;i++){
1066: if (xx[i] < 0){
1067: TaoFunctionReturn(1);
1068: } else if (dx[i]<0){
1069: stepmax1=TaoMin(stepmax1,-xx[i]/dx[i]);
1070: }
1071: }
1072: info = vv->RestoreArray(&dx,&nn2);CHKERRQ(info);
1073: *step1=stepmax1;
1074: } else {
1075: *step1=1.0;
1076: }
1077: info = this->RestoreArray(&xx,&nn1);CHKERRQ(info);
1078: TaoFunctionReturn(0);
1079: }
1081: /*@C
1082: StepBoundInfo - Calculates the largest multiple of a vector that can be added to
1083: this vector while keeping each element of this vector nonnegative.
1084:
1085: Input Parameters:
1086: + txl - the lower bounds on this vector
1087: . txu - the upper bounds on this vector
1088: - tdx - the step direction for this vector
1090: Output Parameters:
1091: + boundmin - the step to closest bound i.e min(a1, ..., an);
1092: . wolfemin - the step to closest bound not equal i.e min(b1, ..., bn);
1093: - boundmax - the step to farthest bound i.e. max(c1, ..., cn);
1095: Where:
1096: if tdx[i] > 0; ai = (txu[i] - this[i])/tdx[i] ; bi=ai, ci=ai;
1097: if tdx[i] < 0; ai = (txl[i] - this[i])/tdx[i] ; bi=ai, ci=ai
1098: if tdx[i] == 0 && txl[i] < x[i] < txu[i] ; ai=TAO_INFINITY, bi=ai, ci=ai;
1099: if tdx[i] == 0 && (txl[i] == x[i] || txu[i] == x[i]) ; ai= 0, bi=TAO_INFINITY, ci=0;
1100:
1101: Note:
1102: If there is no maximum steplength, the output scalar may be set
1103: to TAO_INFINITY.
1105: Level: advanced
1106: @*/
1107: int TaoVec::StepBoundInfo(TaoVec* XL ,TaoVec* XU, TaoVec*S, double *bmin1,double *bmin2, double *bmax){
1108: TaoFunctionBegin;
1109: SETERRQ(56,"Operation not defined");
1110: /* TaoFunctionReturn(1); */
1111: }
1115: /*@C
1116: View - Views the contents of the vector.
1117:
1118: Input Parameters: none
1119:
1120: Level: intermediate
1121: @*/
1122: int TaoVec::View(){
1123: TaoFunctionBegin;
1124: SETERRQ(56,"Operation not defined");
1125: /* TaoFunctionReturn(1); */
1126: }
1130: /*@C
1131: BoundGradientProjection - Projects this vector according to this definition.
1132: If XX[i]==XXL[i], then this[i] = min(G[i],0);
1133: If XX[i]==XXU[i], then this[i] = max(G[i],0);
1134: else this[i] = G[i];
1136: Input Parameters:
1137: . GG,XX,XXL,XXU - the vectors.
1139: Level: advanced
1140: @*/
1141: int TaoVec::BoundGradientProjection(TaoVec* gg,TaoVec* xxll,TaoVec* xx, TaoVec* xxuu){
1142: int info;
1143: TaoInt nn1,nn2,nn3,nn4,nn5,i;
1144: TaoScalar *xptr,*xlptr,*xuptr,*gptr,*gpptr;
1146: TaoFunctionBegin;
1147: info = this->GetArray(&gpptr,&nn1);CHKERRQ(info);
1148: if (this != gg){
1149: info = gg->GetArray(&gptr,&nn2);CHKERRQ(info);
1150: } else {
1151: gptr=gpptr; nn2=nn1;
1152: }
1153: info = xxll->GetArray(&xlptr,&nn3);CHKERRQ(info);
1154: info = xx->GetArray(&xptr,&nn4);CHKERRQ(info);
1155: info = xxuu->GetArray(&xuptr,&nn5);CHKERRQ(info);
1156: if (nn1!=nn2 || nn2!=nn3 || nn3!=nn4 || nn4!=nn5) {TaoFunctionReturn(1);}
1158: for (i=0; i<nn1; i++){
1160: gpptr[i] = gptr[i];
1161: if (gpptr[i]>0 && xptr[i]<=xlptr[i]){
1162: gpptr[i] = 0;
1163: } else if (gpptr[i]<0 && xptr[i]>=xuptr[i]){
1164: gpptr[i] = 0;
1165: }
1166: }
1167: info = this->RestoreArray(&gpptr,&nn1);CHKERRQ(info);
1168: if (this!=gg){
1169: info = gg->RestoreArray(&gptr,&nn2);CHKERRQ(info);
1170: }
1171: info = xxll->RestoreArray(&xlptr,&nn3);CHKERRQ(info);
1172: info = xx->RestoreArray(&xptr,&nn4);CHKERRQ(info);
1173: info = xxuu->RestoreArray(&xuptr,&nn5);CHKERRQ(info);
1175: TaoFunctionReturn(0);
1176: }
1180: /*@C
1181: GetArray - Sets a pointer to the first element in the vector array.
1183: Output Parameters:
1184: + dptr - pointer an the array of numbers
1185: - n - the length of the array
1187: Note:
1188: This operation may not be defined the same for all vector types.
1190: Level: intermediate
1192: .seealso TaoVec::RestoreArray()
1193: @*/
1194: int TaoVec::GetArray(TaoScalar **dptr, TaoInt *n){
1195: TaoFunctionBegin;
1196: SETERRQ(56,"Operation not defined");
1197: /* TaoFunctionReturn(1); */
1198: }
1202: /*@C
1203: RestoreArray - Returns a pointer to the first element in the vector array.
1205: Input Parameters:
1206: + dptr - pointer an the array of numbers obtained from TaoVec::GetArray
1207: - n - the length of the array
1209: Note:
1210: This routine is not used within TAO solvers. Rather, it is used to
1211: implement other routines.
1213: Level: intermediate
1215: .seealso TaoVec::GetArray()
1216: @*/
1217: int TaoVec::RestoreArray(TaoScalar **dptr, TaoInt *n){
1218: TaoFunctionBegin;
1219: *dptr=NULL;
1220: *n=0;
1221: TaoFunctionReturn(0);
1222: }
1226: /*@C
1227: CreateIndexSet - Creates an index set that may be used to describe sets of
1228: elements of this vector.
1230: Output Parameters:
1231: . ss - a new index set
1233: Level: advanced
1234: @*/
1235: int TaoVec::CreateIndexSet(TaoIndexSet **ss){
1236: TaoFunctionBegin;
1237: SETERRQ(56,"Operation not defined");
1238: /* TaoFunctionReturn(1); */
1239: }
1241: inline static double fischer(double a, double b)
1242: {
1243: // Method suggested by Bob Vanderbei
1244: if (a + b <= 0) {
1245: return sqrt(a*a + b*b) - (a + b);
1246: }
1247: return -2.0*a*b / (sqrt(a*a + b*b) + (a + b));
1248: }
1252: /*@C
1253: Fischer - Evaluates the Fischer-Burmeister function for complementarity
1254: problems.
1255:
1256: Input Parameters:
1257: + xx - current point
1258: . ff - function evaluated at x
1259: . ll - lower bounds
1260: - uu - upper bounds
1262: Notes:
1263: The Fischer-Burmeister function is defined as
1264: $ phi(a,b) := sqrt(a*a + b*b) - a - b
1265: and is used reformulate a complementarity problem as a semismooth
1266: system of equations.
1268: The result of this function is done by cases:
1269: + l[i] == -infinity, u[i] == infinity -- fb[i] = -f[i]
1270: . l[i] == -infinity, u[i] finite -- fb[i] = phi(u[i]-x[i], -f[i])
1271: . l[i] finite, u[i] == infinity -- fb[i] = phi(x[i]-l[i], f[i])
1272: . l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u]))
1273: - otherwise l[i] == u[i] -- fb[i] = l[i] - x[i]
1275: Level: advanced
1277: .seealso TaoMat::D_Fischer()
1278: @*/
1279: int TaoVec::Fischer(TaoVec* xx, TaoVec* ff, TaoVec* ll, TaoVec* uu)
1280: {
1281: int info;
1282: TaoInt nn1,nn2,nn3,nn4,nn5,i;
1283: TaoScalar *v,*x,*f,*l,*u;
1285: TaoFunctionBegin;
1286: info = this->GetArray(&v,&nn1); CHKERRQ(info);
1287: info = xx->GetArray(&x,&nn2); CHKERRQ(info);
1288: info = ff->GetArray(&f,&nn3); CHKERRQ(info);
1289: info = ll->GetArray(&l,&nn4); CHKERRQ(info);
1290: info = uu->GetArray(&u,&nn5); CHKERRQ(info);
1291: if (nn1!=nn2 || nn2!=nn3 || nn3!=nn4 || nn4!=nn5) {
1292: TaoFunctionReturn(1);
1293: }
1295: for (i=0;i<nn1;i++) {
1297: if ((l[i] <= -TAO_INFINITY) && (u[i] >= TAO_INFINITY)) {
1298: v[i] = -f[i];
1299: }
1300: else if (l[i] <= -TAO_INFINITY) {
1301: v[i] = -fischer(u[i] - x[i], -f[i]);
1302: }
1303: else if (u[i] >= TAO_INFINITY) {
1304: v[i] = fischer(x[i] - l[i], f[i]);
1305: }
1306: else if (l[i] == u[i]) {
1307: v[i] = l[i] - x[i];
1308: }
1309: else {
1310: v[i] = fischer(u[i] - x[i], -f[i]);
1311: v[i] = fischer(x[i] - l[i], v[i]);
1312: }
1314: }
1316: info = this->RestoreArray(&v,&nn1);CHKERRQ(info);
1317: info = xx->RestoreArray(&x,&nn2);CHKERRQ(info);
1318: info = ff->RestoreArray(&f,&nn3);CHKERRQ(info);
1319: info = ll->RestoreArray(&l,&nn4);CHKERRQ(info);
1320: info = uu->RestoreArray(&u,&nn5);CHKERRQ(info);
1322: TaoFunctionReturn(0);
1323: }
1325: inline static double sfischer(double a, double b, double c)
1326: {
1327: // Method suggested by Bob Vanderbei
1328: if (a + b <= 0) {
1329: return sqrt(a*a + b*b + 2.0*c*c) - (a + b);
1330: }
1331: return 2.0*(c*c - a*b) / (sqrt(a*a + b*b + 2.0*c*c) + (a + b));
1332: }
1336: /*@C
1337: SFischer - Evaluates the Smoothed Fischer-Burmeister function for
1338: complementarity problems.
1340: Input Parameters:
1341: + xx - current point
1342: . ff - function evaluated at x
1343: . ll - lower bounds
1344: . uu - upper bounds
1345: - mu - smoothing parameter
1347: Notes:
1348: The Smoothed Fischer-Burmeister function is defined as
1349: $ phi(a,b) := sqrt(a*a + b*b + 2*mu*mu) - a - b
1350: and is used reformulate a complementarity problem as a semismooth
1351: system of equations.
1353: The result of this function is done by cases:
1354: + l[i] == -infinity, u[i] == infinity -- fb[i] = -f[i] - 2*mu*x[i]
1355: . l[i] == -infinity, u[i] finite -- fb[i] = phi(u[i]-x[i], -f[i], mu)
1356: . l[i] finite, u[i] == infinity -- fb[i] = phi(x[i]-l[i], f[i], mu)
1357: . l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u], mu), mu)
1358: - otherwise l[i] == u[i] -- fb[i] = l[i] - x[i]
1360: Level: advanced
1362: .seealso TaoMat::SD_Fischer()
1363: @*/
1364: int TaoVec::SFischer(TaoVec* xx, TaoVec* ff, TaoVec* ll, TaoVec* uu, double mu)
1365: {
1367: int info;
1368: TaoInt nn1, nn2, nn3, nn4, nn5,i;
1369: TaoScalar *v, *x, *f, *l, *u;
1371: TaoFunctionBegin;
1373: if ((mu >= -TAO_EPSILON) && (mu <= TAO_EPSILON)) {
1374: Fischer(xx, ff, ll, uu);
1375: }
1376: else {
1377: info = this->GetArray(&v, &nn1); CHKERRQ(info);
1378: info = xx->GetArray(&x, &nn2); CHKERRQ(info);
1379: info = ff->GetArray(&f, &nn3); CHKERRQ(info);
1380: info = ll->GetArray(&l, &nn4); CHKERRQ(info);
1381: info = uu->GetArray(&u, &nn5); CHKERRQ(info);
1383: if (nn1!=nn2 || nn2!=nn3 || nn3!=nn4 || nn4!=nn5) {
1384: TaoFunctionReturn(1);
1385: }
1387: for (i = 0; i < nn1; ++i) {
1388: if ((l[i] <= -TAO_INFINITY) && (u[i] >= TAO_INFINITY)) {
1389: v[i] = -f[i] - mu*x[i];
1390: }
1391: else if (l[i] <= -TAO_INFINITY) {
1392: v[i] = -sfischer(u[i] - x[i], -f[i], mu);
1393: }
1394: else if (u[i] >= TAO_INFINITY) {
1395: v[i] = sfischer(x[i] - l[i], f[i], mu);
1396: }
1397: else if (l[i] == u[i]) {
1398: v[i] = l[i] - x[i];
1399: }
1400: else {
1401: v[i] = sfischer(u[i] - x[i], -f[i], mu);
1402: v[i] = sfischer(x[i] - l[i], v[i], mu);
1403: }
1404: }
1406: info = this->RestoreArray(&v, &nn1); CHKERRQ(info);
1407: info = xx->RestoreArray(&x, &nn2); CHKERRQ(info);
1408: info = ff->RestoreArray(&f, &nn3); CHKERRQ(info);
1409: info = ll->RestoreArray(&l, &nn4); CHKERRQ(info);
1410: info = uu->RestoreArray(&u, &nn5); CHKERRQ(info);
1411: }
1412: TaoFunctionReturn(0);
1413: }