Actual source code: tvecsingle.c
1: #include "tvecsingle.h"
2: #include "tao_general.h"
3: #include "stdio.h"
5: TaoVecFloatArray::TaoVecFloatArray(int nn):TaoVec(),n(nn){
6: v=new float[nn];
7: fallocated=1;
8: return;
9: }
11: TaoVecFloatArray::TaoVecFloatArray(int nn, float *ff):TaoVec(),n(nn){
12: v=ff;
13: fallocated=0;
14: return;
15: }
17: int TaoVecFloatArray::Clone( TaoVec** tv ){
19: *tv = new TaoVecFloatArray(this->n);
20: int info = (*tv)->CopyFrom(this);CHKERRQ(info);
21: return 0;
22: }
24: int TaoVecFloatArray::GetFloats(float **fptr, int *nn){
25: *fptr=v;
26: *nn=n;
27: return 0;
28: }
30: int TaoVecFloatArray::RestoreFloats(float **fptr, int *nn){
31: *fptr=0;
32: *nn=0;
33: return 0;
34: }
36: int TaoVecFloatArray::GetDimension(int *nn){
37: *nn=n;
38: return 0;
39: }
41: int TaoVecFloatArray::Compatible(TaoVec *tv, TaoTruth *flag){
42: int nn,info;
43: float *fptr;
44: TaoVecFloatArray* vv = (TaoVecFloatArray*)(tv);
46: info = vv->GetData(&fptr,&nn);
47: if (info==0 && nn == n) *flag=TAO_TRUE;
48: else *flag=TAO_FALSE;
49: return 0;
50: }
52: int TaoVecFloatArray::View(){
53: for (int i=0;i<n;++i)
54: printf(" %4.2e \n ",v[i]);
55: return 0;
56: }
61: int TaoVecFloatArray::SetToConstant( double c ){
62: int i,nn,info;
63: float *tptr;
65: TaoFunctionBegin;
66: info = this->GetFloats(&tptr,&nn);CHKERRQ(info);
67: for (i=0;i<nn;i++){ tptr[i]=c; }
68: info = this->RestoreFloats(&tptr,&nn);CHKERRQ(info);
69: TaoFunctionReturn(0);
70: }
74: int TaoVecFloatArray::SetToZero(){
75: int i,nn,info;
76: float *tptr;
78: TaoFunctionBegin;
79: info = this->GetFloats(&tptr,&nn);CHKERRQ(info);
80: for (i=0;i<nn;i++){ tptr[i]=0; }
81: info = this->RestoreFloats(&tptr,&nn);CHKERRQ(info);
82: TaoFunctionReturn(0);
83: }
87: int TaoVecFloatArray::CopyFrom( TaoVec* tv ){
88: int i,nn1,nn2,info;
89: float *tptr1,*tptr2;
90: TaoVecFloatArray* vv = (TaoVecFloatArray*)(tv);
92: TaoFunctionBegin;
93: info = this->GetFloats(&tptr1,&nn1);CHKERRQ(info);
94: info = vv->GetFloats(&tptr2,&nn2);CHKERRQ(info);
96: for (i=0;i<nn1;i++){ tptr1[i]=tptr2[i]; }
98: info = vv->RestoreFloats(&tptr2,&nn2);CHKERRQ(info);
99: info = this->RestoreFloats(&tptr1,&nn1);CHKERRQ(info);
100: TaoFunctionReturn(0);
101: }
105: int TaoVecFloatArray::ScaleCopyFrom( double a, TaoVec* tv ){
106: int i,nn1,nn2,info;
107: float *tptr1,*tptr2;
108: TaoVecFloatArray* vv = (TaoVecFloatArray*)(tv);
110: TaoFunctionBegin;
111: info = this->GetFloats(&tptr1,&nn1);CHKERRQ(info);
112: info = vv->GetFloats(&tptr2,&nn2);CHKERRQ(info);
113: for (i=0;i<nn1;i++){ tptr1[i]=a*tptr2[i]; }
114: info = vv->RestoreFloats(&tptr2,&nn2);CHKERRQ(info);
115: info = this->RestoreFloats(&tptr1,&nn1);CHKERRQ(info);
116: TaoFunctionReturn(0);
117: }
121: int TaoVecFloatArray::NormInfinity(double *vnorm){
122: int i,nn,info;
123: float dd=0, *vv;
125: TaoFunctionBegin;
126: info = this->GetFloats(&vv,&nn);CHKERRQ(info);
127: for (i=0;i<nn;i++){
128: if (vv[i]<0 && vv[i]<-dd) dd=-vv[i];
129: else if (vv[i]>0 && vv[i]>dd) dd=vv[i];
130: }
131: info = this->RestoreFloats(&vv,&nn);CHKERRQ(info);
132: *vnorm=dd;
133: TaoFunctionReturn(0);
134: }
138: int TaoVecFloatArray::Norm1(double *vnorm){
139: int i,nn,info;
140: float dd=0, *vv;
142: TaoFunctionBegin;
143: info = this->GetFloats(&vv,&nn);CHKERRQ(info);
144: for (i=0;i<nn;i++){
145: if (vv[i]<0) dd-=vv[i];
146: else dd+=vv[i];
147: }
148: info = this->RestoreFloats(&vv,&nn);CHKERRQ(info);
149: TaoFunctionReturn(0);
150: }
154: int TaoVecFloatArray::Norm2(double *vnorm){
155: int i,nn,info;
156: float dd=0, *vv;
158: TaoFunctionBegin;
159: info = this->GetFloats(&vv,&nn);CHKERRQ(info);
160: for (i=0;i<nn;i++) dd+=vv[i]*vv[i];
161: *vnorm=sqrt(dd);
162: info = this->RestoreFloats(&vv,&nn);CHKERRQ(info);
163: TaoFunctionReturn(0);
164: }
168: int TaoVecFloatArray::Norm2squared(double *vnorm2){
169: int i,nn,info;
170: float dd=0, *vv;
172: TaoFunctionBegin;
173: info = this->GetFloats(&vv,&nn);CHKERRQ(info);
174: for (i=0;i<nn;i++) dd+=vv[i]*vv[i];
175: *vnorm2=dd;
176: info = this->RestoreFloats(&vv,&nn);CHKERRQ(info);
177: TaoFunctionReturn(0);
178: }
182: int TaoVecFloatArray::Scale( double alpha ){
183: int i,nn,info;
184: float *vv;
186: TaoFunctionBegin;
187: info = this->GetFloats(&vv,&nn);CHKERRQ(info);
188: for (i=0;i<nn;i++) vv[i]*=alpha;
189: info = this->RestoreFloats(&vv,&nn);CHKERRQ(info);
190: TaoFunctionReturn(0);
191: }
195: int TaoVecFloatArray::Axpy( double alpha, TaoVec* tv ){
196: int i,nn1,nn2,info;
197: float *tptr1,*tptr2;
198: TaoVecFloatArray* xx = (TaoVecFloatArray*)(tv);
200: TaoFunctionBegin;
201: info = this->GetFloats(&tptr1,&nn1);CHKERRQ(info);
202: info = xx->GetFloats(&tptr2,&nn2);CHKERRQ(info);
203: for (i=0;i<nn1;i++){ tptr1[i]+= alpha * tptr2[i]; }
204: info = this->RestoreFloats(&tptr1,&nn1);CHKERRQ(info);
205: info = xx->RestoreFloats(&tptr2,&nn2);CHKERRQ(info);
206: TaoFunctionReturn(0);
207: }
211: int TaoVecFloatArray::Axpby( double alpha, TaoVec* tv, double beta ){
212: int i,nn1,nn2,info;
213: float *tptr1,*tptr2;
214: TaoVecFloatArray* xx = (TaoVecFloatArray*)(tv);
216: TaoFunctionBegin;
217: info = this->GetFloats(&tptr1,&nn1);CHKERRQ(info);
218: info = xx->GetFloats(&tptr2,&nn2);CHKERRQ(info);
219: for (i=0;i<nn1;i++){ tptr1[i] = beta * tptr1[i] + alpha * tptr2[i]; }
220: info = this->RestoreFloats(&tptr1,&nn1);CHKERRQ(info);
221: info = xx->RestoreFloats(&tptr2,&nn2);CHKERRQ(info);
222: TaoFunctionReturn(0);
223: }
227: int TaoVecFloatArray::Aypx( double alpha, TaoVec* tv ){
228: int i,nn1,nn2,info;
229: float *tptr1,*tptr2;
230: TaoVecFloatArray* xx = (TaoVecFloatArray*)(tv);
232: TaoFunctionBegin;
233: info = this->GetFloats(&tptr1,&nn1);CHKERRQ(info);
234: info = xx->GetFloats(&tptr2,&nn2);CHKERRQ(info);
235: for (i=0;i<nn1;i++){ tptr1[i] = alpha * tptr1[i] + tptr2[i]; }
236: info = this->RestoreFloats(&tptr1,&nn1);CHKERRQ(info);
237: info = xx->RestoreFloats(&tptr2,&nn2);CHKERRQ(info);
238: TaoFunctionReturn(0);
239: }
240:
243: int TaoVecFloatArray::AddConstant( double alpha ){
244: int i,nn,info;
245: float *vv;
247: TaoFunctionBegin;
248: info = this->GetFloats(&vv,&nn);CHKERRQ(info);
249: for (i=0;i<nn;i++) vv[i]+=alpha;
250: info = this->RestoreFloats(&vv,&nn);CHKERRQ(info);
251: TaoFunctionReturn(0);
252: }
256: int TaoVecFloatArray::Dot( TaoVec* tv, double *vDotv ){
257: int i,nn1,nn2,info;
258: float dd=0,*tptr1,*tptr2;
259: TaoVecFloatArray* vv = (TaoVecFloatArray*)(tv);
261: TaoFunctionBegin;
262: info = this->GetFloats(&tptr1,&nn1);CHKERRQ(info);
263: info = vv->GetFloats(&tptr2,&nn2);CHKERRQ(info);
264: for (i=0;i<nn1;i++) dd+=tptr1[i]*tptr2[i];
265: info = this->RestoreFloats(&tptr1,&nn1);CHKERRQ(info);
266: info = vv->RestoreFloats(&tptr2,&nn2);CHKERRQ(info);
267: *vDotv=dd;
268: TaoFunctionReturn(0);
269: }
273: int TaoVecFloatArray::Negate(){
274: int i,nn,info;
275: float *vv;
277: TaoFunctionBegin;
278: info = this->GetFloats(&vv,&nn);CHKERRQ(info);
279: for (i=0;i<nn;i++) vv[i]=-vv[i];
280: info = this->RestoreFloats(&vv,&nn);CHKERRQ(info);
281: TaoFunctionReturn(0);
282: }
286: int TaoVecFloatArray::Reciprocal(){
287: int i,nn,info;
288: float *vv;
290: TaoFunctionBegin;
291: info = this->GetFloats(&vv,&nn);CHKERRQ(info);
292: for (i=0;i<nn;i++){
293: if (vv[i]!=0) vv[i]= 1.0/vv[i];
294: else vv[i]=TAO_INFINITY;
295: }
296: info = this->RestoreFloats(&vv,&nn);CHKERRQ(info);
297: TaoFunctionReturn(0);
298: }
302: int TaoVecFloatArray::PointwiseMultiply( TaoVec* tv, TaoVec* tw ){
303: int i,nn1,nn2,nn3,info;
304: float *tptr1,*tptr2,*tptr3;
305: TaoVecFloatArray* vv = (TaoVecFloatArray*)(tv);
306: TaoVecFloatArray* ww = (TaoVecFloatArray*)(tw);
308: TaoFunctionBegin;
309: info = this->GetFloats(&tptr1,&nn1);CHKERRQ(info);
310: info = vv->GetFloats(&tptr2,&nn2);CHKERRQ(info);
311: info = ww->GetFloats(&tptr3,&nn3);CHKERRQ(info);
312: for (i=0;i<nn1;i++) tptr1[i]=tptr2[i] * tptr3[i];
313: info = this->RestoreFloats(&tptr1,&nn1);CHKERRQ(info);
314: info = vv->RestoreFloats(&tptr2,&nn2);CHKERRQ(info);
315: info = ww->RestoreFloats(&tptr3,&nn3);CHKERRQ(info);
317: TaoFunctionReturn(0);
318: }
322: int TaoVecFloatArray::PointwiseDivide( TaoVec* tv , TaoVec* tw){
323: int i,nn1,nn2,nn3,info;
324: float *tptr1,*tptr2,*tptr3;
325: TaoVecFloatArray* vv = (TaoVecFloatArray*)(tv);
326: TaoVecFloatArray* ww = (TaoVecFloatArray*)(tw);
328: TaoFunctionBegin;
329: info = this->GetFloats(&tptr1,&nn1);CHKERRQ(info);
330: info = vv->GetFloats(&tptr2,&nn2);CHKERRQ(info);
331: info = ww->GetFloats(&tptr3,&nn3);CHKERRQ(info);
333: for (i=0;i<nn1;i++) tptr1[i]=tptr2[i] / tptr3[i];
334: info = this->RestoreFloats(&tptr1,&nn1);CHKERRQ(info);
335: info = vv->RestoreFloats(&tptr2,&nn2);CHKERRQ(info);
336: info = ww->RestoreFloats(&tptr3,&nn3);CHKERRQ(info);
338: TaoFunctionReturn(0);
339: }
343: int TaoVecFloatArray::Median( TaoVec* tv, TaoVec* tw, TaoVec* tx){
344: int i,nn1,nn2,nn3,nn4,info;
345: float *tptr1,*tptr2,*tptr3,*tptr4;
346: TaoVecFloatArray* vv = (TaoVecFloatArray*)(tv);
347: TaoVecFloatArray* ww = (TaoVecFloatArray*)(tw);
348: TaoVecFloatArray* xx = (TaoVecFloatArray*)(tx);
350: TaoFunctionBegin;
351: info = this->GetFloats(&tptr1,&nn1);CHKERRQ(info);
352: info = vv->GetFloats(&tptr2,&nn2);CHKERRQ(info);
353: info = ww->GetFloats(&tptr3,&nn3);CHKERRQ(info);
354: info = xx->GetFloats(&tptr4,&nn4);CHKERRQ(info);
356: for (i=0;i<nn1;i++){
357: if (tptr2[i]<=tptr3[i] && tptr3[i] <= tptr4[i]){
358: tptr1[i]=tptr3[i];
359: } else if (tptr4[i]<=tptr3[i] && tptr3[i] <= tptr2[i]){
360: tptr1[i]=tptr3[i];
361: } else if (tptr3[i]<=tptr2[i] && tptr2[i] <= tptr4[i]){
362: tptr1[i]=tptr2[i];
363: } else if (tptr4[i]<=tptr2[i] && tptr2[i] <= tptr3[i]){
364: tptr1[i]=tptr2[i];
365: } else {
366: tptr1[i]=tptr4[i];
367: }
368: }
369: info = this->RestoreFloats(&tptr1,&nn1);CHKERRQ(info);
370: info = vv->RestoreFloats(&tptr2,&nn2);CHKERRQ(info);
371: info = ww->RestoreFloats(&tptr3,&nn3);CHKERRQ(info);
372: info = xx->RestoreFloats(&tptr4,&nn4);CHKERRQ(info);
374: TaoFunctionReturn(0);
375: }
379: int TaoVecFloatArray::PointwiseMinimum( TaoVec* tv, TaoVec* tw){
380: int i,nn1,nn2,nn3,info;
381: float *tptr1,*tptr2,*tptr3;
382: TaoVecFloatArray* vv = (TaoVecFloatArray*)(tv);
383: TaoVecFloatArray* ww = (TaoVecFloatArray*)(tw);
385: TaoFunctionBegin;
386: info = this->GetFloats(&tptr1,&nn1);CHKERRQ(info);
387: info = vv->GetFloats(&tptr2,&nn2);CHKERRQ(info);
388: info = ww->GetFloats(&tptr3,&nn3);CHKERRQ(info);
390: for (i=0;i<nn1;i++) tptr1[i] = TaoMin( tptr2[i] , tptr3[i]);
392: info = this->RestoreFloats(&tptr1,&nn1);CHKERRQ(info);
393: info = vv->RestoreFloats(&tptr2,&nn2);CHKERRQ(info);
394: info = ww->RestoreFloats(&tptr3,&nn3);CHKERRQ(info);
396: TaoFunctionReturn(0);
397: }
401: int TaoVecFloatArray::PointwiseMaximum( TaoVec* tv, TaoVec* tw){
402: int i,nn1,nn2,nn3,info;
403: float *tptr1,*tptr2,*tptr3;
404: TaoVecFloatArray* vv = (TaoVecFloatArray*)(tv);
405: TaoVecFloatArray* ww = (TaoVecFloatArray*)(tw);
407: TaoFunctionBegin;
408: info = this->GetFloats(&tptr1,&nn1);CHKERRQ(info);
409: info = vv->GetFloats(&tptr2,&nn2);CHKERRQ(info);
410: info = ww->GetFloats(&tptr3,&nn3);CHKERRQ(info);
412: for (i=0;i<nn1;i++) tptr1[i] = TaoMax( tptr2[i] , tptr3[i]);
414: info = this->RestoreFloats(&tptr1,&nn1);CHKERRQ(info);
415: info = vv->RestoreFloats(&tptr2,&nn2);CHKERRQ(info);
416: info = ww->RestoreFloats(&tptr3,&nn3);CHKERRQ(info);
418: TaoFunctionReturn(0);
419: }
423: int TaoVecFloatArray::Waxpby ( double a, TaoVec* tx, double b, TaoVec* ty){
424: int i,nn1,nn2,nn3,info;
425: float *tptr1,*tptr2,*tptr3;
426: TaoVecFloatArray* xx = (TaoVecFloatArray*)(tx);
427: TaoVecFloatArray* yy = (TaoVecFloatArray*)(ty);
429: TaoFunctionBegin;
430: info = this->GetFloats(&tptr1,&nn1);CHKERRQ(info);
431: info = xx->GetFloats(&tptr2,&nn2);CHKERRQ(info);
432: info = yy->GetFloats(&tptr3,&nn3);CHKERRQ(info);
433: if (nn1!=nn2 || nn2!=nn3) {TaoFunctionReturn(1);}
434: for (i=0;i<nn1;i++){ tptr1[i] = a * tptr2[i] + b * tptr3[i]; }
435: info = this->RestoreFloats(&tptr1,&nn1);CHKERRQ(info);
436: info = xx->RestoreFloats(&tptr2,&nn2);CHKERRQ(info);
437: info = yy->RestoreFloats(&tptr3,&nn3);CHKERRQ(info);
438: TaoFunctionReturn(0);
439: }
443: int TaoVecFloatArray::AbsoluteValue(){
444: int i,nn,info;
445: float *vv;
447: TaoFunctionBegin;
448: info = this->GetFloats(&vv,&nn);CHKERRQ(info);
449: for (i=0;i<nn;i++){
450: vv[i]= TaoAbsScalar(vv[i]);
451: }
452: info = this->RestoreFloats(&vv,&nn);CHKERRQ(info);
453: TaoFunctionReturn(0);
454: }
458: int TaoVecFloatArray::MinElement(double*val){
459: int i,nn,info;
460: float dd=TAO_INFINITY,*vv;
462: TaoFunctionBegin;
463: info = this->GetFloats(&vv,&nn);CHKERRQ(info);
464: for (i=0;i<nn;i++){
465: if (vv[i]<dd) dd=vv[i];
466: }
467: info = this->RestoreFloats(&vv,&nn);CHKERRQ(info);
468: *val = dd;
469: TaoFunctionReturn(0);
470: }
474: int TaoVecFloatArray::StepMax(TaoVec*tv,double*step1){
475: int i,nn1,nn2,info;
476: float *xx,*dx;
477: float stepmax1=TAO_INFINITY;
478: TaoVecFloatArray* vv = (TaoVecFloatArray*)(tv);
480: TaoFunctionBegin;
481: info = this->GetFloats(&xx,&nn1);CHKERRQ(info);
482: info = vv->GetFloats(&dx,&nn2);CHKERRQ(info);
483:
484: for (i=0;i<nn1;i++){
485: if (xx[i] < 0){
486: TaoFunctionReturn(1);
487: } else if (dx[i]<0){
488: stepmax1=TaoMin(stepmax1,-xx[i]/dx[i]);
489: }
490: }
491: *step1=stepmax1;
493: info = vv->RestoreFloats(&dx,&nn2);CHKERRQ(info);
494: info = this->RestoreFloats(&xx,&nn1);CHKERRQ(info);
495: TaoFunctionReturn(0);
496: }
500: int TaoVecFloatArray::StepMax2(TaoVec*tv,TaoVec*txl,TaoVec*txu, double*step2){
501: int i,nn1,nn2,info;
502: float *xx,*dx,*xl,*xu;
503: float stepmax2=0;
504: TaoVecFloatArray* vv = (TaoVecFloatArray*)(tv);
505: TaoVecFloatArray* xxll = (TaoVecFloatArray*)(txl);
506: TaoVecFloatArray* xxuu = (TaoVecFloatArray*)(txu);
508: TaoFunctionBegin;
509: info = this->GetFloats(&xx,&nn1);CHKERRQ(info);
510: info = vv->GetFloats(&dx,&nn2);CHKERRQ(info);
511: info = xxll->GetFloats(&xl,&nn2);CHKERRQ(info);
512: info = xxuu->GetFloats(&xu,&nn2);CHKERRQ(info);
513:
514: for (i=0;i<nn1;i++){
515: if (dx[i] > 0){
516: stepmax2=TaoMax(stepmax2,(xu[i]-xx[i])/dx[i]);
517: } else if (dx[i]<0){
518: stepmax2=TaoMax(stepmax2,(xl[i]-xx[i])/dx[i]);
519: }
520: }
521: info = this->RestoreFloats(&xx,&nn1);CHKERRQ(info);
522: info = vv->RestoreFloats(&dx,&nn2);CHKERRQ(info);
523: info = xxll->RestoreFloats(&xl,&nn2);CHKERRQ(info);
524: info = xxuu->RestoreFloats(&xu,&nn2);CHKERRQ(info);
525:
526: *step2=stepmax2;
527: TaoFunctionReturn(0);
528: }
533: int TaoVecFloatArray::BoundGradientProjection(TaoVec*tg,TaoVec*txl,TaoVec*tx, TaoVec*txu){
534: int i,nn1,nn2,nn3,nn4,nn5,info;
535: float *xptr,*xlptr,*xuptr,*gptr,*gpptr;
536: TaoVecFloatArray* gg = (TaoVecFloatArray*)(tg);
537: TaoVecFloatArray* xxll = (TaoVecFloatArray*)(txl);
538: TaoVecFloatArray* xx = (TaoVecFloatArray*)(tx);
539: TaoVecFloatArray* xxuu = (TaoVecFloatArray*)(txu);
541: TaoFunctionBegin;
542: info = this->GetFloats(&gpptr,&nn1);CHKERRQ(info);
543: info = gg->GetFloats(&gptr,&nn2);CHKERRQ(info);
544: info = xxll->GetFloats(&xlptr,&nn3);CHKERRQ(info);
545: info = xx->GetFloats(&xptr,&nn4);CHKERRQ(info);
546: info = xxuu->GetFloats(&xuptr,&nn5);CHKERRQ(info);
547: if (nn1!=nn2 || nn2!=nn3 || nn3!=nn4 || nn4!=nn5) {TaoFunctionReturn(1);}
549: for (i=0; i<nn1; i++){
551: gpptr[i] = gptr[i];
552: if (gpptr[i]>0 && xptr[i]<=xlptr[i]){
553: gpptr[i] = 0;
554: } else if (gpptr[i]<0 && xptr[i]>=xuptr[i]){
555: gpptr[i] = 0;
556: }
557: }
558: info = this->RestoreFloats(&gpptr,&nn1);CHKERRQ(info);
559: info = gg->RestoreFloats(&gptr,&nn2);CHKERRQ(info);
560: info = xxll->RestoreFloats(&xlptr,&nn3);CHKERRQ(info);
561: info = xx->RestoreFloats(&xptr,&nn4);CHKERRQ(info);
562: info = xxuu->RestoreFloats(&xuptr,&nn5);CHKERRQ(info);
564: TaoFunctionReturn(0);
565: }
567: inline static float fischer(float a, float b)
568: {
570: #ifdef TODD
572: if (TaoAbsScalar(a) > TaoAbsScalar(b)) {
573: return sqrt(a*a + b*b) - a - b;
574: }
575: return sqrt(a*a + b*b) - b - a;
577: #else
579: // Method suggested by Bob Vanderbei
581: if (a + b <= 0) {
582: return sqrt(a*a + b*b) - (a + b);
583: }
584: return -2.0*a*b / (sqrt(a*a + b*b) + (a + b));
586: #endif
588: }
592: int TaoVecFloatArray::Fischer(TaoVec *tx, TaoVec *tf, TaoVec *tl, TaoVec *tu)
593: {
594: int i,nn1,nn2,nn3,nn4,nn5,info;
595: float *vv,*x,*f,*l,*u;
597: TaoVecFloatArray* xx = (TaoVecFloatArray*)(tx);
598: TaoVecFloatArray* ff = (TaoVecFloatArray*)(tf);
599: TaoVecFloatArray* ll = (TaoVecFloatArray*)(tl);
600: TaoVecFloatArray* uu = (TaoVecFloatArray*)(tu);
602: TaoFunctionBegin;
603: info = this->GetFloats(&vv,&nn1);CHKERRQ(info);
604: info = xx->GetFloats(&x,&nn2);CHKERRQ(info);
605: info = ff->GetFloats(&f,&nn3);CHKERRQ(info);
606: info = ll->GetFloats(&l,&nn4);CHKERRQ(info);
607: info = uu->GetFloats(&u,&nn5);CHKERRQ(info);
608: if (nn1!=nn2 || nn2!=nn3 || nn3!=nn4 || nn4!=nn5) {
609: TaoFunctionReturn(1);
610: }
612: for (i=0;i<nn1;i++) {
614: if ((l[i] <= -TAO_INFINITY) && (u[i] >= TAO_INFINITY)) {
615: vv[i] = -f[i];
616: }
617: else if (l[i] <= -TAO_INFINITY) {
618: vv[i] = -fischer(u[i] - x[i], -f[i]);
619: }
620: else if (u[i] >= TAO_INFINITY) {
621: vv[i] = fischer(x[i] - l[i], f[i]);
622: }
623: else if (l[i] == u[i]) {
624: vv[i] = l[i] - x[i];
625: }
626: else {
627: vv[i] = fischer(u[i] - x[i], -f[i]);
628: vv[i] = fischer(x[i] - l[i], vv[i]);
629: }
631: }
633: info = this->RestoreFloats(&vv,&nn1);CHKERRQ(info);
634: info = xx->RestoreFloats(&x,&nn2);CHKERRQ(info);
635: info = ff->RestoreFloats(&f,&nn3);CHKERRQ(info);
636: info = ll->RestoreFloats(&l,&nn4);CHKERRQ(info);
637: info = uu->RestoreFloats(&u,&nn5);CHKERRQ(info);
639: TaoFunctionReturn(0);
640: }
642: inline static float sfischer(float a, float b, float c)
643: {
645: #ifdef TODD
647: if (TaoAbsScalar(a) > TaoAbsScalar(b)) {
648: return sqrt(a*a + b*b + 2.0*c*c) - a - b;
649: }
650: return sqrt(a*a + b*b + 2.0*c*c) - b - a;
652: #else
654: // Method suggested by Bob Vanderbei
656: if (a + b <= 0) {
657: return sqrt(a*a + b*b + 2.0*c*c) - (a + b);
658: }
659: return 2.0*(c*c - a*b) / (sqrt(a*a + b*b + 2.0*c*c) + (a + b));
661: #endif
663: }
667: int TaoVecFloatArray::SFischer(TaoVec *tx, TaoVec *tf, TaoVec *tl, TaoVec *tu, double mu)
668: {
669: int i, nn1, nn2, nn3, nn4, nn5, info;
670: float *vv, *x, *f, *l, *u;
672: TaoVecFloatArray *xx = (TaoVecFloatArray *)(tx);
673: TaoVecFloatArray *ff = (TaoVecFloatArray *)(tf);
674: TaoVecFloatArray *ll = (TaoVecFloatArray *)(tl);
675: TaoVecFloatArray *uu = (TaoVecFloatArray *)(tu);
677: TaoFunctionBegin;
679: if ((mu >= -TAO_EPSILON) && (mu <= TAO_EPSILON)) {
680: Fischer(tx, tf, tl, tu);
681: }
682: else {
683: info = this->GetFloats(&vv, &nn1); CHKERRQ(info);
684: info = xx->GetFloats(&x, &nn2); CHKERRQ(info);
685: info = ff->GetFloats(&f, &nn3); CHKERRQ(info);
686: info = ll->GetFloats(&l, &nn4); CHKERRQ(info);
687: info = uu->GetFloats(&u, &nn5); CHKERRQ(info);
689: if (nn1!=nn2 || nn2!=nn3 || nn3!=nn4 || nn4!=nn5) {
690: TaoFunctionReturn(1);
691: }
693: for (i = 0; i < nn1; ++i) {
694: if ((l[i] <= -TAO_INFINITY) && (u[i] >= TAO_INFINITY)) {
695: vv[i] = -f[i] - mu*x[i];
696: }
697: else if (l[i] <= -TAO_INFINITY) {
698: vv[i] = -sfischer(u[i] - x[i], -f[i], mu);
699: }
700: else if (u[i] >= TAO_INFINITY) {
701: vv[i] = sfischer(x[i] - l[i], f[i], mu);
702: }
703: else if (l[i] == u[i]) {
704: vv[i] = l[i] - x[i];
705: }
706: else {
707: vv[i] = sfischer(u[i] - x[i], -f[i], mu);
708: vv[i] = sfischer(x[i] - l[i], vv[i], mu);
709: }
710: }
712: info = this->RestoreFloats(&vv, &nn1); CHKERRQ(info);
713: info = xx->RestoreFloats(&x, &nn2); CHKERRQ(info);
714: info = ff->RestoreFloats(&f, &nn3); CHKERRQ(info);
715: info = ll->RestoreFloats(&l, &nn4); CHKERRQ(info);
716: info = uu->RestoreFloats(&u, &nn5); CHKERRQ(info);
717: }
718: TaoFunctionReturn(0);
719: }