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