Actual source code: taovec.c
1: #include "tao_general.h" /*I "tao_solver.h" I*/
2: #include taovec.h
7: /*@C
8: TaoVecDestroy - Destroys the TaoVec object.
10: Input Parameter:
11: . vv - the vector
13: Level: beginner
14: @*/
15: int TaoVecDestroy( TaoVec* vv){
16: TaoFunctionBegin;
17: if (vv!=TAO_NULL && vv!=0){
18: delete vv;
19: }
20: vv=TAO_NULL;
21: TaoFunctionReturn(0);
22: }
26: /*@C
27: CloneVecs - Creates an array of pointers to new TaoVec objects. The new
28: objects have the same structure as this one.
30: Input Parameter:
31: . nn - number of new vectors
33: Output Parameter:
34: . tvs - pointer to array TaoVec pointers
36: .seealso TaoVec::Clone()
38: Level: intermediate
39: @*/
40: int TaoVec::CloneVecs(int nn, TaoVec***tvs){
41: int i,info;
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(int nn, TaoVec**ntv){
66: int i,info;
67: TaoFunctionBegin;
68: for (i=0;i<nn;i++){
69: info = TaoVecDestroy( ntv[i] );CHKERRQ(info);
70: }
71: delete[] ntv;
72: TaoFunctionReturn(0);
73: }
77: /*@C
78: Clone - Creates a new TaoVec object with the same structure as this one.
80: Input:
81: . vv - address of a pointer to a TaoVec
83: Output Parameter:
84: . vv - address of a pointer to new TaoVec object
86: .seealso TaoVec::CloneVecs()
88: Level: intermediate
89: @*/
90: int TaoVec::Clone( TaoVec* *vv ){
91: TaoFunctionBegin;
92: SETERRQ(56,"Operation not defined");
93: /* TaoFunctionReturn(1); */
94: }
98: /*@C
99: Compatible - Determines whether this vector belongs to the same space as another,
100: and operations such as inner product and sum are well defined.
102: Input Parameter:
103: . vv - TAO vector to which to the comparison is made
105: Output Value:
106: . flag - TAO_TRUE if the two vectors are Compatible and TAO_FALSE otherwise.
108: Level: advanced
110: .seealso TaoVec::GetDimension()
111: @*/
112: int TaoVec::Compatible(TaoVec* vv, TaoTruth *flag){
113: TaoFunctionBegin;
114: if (!flag){
115: TaoFunctionReturn(1);
116: }
117: *flag=TAO_FALSE;
118: TaoFunctionReturn(0);
119: }
123: /*@C
124: SetToConstant - Sets each element of this vector equal to the specified constant.
126: Input Parameter:
127: . c - a constant
129: Level: intermediate
131: .seealso TaoVec::Scale()
132: @*/
133: int TaoVec::SetToConstant( double c ){
134: int i,nn,info;
135: TaoScalar *tptr;
137: TaoFunctionBegin;
138: info = this->GetArray(&tptr,&nn);CHKERRQ(info);
139: for (i=0;i<nn;i++){ tptr[i]=c; }
140: info = this->RestoreArray(&tptr,&nn);CHKERRQ(info);
141: TaoFunctionReturn(0);
142: }
146: /*@C
147: SetToZero - Sets each element of this vector equal to zero.
149: Input Parameters: none
151: Level: intermediate
153: .seealso TaoVec::SetToConstant()
154: @*/
155: int TaoVec::SetToZero(){
156: int i,nn,info;
157: TaoScalar *tptr;
159: TaoFunctionBegin;
160: info = this->GetArray(&tptr,&nn);CHKERRQ(info);
161: for (i=0;i<nn;i++){ tptr[i]=0; }
162: info = this->RestoreArray(&tptr,&nn);CHKERRQ(info);
163: TaoFunctionReturn(0);
164: }
168: /*@C
169: CopyFrom - Copies the contents of one vector into this vector.
171: Input Parameter:
172: . vv - A TaoVec from which the contents will be copied.
174: Level: intermediate
176: .seealso TaoVec::Axpy()
177: @*/
178: int TaoVec::CopyFrom( TaoVec* vv ){
179: int i,nn1,nn2,info;
180: TaoScalar *tptr1,*tptr2;
182: TaoFunctionBegin;
183: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
184: if (vv!=this){
185: info = vv->GetArray(&tptr2,&nn2);CHKERRQ(info);
186: if (nn1!=nn2) {TaoFunctionReturn(1);}
187: for (i=0;i<nn1;i++){ tptr1[i]=tptr2[i]; }
188: info = vv->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
189: }
190: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
191: TaoFunctionReturn(0);
192: }
196: /*@C
197: ScaleCopyFrom - Copies the contents of one vector into this vector and scales it.
199: Input Parameter:
200: + a - the scalar
201: - vv - A TaoVec from which the contents will be copied.
203: Level: intermediate
205: .seealso TaoVec::Axpy(), TaoVec::Aypx()
206: @*/
207: int TaoVec::ScaleCopyFrom( double a, TaoVec* vv ){
208: int i,nn1,nn2,info;
209: TaoScalar *tptr1,*tptr2;
211: TaoFunctionBegin;
212: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
213: if (vv!=this){
214: info = vv->GetArray(&tptr2,&nn2);CHKERRQ(info);
215: if (nn1!=nn2) {TaoFunctionReturn(1);}
216: for (i=0;i<nn1;i++){ tptr1[i]=a*tptr2[i]; }
217: info = vv->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
218: } else {
219: this->Scale(a);
220: }
221: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
222: TaoFunctionReturn(0);
223: }
227: /*@C
228: NormInfinity - Computes the infinity norm of this vector.
230: Ouput Parameter:
231: . vnorm - the infinity norm of this vector
233: Level: intermediate
235: .seealso TaoVec::Norm1(), TaoVec::Norm2()
236: @*/
237: int TaoVec::NormInfinity(double *vnorm){
238: int i,nn,info;
239: TaoScalar dd=0, *v;
241: TaoFunctionBegin;
242: info = this->GetArray(&v,&nn);CHKERRQ(info);
243: for (i=0;i<nn;i++){
244: if (v[i]<0 && v[i]<-dd) dd=-v[i];
245: else if (v[i]>0 && v[i]>dd) dd=v[i];
246: }
247: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
248: *vnorm=dd;
249: TaoFunctionReturn(0);
250: }
254: /*@C
255: Norm1 - Computes the one-norm of this vector.
257: Ouput Parameter:
258: . vnorm - the one-norm of this vector
260: Level: intermediate
262: .seealso TaoVec::NormInfinity(), TaoVec::Norm2()
263: @*/
264: int TaoVec::Norm1(double *vnorm){
265: int i,nn,info;
266: TaoScalar dd=0, *v;
268: TaoFunctionBegin;
269: info = this->GetArray(&v,&nn);CHKERRQ(info);
270: for (i=0;i<nn;i++){
271: if (v[i]<0) dd-=v[i];
272: else dd+=v[i];
273: }
274: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
275: TaoFunctionReturn(0);
276: }
280: /*@C
281: Norm2 - Compute the two-norm of this vector.
283: Ouput Parameter:
284: . vnorm - the two-norm of this vector
286: Level: intermediate
288: .seealso TaoVec::Norm1(), TaoVec::NormInfinity(), TaoVec::Norm2squared()
289: @*/
290: int TaoVec::Norm2(double *vnorm){
291: int i,nn,info;
292: TaoScalar dd=0, *v;
294: TaoFunctionBegin;
295: info = this->GetArray(&v,&nn);CHKERRQ(info);
296: for (i=0;i<nn;i++) dd+=v[i]*v[i];
297: *vnorm=sqrt(dd);
298: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
299: TaoFunctionReturn(0);
300: }
304: /*@C
305: Norm2squared - Computes the square of the two norm of this vector.
307: Ouput Parameter:
308: . vnorm2 - the square of the two norm of this vector
310: Level: intermediate
312: .seealso TaoVec::Norm2()
313: @*/
314: int TaoVec::Norm2squared(double *vnorm2){
315: int i,nn,info;
316: TaoScalar dd=0, *v;
318: TaoFunctionBegin;
319: info = this->GetArray(&v,&nn);CHKERRQ(info);
320: for (i=0;i<nn;i++) dd+=v[i]*v[i];
321: *vnorm2=dd;
322: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
323: TaoFunctionReturn(0);
324: }
328: /*@C
329: Scale - Multiplies this vector by a scalar.
331: Input Parameter:
332: . alpha - the scalar
334: Level: intermediate
336: .seealso TaoVec::SetToConstant(), TaoVec::Aypx()
337: @*/
338: int TaoVec::Scale( double alpha ){
339: int i,nn,info;
340: TaoScalar *v;
342: TaoFunctionBegin;
343: info = this->GetArray(&v,&nn);CHKERRQ(info);
344: for (i=0;i<nn;i++) v[i]*=alpha;
345: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
346: TaoFunctionReturn(0);
347: }
351: /*@C
352: Axpy - Adds a scalar multiple of a vector to this vector. (this += alpha * xx)
354: Input Parameter:
355: + alpha - the scalar
356: - xx - the vector
358: Level: intermediate
360: .seealso TaoVec::CopyFrom(), TaoVec::Aypx()
361: @*/
362: int TaoVec::Axpy( double alpha, TaoVec* xx ){
363: int i,nn1,nn2,info;
364: TaoScalar *tptr1,*tptr2;
366: TaoFunctionBegin;
367: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
368: info = xx->GetArray(&tptr2,&nn2);CHKERRQ(info);
369: if (nn1!=nn2) {TaoFunctionReturn(1);}
370: for (i=0;i<nn1;i++){ tptr1[i]+= alpha * tptr2[i]; }
371: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
372: info = xx->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
373: TaoFunctionReturn(0);
374: }
378: /*@C
379: Axpby - Adds a scalar multiple of a vector to a multiple of this vector. (this=alpha*xx + beta*this)
381: Input Parameter:
382: + alpha - the scalar of tx
383: . xx - the vector
384: - beta - the scalar multiple of this vector
386: Level: intermediate
388: .seealso TaoVec::Axpy(), TaoVec::Aypx()
389: @*/
390: int TaoVec::Axpby( double alpha, TaoVec* xx, double beta ){
391: int i,nn1,nn2,info;
392: TaoScalar *tptr1,*tptr2;
394: TaoFunctionBegin;
395: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
396: info = xx->GetArray(&tptr2,&nn2);CHKERRQ(info);
397: if (nn1!=nn2) {TaoFunctionReturn(1);}
398: for (i=0;i<nn1;i++){ tptr1[i] = beta * tptr1[i] + alpha * tptr2[i]; }
399: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
400: info = xx->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
401: TaoFunctionReturn(0);
402: }
406: /*@C
407: Aypx - Adds a vector to a scalar multiple of this vector. (this=alpha*this+xx)
409: Input Parameter:
410: + alpha - the scalar
411: - xx - the vector
413: Level: intermediate
415: .seealso TaoVec::Scale(), TaoVec::Axpy()
416: @*/
417: int TaoVec::Aypx( double alpha, TaoVec* xx ){
418: int i,nn1,nn2,info;
419: TaoScalar *tptr1,*tptr2;
421: TaoFunctionBegin;
422: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
423: info = xx->GetArray(&tptr2,&nn2);CHKERRQ(info);
424: if (nn1!=nn2) {TaoFunctionReturn(1);}
425: for (i=0;i<nn1;i++){ tptr1[i] = alpha * tptr1[i] + tptr2[i]; }
426: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
427: info = xx->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
428: TaoFunctionReturn(0);
429: }
430:
433: /*@C
434: AddConstant - Adds a constant to each element of this vector.
436: Input Parameter:
437: . alpha - the scalar
439: Level: intermediate
441: .seealso TaoVec::SetToConstant(), TaoVec::Axpy()
442: @*/
443: int TaoVec::AddConstant( double alpha ){
444: int i,nn,info;
445: TaoScalar *v;
447: TaoFunctionBegin;
448: info = this->GetArray(&v,&nn);CHKERRQ(info);
449: for (i=0;i<nn;i++) v[i]+=alpha;
450: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
451: TaoFunctionReturn(0);
452: }
456: /*@C
457: Dot - Computes the inner product of this vector with another vector.
459: Input Parameter:
460: . vv - another TaoVec object
462: Output Parameter:
463: . vDotv - the inner product of the two vectors
465: Level: intermediate
467: .seealso TaoVec::Norm()
468: @*/
469: int TaoVec::Dot( TaoVec* vv, double *vDotv ){
470: int i,nn1,nn2,info;
471: TaoScalar dd=0,*tptr1,*tptr2;
473: TaoFunctionBegin;
474: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
475: info = vv->GetArray(&tptr2,&nn2);CHKERRQ(info);
476: if (nn1!=nn2) {TaoFunctionReturn(1);}
477: for (i=0;i<nn1;i++) dd+=tptr1[i]*tptr2[i];
478: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
479: info = vv->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
480: *vDotv=dd;
481: TaoFunctionReturn(0);
482: }
486: /*@C
487: Negate - Multiplies the elements of this vector by negative one.
489: Input Parameters: none
491: Level: intermediate
493: .seealso TaoVec::Scale()
494: @*/
495: int TaoVec::Negate(){
496: int i,nn,info;
497: TaoScalar *v;
499: TaoFunctionBegin;
500: info = this->GetArray(&v,&nn);CHKERRQ(info);
501: for (i=0;i<nn;i++) v[i]=-v[i];
502: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
503: TaoFunctionReturn(0);
504: }
508: /*@C
509: Reciprocal - Sets each element of this vector to its Reciprocal.
511: Input Parameters: none
513: Level: intermediate
515: .seealso TaoVec::PointwiseDivide()
516: @*/
517: int TaoVec::Reciprocal(){
518: int i,nn,info;
519: TaoScalar *v;
521: TaoFunctionBegin;
522: info = this->GetArray(&v,&nn);CHKERRQ(info);
523: for (i=0;i<nn;i++){
524: if (v[i]!=0) v[i]= 1.0/v[i];
525: else v[i]=TAO_INFINITY;
526: }
527: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
528: TaoFunctionReturn(0);
529: }
533: /*@C
534: GetDimension - Gets the dimension of the vector space where this vector belongs.
536: Output Parameter:
537: . n - the dimension of the vector space
539: Level: intermediate
541: .seealso TaoVec::Compatible()
542: @*/
543: int TaoVec::GetDimension(int *n){
544: TaoFunctionBegin;
545: SETERRQ(56,"Operation not defined");
546: /* TaoFunctionReturn(1); */
547: }
551: /*@C
552: PointwiseMultiply - Computes the componentwise multiplication of two vectors
553: and stores the result in this vector.
555: Input Parameters:
556: . vv, ww - the two vectors
558: Level: intermediate
560: .seealso TaoVec::PointwiseDivide()
561: @*/
562: int TaoVec::PointwiseMultiply( TaoVec* vv, TaoVec* ww ){
563: int i,nn1,nn2,nn3,info;
564: TaoScalar *tptr1,*tptr2,*tptr3;
566: TaoFunctionBegin;
567: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
568: if (this!=vv){
569: info = vv->GetArray(&tptr2,&nn2);CHKERRQ(info);
570: } else {
571: tptr2=tptr1; nn2=nn1;
572: }
573: if (this!=ww && vv!=ww){
574: info = ww->GetArray(&tptr3,&nn3);CHKERRQ(info);
575: } else if (vv==ww){
576: tptr3=tptr2; nn3=nn2;
577: } else {
578: tptr3=tptr1; nn3=nn1;
579: }
580: if (nn1!=nn2 || nn2!=nn3) {TaoFunctionReturn(1);}
581: for (i=0;i<nn1;i++) tptr1[i]=tptr2[i] * tptr3[i];
582: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
583: if (this!=vv){
584: info = vv->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
585: }
586: if (this!=ww && vv!=ww){
587: info = ww->RestoreArray(&tptr3,&nn3);CHKERRQ(info);
588: }
590: TaoFunctionReturn(0);
591: }
595: /*@C
596: PointwiseDivide - Computes the componentwise division of two vectors
597: and stores the result in this vector.
599: Input Parameters:
600: + vv - the vector of numerators
601: - ww - the vector of denominators
603: Level: intermediate
605: .seealso TaoVec::PointwiseMultiply()
606: @*/
607: int TaoVec::PointwiseDivide( TaoVec* vv , TaoVec* ww){
608: int i,nn1,nn2,nn3,info;
609: TaoScalar *tptr1,*tptr2,*tptr3;
611: TaoFunctionBegin;
612: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
613: if (this!=vv){
614: info = vv->GetArray(&tptr2,&nn2);CHKERRQ(info);
615: } else {
616: tptr2=tptr1; nn2=nn1;
617: }
618: if (this!=ww && vv!=ww){
619: info = ww->GetArray(&tptr3,&nn3);CHKERRQ(info);
620: } else if (vv==ww){
621: tptr3=tptr2; nn3=nn2;
622: } else {
623: tptr3=tptr1; nn3=nn1;
624: }
625: if (nn1!=nn2 || nn2!=nn3) {TaoFunctionReturn(1);}
626: for (i=0;i<nn1;i++) tptr1[i]=tptr2[i] / tptr3[i];
627: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
628: if (this!=vv){
629: info = vv->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
630: }
631: if (this!=ww && vv!=ww){
632: info = ww->RestoreArray(&tptr3,&nn3);CHKERRQ(info);
633: }
634: TaoFunctionReturn(0);
635: }
639: /*@C
640: Median - Computes the componentwise median of three vectors
641: and stores the result in this vector.
643: Input Parameters:
644: . vv, ww, xx - the three vectors
646: Level: intermediate
648: @*/
649: int TaoVec::Median( TaoVec* vv, TaoVec* ww, TaoVec* xx){
650: int i,nn1,nn2,nn3,nn4,info;
651: TaoScalar *tptr1,*tptr2,*tptr3,*tptr4;
653: TaoFunctionBegin;
654: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
655: if (this!=vv){
656: info = vv->GetArray(&tptr2,&nn2);CHKERRQ(info);
657: } else {
658: tptr2=tptr1; nn2=nn1;
659: }
660: if (this!=ww){
661: info = ww->GetArray(&tptr3,&nn3);CHKERRQ(info);
662: } else {
663: tptr3=tptr1; nn3=nn1;
664: }
665: if (this!=xx){
666: info = xx->GetArray(&tptr4,&nn4);CHKERRQ(info);
667: } else {
668: tptr4=tptr1; nn4=nn1;
669: }
671: if (nn1!=nn2 || nn2!=nn3 || nn3!=nn4) {TaoFunctionReturn(1);}
672: for (i=0;i<nn1;i++){
673: if (tptr2[i]<=tptr3[i] && tptr3[i] <= tptr4[i]){
674: tptr1[i]=tptr3[i];
675: } else if (tptr4[i]<=tptr3[i] && tptr3[i] <= tptr2[i]){
676: tptr1[i]=tptr3[i];
677: } else if (tptr3[i]<=tptr2[i] && tptr2[i] <= tptr4[i]){
678: tptr1[i]=tptr2[i];
679: } else if (tptr4[i]<=tptr2[i] && tptr2[i] <= tptr3[i]){
680: tptr1[i]=tptr2[i];
681: } else {
682: tptr1[i]=tptr4[i];
683: }
684: }
685: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
686: if (this!=vv){
687: info = vv->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
688: }
689: if (this!=ww){
690: info = ww->RestoreArray(&tptr3,&nn3);CHKERRQ(info);
691: }
692: if (this!=xx){
693: info = xx->RestoreArray(&tptr4,&nn4);CHKERRQ(info);
694: }
695: TaoFunctionReturn(0);
696: }
700: /*@C
701: PointwiseMinimum - Computes the componentwise minimum of two vectors
702: and stores the result in this vector.
704: Input Parameters:
705: . vv, ww - the two vectors
707: Level: intermediate
709: .seealso TaoVec::PointwiseMaximum()
710: @*/
711: int TaoVec::PointwiseMinimum( TaoVec* vv, TaoVec* ww){
712: int i,nn1,nn2,nn3,info;
713: TaoScalar *tptr1,*tptr2,*tptr3;
715: TaoFunctionBegin;
716: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
717: if (vv!=this){
718: info = vv->GetArray(&tptr2,&nn2);CHKERRQ(info);
719: } else {
720: tptr2=tptr1; nn2=nn1;
721: }
722: if (ww!=this){
723: info = ww->GetArray(&tptr3,&nn3);CHKERRQ(info);
724: } else {
725: tptr3=tptr1; nn3=nn1;
726: }
728: if (nn1!=nn2 || nn2!=nn3) {TaoFunctionReturn(1);}
729: for (i=0;i<nn1;i++) tptr1[i] = TaoMin( tptr2[i] , tptr3[i]);
730: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
731: if (vv!=this){
732: info = vv->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
733: }
734: if (ww!=this){
735: info = ww->RestoreArray(&tptr3,&nn3);CHKERRQ(info);
736: }
737: TaoFunctionReturn(0);
738: }
742: /*@C
743: PointwiseMaximum - Computes the componentwise minimum of two vectors
744: and stores the result in this vector.
746: Input Parameters:
747: . vv, ww - the two vectors
749: Level: intermediate
751: .seealso TaoVec::PointwiseMinimum()
752: @*/
753: int TaoVec::PointwiseMaximum( TaoVec* vv, TaoVec* ww){
754: int i,nn1,nn2,nn3,info;
755: TaoScalar *tptr1,*tptr2,*tptr3;
757: TaoFunctionBegin;
758: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
759: if (vv!=this){
760: info = vv->GetArray(&tptr2,&nn2);CHKERRQ(info);
761: } else {
762: tptr2=tptr1; nn2=nn1;
763: }
764: if (ww!=this){
765: info = ww->GetArray(&tptr3,&nn3);CHKERRQ(info);
766: } else {
767: tptr3=tptr1; nn3=nn1;
768: }
769: if (nn1!=nn2 || nn2!=nn3) {TaoFunctionReturn(1);}
770: for (i=0;i<nn1;i++) tptr1[i] = TaoMax( tptr2[i] , tptr3[i]);
771: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
772: if (vv!=this){
773: info = vv->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
774: }
775: if (ww!=this){
776: info = ww->RestoreArray(&tptr3,&nn3);CHKERRQ(info);
777: }
778: TaoFunctionReturn(0);
779: }
783: /*@C
784: Waxpby - Sums two scaled vectors and stores the result in this vector. (this=alpha*xx+beta*yy)
786: Input Parameters:
787: + a - the multiple of the first vector
788: . xx - the first vector
789: . b - the multiple of the second vector
790: - yy - the second vector
792: Level: intermediate
794: .seealso TaoVec::Axpy(), TaoVec::Axpby(), TaoVec::Aypx();
795: @*/
796: int TaoVec::Waxpby ( double a, TaoVec* xx, double b, TaoVec* yy){
797: int i,nn1,nn2,nn3,info;
798: TaoScalar *tptr1,*tptr2,*tptr3;
800: TaoFunctionBegin;
801: info = this->GetArray(&tptr1,&nn1);CHKERRQ(info);
802: info = xx->GetArray(&tptr2,&nn2);CHKERRQ(info);
803: info = yy->GetArray(&tptr3,&nn3);CHKERRQ(info);
804: if (nn1!=nn2 || nn2!=nn3) {TaoFunctionReturn(1);}
805: for (i=0;i<nn1;i++){ tptr1[i] = a * tptr2[i] + b * tptr3[i]; }
806: info = this->RestoreArray(&tptr1,&nn1);CHKERRQ(info);
807: info = xx->RestoreArray(&tptr2,&nn2);CHKERRQ(info);
808: info = yy->RestoreArray(&tptr3,&nn3);CHKERRQ(info);
809: TaoFunctionReturn(0);
810: }
814: /*@C
815: AbsoluteValue - Sets each element of this vector equal to its magnitude.
817: Input Parameters: none
819: Level: intermediate
820: @*/
821: int TaoVec::AbsoluteValue(){
822: int i,nn,info;
823: TaoScalar *v;
825: TaoFunctionBegin;
826: info = this->GetArray(&v,&nn);CHKERRQ(info);
827: for (i=0;i<nn;i++){
828: v[i]= TaoAbsScalar(v[i]);
829: }
830: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
831: TaoFunctionReturn(0);
832: }
836: /*@C
837: MinElement - Finds the smallest element of this vector.
839: Output Parameter:
840: . val - the smallest value in the vector
842: Level: intermediate
843: @*/
844: int TaoVec::MinElement(double *val){
845: int i,nn,info;
846: TaoScalar dd=TAO_INFINITY,*v;
848: TaoFunctionBegin;
849: info = this->GetArray(&v,&nn);CHKERRQ(info);
850: for (i=0;i<nn;i++){
851: if (v[i]<dd) dd=v[i];
852: }
853: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
854: *val = dd;
855: TaoFunctionReturn(0);
856: }
860: /*@C
861: Sign - Replace each element of this vector with -1, 0, or 1, depending on its sign.
863: Level: intermediate
864: @*/
865: int TaoVec::Sign(){
866: int i,nn,info;
867: TaoScalar *v;
869: TaoFunctionBegin;
870: info = this->GetArray(&v,&nn);CHKERRQ(info);
871: for (i=0;i<nn;i++){
872: if (v[i]<0){
873: v[i]=-1;
874: } else if (v[i]>0){
875: v[i]=1.0;
876: } else {
877: v[i]=0.0;
878: }
879: }
880: info = this->RestoreArray(&v,&nn);CHKERRQ(info);
881: TaoFunctionReturn(0);
882: }
886: /*@C
887: SetReducedVec - Sets the reduced space of the vector that this
888: vector should represent. The index set describes which
889: elements of the vector should be used. This routine also
890: copies the appropriate elements from the full space vector
891: into the reduced space vector.
893: Input Parameters:
894: + vv - a vector
895: - ss - an index set
897: Level: advanced
899: .seealso TaoVec::CreateReducedVec(), TaoVec::ReducedCopyFromFull(), TaoVec::ReducedXPY(), TaoSelectSubset()
900: @*/
901: int TaoVec::SetReducedVec(TaoVec* vv,TaoIndexSet* ss){
902: TaoFunctionBegin;
903: SETERRQ(56,"Operation not defined");
904: /* TaoFunctionReturn(1); */
905: }
909: /*@C
910: ReducedCopyFromFull - Copies the appropriate elements of the vector into
911: this reduced vector.
912:
913: Input Parameters:
914: + ss - an index set
915: - vv - a full vector
917: Level: advanced
919: .seealso TaoVec::CreateReducedVec(), TaoVec::ReducedXPY()
920: @*/
921: int TaoVec::ReducedCopyFromFull(TaoVec* vv, TaoIndexSet* ss){
922: TaoFunctionBegin;
923: SETERRQ(56,"Operation not defined");
924: /* TaoFunctionReturn(1); */
925: }
929: /*@C
930: ReducedXPY - Adds a reduced vector to the appropriate elements of this vector.
931:
932: Input Parameters:
933: + vv - the reduced vector
934: - ss - the index set identifying which elements of this vector should be supplemented
936: Level: advanced
938: .seealso TaoVec::CreateReducedVec(), TaoVec::ReducedCopyFromFull()
939: @*/
940: int TaoVec::ReducedXPY(TaoVec* vv, TaoIndexSet* ss){
941: TaoFunctionBegin;
942: SETERRQ(56,"Operation not defined");
943: /* TaoFunctionReturn(1); */
944: }
948: /*@C
949: StepMax - Calculates the largest multiple of a vector that can be added to
950: this vector while keeping each element of this vector nonnegative.
951:
952: Input Parameters:
953: . vv - the step direction
955: Input Parameters:
956: . step1 - the maximum stepsize
958: Note:
959: This vector should contain all positive elements.
961: Note:
962: If there is no maximum steplength, the output scalar may be set
963: to TAO_INFINITY.
965: Level: advanced
966: @*/
967: int TaoVec::StepMax(TaoVec* vv,double *step1){
968: int i,nn1,nn2,info;
969: TaoScalar *xx,*dx;
970: double stepmax1=TAO_INFINITY;
972: TaoFunctionBegin;
973: info = this->GetArray(&xx,&nn1);CHKERRQ(info);
974: if (this!=vv){
975: info = vv->GetArray(&dx,&nn2);CHKERRQ(info);
976: if (nn1!=nn2) {TaoFunctionReturn(1);}
977: for (i=0;i<nn1;i++){
978: if (xx[i] < 0){
979: TaoFunctionReturn(1);
980: } else if (dx[i]<0){
981: stepmax1=TaoMin(stepmax1,-xx[i]/dx[i]);
982: }
983: }
984: info = vv->RestoreArray(&dx,&nn2);CHKERRQ(info);
985: *step1=stepmax1;
986: } else {
987: *step1=1.0;
988: }
989: info = this->RestoreArray(&xx,&nn1);CHKERRQ(info);
990: TaoFunctionReturn(0);
991: }
993: /*@C
994: StepBoundInfo - Calculates the largest multiple of a vector that can be added to
995: this vector while keeping each element of this vector nonnegative.
996:
997: Input Parameters:
998: + txl - the lower bounds on this vector
999: . txu - the upper bounds on this vector
1000: - tdx - the step direction for this vector
1002: Output Parameters:
1003: + boundmin - the step to closest bound i.e min(a1, ..., an);
1004: . wolfemin - the step to closest bound not equal i.e min(b1, ..., bn);
1005: - boundmax - the step to farthest bound i.e. max(c1, ..., cn);
1007: Where:
1008: if tdx[i] > 0; ai = (txu[i] - this[i])/tdx[i] ; bi=ai, ci=ai;
1009: if tdx[i] < 0; ai = (txl[i] - this[i])/tdx[i] ; bi=ai, ci=ai
1010: if tdx[i] == 0 && txl[i] < x[i] < txu[i] ; ai=TAO_INFINITY, bi=ai, ci=ai;
1011: if tdx[i] == 0 && (txl[i] == x[i] || txu[i] == x[i]) ; ai= 0, bi=TAO_INFINITY, ci=0;
1012:
1013: Note:
1014: If there is no maximum steplength, the output scalar may be set
1015: to TAO_INFINITY.
1017: Level: advanced
1018: @*/
1019: int TaoVec::StepBoundInfo(TaoVec* XL ,TaoVec* XU, TaoVec*S, double *bmin1,double *bmin2, double *bmax){
1020: TaoFunctionBegin;
1021: SETERRQ(56,"Operation not defined");
1022: /* TaoFunctionReturn(1); */
1023: }
1027: /*@C
1028: View - Views the contents of the vector.
1029:
1030: Input Parameters: none
1031:
1032: Level: intermediate
1033: @*/
1034: int TaoVec::View(){
1035: TaoFunctionBegin;
1036: SETERRQ(56,"Operation not defined");
1037: /* TaoFunctionReturn(1); */
1038: }
1042: /*@C
1043: BoundGradientProjection - Projects this vector according to this definition.
1044: If XX[i]==XXL[i], then this[i] = min(G[i],0);
1045: If XX[i]==XXU[i], then this[i] = max(G[i],0);
1046: else this[i] = G[i];
1048: Input Parameters:
1049: . GG,XX,XXL,XXU - the vectors.
1051: Level: advanced
1052: @*/
1053: int TaoVec::BoundGradientProjection(TaoVec* gg,TaoVec* xxll,TaoVec* xx, TaoVec* xxuu){
1054: int i,nn1,nn2,nn3,nn4,nn5,info;
1055: TaoScalar *xptr,*xlptr,*xuptr,*gptr,*gpptr;
1057: TaoFunctionBegin;
1058: info = this->GetArray(&gpptr,&nn1);CHKERRQ(info);
1059: if (this != gg){
1060: info = gg->GetArray(&gptr,&nn2);CHKERRQ(info);
1061: } else {
1062: gptr=gpptr; nn2=nn1;
1063: }
1064: info = xxll->GetArray(&xlptr,&nn3);CHKERRQ(info);
1065: info = xx->GetArray(&xptr,&nn4);CHKERRQ(info);
1066: info = xxuu->GetArray(&xuptr,&nn5);CHKERRQ(info);
1067: if (nn1!=nn2 || nn2!=nn3 || nn3!=nn4 || nn4!=nn5) {TaoFunctionReturn(1);}
1069: for (i=0; i<nn1; i++){
1071: gpptr[i] = gptr[i];
1072: if (gpptr[i]>0 && xptr[i]<=xlptr[i]){
1073: gpptr[i] = 0;
1074: } else if (gpptr[i]<0 && xptr[i]>=xuptr[i]){
1075: gpptr[i] = 0;
1076: }
1077: }
1078: info = this->RestoreArray(&gpptr,&nn1);CHKERRQ(info);
1079: if (this!=gg){
1080: info = gg->RestoreArray(&gptr,&nn2);CHKERRQ(info);
1081: }
1082: info = xxll->RestoreArray(&xlptr,&nn3);CHKERRQ(info);
1083: info = xx->RestoreArray(&xptr,&nn4);CHKERRQ(info);
1084: info = xxuu->RestoreArray(&xuptr,&nn5);CHKERRQ(info);
1086: TaoFunctionReturn(0);
1087: }
1091: /*@C
1092: GetArray - Sets a pointer to the first element in the vector array.
1094: Output Parameters:
1095: + dptr - pointer an the array of numbers
1096: - n - the length of the array
1098: Note:
1099: This operation may not be defined the same for all vector types.
1101: Level: intermediate
1103: .seealso TaoVec::RestoreArray()
1104: @*/
1105: int TaoVec::GetArray(TaoScalar **dptr, int *n){
1106: TaoFunctionBegin;
1107: SETERRQ(56,"Operation not defined");
1108: /* TaoFunctionReturn(1); */
1109: }
1113: /*@C
1114: RestoreArray - Returns a pointer to the first element in the vector array.
1116: Input Parameters:
1117: + dptr - pointer an the array of numbers obtained from TaoVec::GetArray
1118: - n - the length of the array
1120: Note:
1121: This routine is not used within TAO solvers. Rather, it is used to
1122: implement other routines.
1124: Level: intermediate
1126: .seealso TaoVec::GetArray()
1127: @*/
1128: int TaoVec::RestoreArray(TaoScalar **dptr, int *n){
1129: TaoFunctionBegin;
1130: *dptr=NULL;
1131: *n=0;
1132: TaoFunctionReturn(0);
1133: }
1137: /*@C
1138: CreateIndexSet - Creates an index set that may be used to describe sets of
1139: elements of this vector.
1141: Output Parameters:
1142: . ss - a new index set
1144: Level: advanced
1145: @*/
1146: int TaoVec::CreateIndexSet(TaoIndexSet **ss){
1147: TaoFunctionBegin;
1148: SETERRQ(56,"Operation not defined");
1149: /* TaoFunctionReturn(1); */
1150: }
1152: inline static double fischer(double a, double b)
1153: {
1155: #ifdef TODD
1157: if (TaoAbsScalar(a) > TaoAbsScalar(b)) {
1158: return sqrt(a*a + b*b) - a - b;
1159: }
1160: return sqrt(a*a + b*b) - b - a;
1162: #else
1164: // Method suggested by Bob Vanderbei
1166: if (a + b <= 0) {
1167: return sqrt(a*a + b*b) - (a + b);
1168: }
1169: return -2.0*a*b / (sqrt(a*a + b*b) + (a + b));
1171: #endif
1173: }
1177: /*@C
1178: Fischer - Evaluates the Fischer-Burmeister function for complementarity
1179: problems.
1180:
1181: Input Parameters:
1182: + xx - current point
1183: . ff - function evaluated at x
1184: . ll - lower bounds
1185: - uu - upper bounds
1187: Notes:
1188: The Fischer-Burmeister function is defined as
1189: $ phi(a,b) := sqrt(a*a + b*b) - a - b
1190: and is used reformulate a complementarity problem as a semismooth
1191: system of equations.
1193: The result of this function is done by cases:
1194: + l[i] == -infinity, u[i] == infinity -- fb[i] = -f[i]
1195: . l[i] == -infinity, u[i] finite -- fb[i] = phi(u[i]-x[i], -f[i])
1196: . l[i] finite, u[i] == infinity -- fb[i] = phi(x[i]-l[i], f[i])
1197: . l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u]))
1198: - otherwise l[i] == u[i] -- fb[i] = l[i] - x[i]
1200: Level: advanced
1202: .seealso TaoMat::D_Fischer()
1203: @*/
1204: int TaoVec::Fischer(TaoVec* xx, TaoVec* ff, TaoVec* ll, TaoVec* uu)
1205: {
1206: int i,nn1,nn2,nn3,nn4,nn5,info;
1207: TaoScalar *v,*x,*f,*l,*u;
1209: TaoFunctionBegin;
1210: info = this->GetArray(&v,&nn1); CHKERRQ(info);
1211: info = xx->GetArray(&x,&nn2); CHKERRQ(info);
1212: info = ff->GetArray(&f,&nn3); CHKERRQ(info);
1213: info = ll->GetArray(&l,&nn4); CHKERRQ(info);
1214: info = uu->GetArray(&u,&nn5); CHKERRQ(info);
1215: if (nn1!=nn2 || nn2!=nn3 || nn3!=nn4 || nn4!=nn5) {
1216: TaoFunctionReturn(1);
1217: }
1219: for (i=0;i<nn1;i++) {
1221: if ((l[i] <= -TAO_INFINITY) && (u[i] >= TAO_INFINITY)) {
1222: v[i] = -f[i];
1223: }
1224: else if (l[i] <= -TAO_INFINITY) {
1225: v[i] = -fischer(u[i] - x[i], -f[i]);
1226: }
1227: else if (u[i] >= TAO_INFINITY) {
1228: v[i] = fischer(x[i] - l[i], f[i]);
1229: }
1230: else if (l[i] == u[i]) {
1231: v[i] = l[i] - x[i];
1232: }
1233: else {
1234: v[i] = fischer(u[i] - x[i], -f[i]);
1235: v[i] = fischer(x[i] - l[i], v[i]);
1236: }
1238: }
1240: info = this->RestoreArray(&v,&nn1);CHKERRQ(info);
1241: info = xx->RestoreArray(&x,&nn2);CHKERRQ(info);
1242: info = ff->RestoreArray(&f,&nn3);CHKERRQ(info);
1243: info = ll->RestoreArray(&l,&nn4);CHKERRQ(info);
1244: info = uu->RestoreArray(&u,&nn5);CHKERRQ(info);
1246: TaoFunctionReturn(0);
1247: }
1249: inline static double sfischer(double a, double b, double c)
1250: {
1252: #ifdef TODD
1254: if (TaoAbsScalar(a) > TaoAbsScalar(b)) {
1255: return sqrt(a*a + b*b + 2.0*c*c) - a - b;
1256: }
1257: return sqrt(a*a + b*b + 2.0*c*c) - b - a;
1259: #else
1261: // Method suggested by Bob Vanderbei
1263: if (a + b <= 0) {
1264: return sqrt(a*a + b*b + 2.0*c*c) - (a + b);
1265: }
1266: return 2.0*(c*c - a*b) / (sqrt(a*a + b*b + 2.0*c*c) + (a + b));
1268: #endif
1270: }
1274: /*@C
1275: SFischer - Evaluates the Smoothed Fischer-Burmeister function for
1276: complementarity problems.
1278: Input Parameters:
1279: + xx - current point
1280: . ff - function evaluated at x
1281: . ll - lower bounds
1282: . uu - upper bounds
1283: - mu - smoothing parameter
1285: Notes:
1286: The Smoothed Fischer-Burmeister function is defined as
1287: $ phi(a,b) := sqrt(a*a + b*b + 2*mu*mu) - a - b
1288: and is used reformulate a complementarity problem as a semismooth
1289: system of equations.
1291: The result of this function is done by cases:
1292: + l[i] == -infinity, u[i] == infinity -- fb[i] = -f[i] - 2*mu*x[i]
1293: . l[i] == -infinity, u[i] finite -- fb[i] = phi(u[i]-x[i], -f[i], mu)
1294: . l[i] finite, u[i] == infinity -- fb[i] = phi(x[i]-l[i], f[i], mu)
1295: . l[i] finite < u[i] finite -- fb[i] = phi(x[i]-l[i], phi(u[i]-x[i], -f[u], mu), mu)
1296: - otherwise l[i] == u[i] -- fb[i] = l[i] - x[i]
1298: Level: advanced
1300: .seealso TaoMat::SD_Fischer()
1301: @*/
1302: int TaoVec::SFischer(TaoVec* xx, TaoVec* ff, TaoVec* ll, TaoVec* uu, double mu)
1303: {
1305: int i, nn1, nn2, nn3, nn4, nn5, info;
1306: TaoScalar *v, *x, *f, *l, *u;
1308: TaoFunctionBegin;
1310: if ((mu >= -TAO_EPSILON) && (mu <= TAO_EPSILON)) {
1311: Fischer(xx, ff, ll, uu);
1312: }
1313: else {
1314: info = this->GetArray(&v, &nn1); CHKERRQ(info);
1315: info = xx->GetArray(&x, &nn2); CHKERRQ(info);
1316: info = ff->GetArray(&f, &nn3); CHKERRQ(info);
1317: info = ll->GetArray(&l, &nn4); CHKERRQ(info);
1318: info = uu->GetArray(&u, &nn5); CHKERRQ(info);
1320: if (nn1!=nn2 || nn2!=nn3 || nn3!=nn4 || nn4!=nn5) {
1321: TaoFunctionReturn(1);
1322: }
1324: for (i = 0; i < nn1; ++i) {
1325: if ((l[i] <= -TAO_INFINITY) && (u[i] >= TAO_INFINITY)) {
1326: v[i] = -f[i] - mu*x[i];
1327: }
1328: else if (l[i] <= -TAO_INFINITY) {
1329: v[i] = -sfischer(u[i] - x[i], -f[i], mu);
1330: }
1331: else if (u[i] >= TAO_INFINITY) {
1332: v[i] = sfischer(x[i] - l[i], f[i], mu);
1333: }
1334: else if (l[i] == u[i]) {
1335: v[i] = l[i] - x[i];
1336: }
1337: else {
1338: v[i] = sfischer(u[i] - x[i], -f[i], mu);
1339: v[i] = sfischer(x[i] - l[i], v[i], mu);
1340: }
1341: }
1343: info = this->RestoreArray(&v, &nn1); CHKERRQ(info);
1344: info = xx->RestoreArray(&x, &nn2); CHKERRQ(info);
1345: info = ff->RestoreArray(&f, &nn3); CHKERRQ(info);
1346: info = ll->RestoreArray(&l, &nn4); CHKERRQ(info);
1347: info = uu->RestoreArray(&u, &nn5); CHKERRQ(info);
1348: }
1349: TaoFunctionReturn(0);
1350: }