Actual source code: pvec.c
1: #include "tao_general.h"
2: #include "private/vecimpl.h" /*I "tao_solver.h" I*/
3: #include "petscksp.h"
9: int VecCompare(Vec V1,Vec V2, PetscTruth *flg){
10: int info;
11: PetscInt n1,n2,N1,N2;
16: info = VecGetSize(V1,&N1);CHKERRQ(info);
17: info = VecGetSize(V2,&N2);CHKERRQ(info);
18: info = VecGetLocalSize(V1,&n1);CHKERRQ(info);
19: info = VecGetLocalSize(V2,&n2);CHKERRQ(info);
20: if (N1==N2 && n1==n2)
21: *flg=PETSC_TRUE;
22: else
23: *flg=PETSC_FALSE;
25: return(0);
26: }
28: /* ---------------------------------------------------------- */
31: int VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
32: {
33: int i,info;
34: PetscInt n,low1,low2,low3,low4,high1,high2,high3,high4;
35: PetscScalar *v1,*v2,*v3,*vmed;
43: if (Vec1==Vec2 || Vec1==Vec3){
44: info=VecCopy(Vec1,VMedian); CHKERRQ(info);
45: return(0);
46: }
47: if (Vec2==Vec3){
48: info=VecCopy(Vec2,VMedian); CHKERRQ(info);
49: return(0);
50: }
58: info = VecGetOwnershipRange(Vec1, &low1, &high1); CHKERRQ(info);
59: info = VecGetOwnershipRange(Vec2, &low2, &high2); CHKERRQ(info);
60: info = VecGetOwnershipRange(Vec3, &low3, &high3); CHKERRQ(info);
61: info = VecGetOwnershipRange(VMedian, &low4, &high4); CHKERRQ(info);
62: if ( low1!= low2 || low1!= low3 || low1!= low4 ||
63: high1!= high2 || high1!= high3 || high1!= high4)
64: SETERRQ(1,"InCompatible vector local lengths");
66: info = VecGetArray(Vec1,&v1); CHKERRQ(info);
67: info = VecGetArray(Vec2,&v2); CHKERRQ(info);
68: info = VecGetArray(Vec3,&v3); CHKERRQ(info);
70: if ( VMedian != Vec1 && VMedian != Vec2 && VMedian != Vec3){
71: info = VecGetArray(VMedian,&vmed); CHKERRQ(info);
72: } else if ( VMedian==Vec1 ){
73: vmed=v1;
74: } else if ( VMedian==Vec2 ){
75: vmed=v2;
76: } else {
77: vmed=v3;
78: }
80: info=VecGetLocalSize(Vec1,&n); CHKERRQ(info);
82: for (i=0;i<n;i++){
83: vmed[i]=PetscMax(PetscMax(PetscMin(v1[i],v2[i]),PetscMin(v1[i],v3[i])),PetscMin(v2[i],v3[i]));
84: }
86: info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
87: info = VecRestoreArray(Vec2,&v2); CHKERRQ(info);
88: info = VecRestoreArray(Vec3,&v2); CHKERRQ(info);
89:
90: if (VMedian!=Vec1 && VMedian != Vec2 && VMedian != Vec3){
91: info = VecRestoreArray(VMedian,&vmed); CHKERRQ(info);
92: }
94: return(0);
95: }
97: /* ----------------------------------------------------------
100: int VecPointwiseMin(Vec Vec1, Vec Vec2, Vec VMin)
101: {
102: int i,n,low1,low2,low3,high1,high2,high3,info;
103: PetscScalar *v1,*v2,*vmin;
110: if (Vec1==Vec2){
111: info=VecCopy(Vec1,VMin); CHKERRQ(info);
112: return(0);
113: }
118: info = VecGetOwnershipRange(Vec1, &low1, &high1); CHKERRQ(info);
119: info = VecGetOwnershipRange(Vec2, &low2, &high2); CHKERRQ(info);
120: info = VecGetOwnershipRange(VMin, &low3, &high3); CHKERRQ(info);
122: if ( low1!= low2 || low1!= low3 || high1!= high2 || high1!= high3)
123: SETERRQ(1,"InCompatible vector local lengths");
125: info = VecGetArray(Vec1,&v1); CHKERRQ(info);
126: info = VecGetArray(Vec2,&v2); CHKERRQ(info);
128: if ( VMin != Vec1 && VMin != Vec2){
129: info = VecGetArray(VMin,&vmin); CHKERRQ(info);
130: } else if ( VMin==Vec1 ){
131: vmin=v1;
132: } else {
133: vmin =v2;
134: }
136: info=VecGetLocalSize(Vec1,&n); CHKERRQ(info);
138: for (i=0; i<n; i++){
139: vmin[i] = PetscMin(v1[i],v2[i]);
140: }
142: info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
143: info = VecRestoreArray(Vec2,&v2); CHKERRQ(info);
145: if (VMin!=Vec1 && VMin != Vec2){
146: info = VecRestoreArray(VMin,&vmin); CHKERRQ(info);
147: }
149: return(0);
150: }
151: */
152: /* ----------------------------------------------------------
155: int VecPointwiseMax(Vec Vec1, Vec Vec2, Vec VMax)
156: {
157: int i,n,low1,low2,low3,high1,high2,high3,info;
158: PetscScalar *v1,*v2,*vmax;
162: if (Vec1==Vec2){
163: info=VecCopy(Vec1,VMax); CHKERRQ(info);
164: return(0);
165: }
174: info = VecGetOwnershipRange(Vec1, &low1, &high1); CHKERRQ(info);
175: info = VecGetOwnershipRange(Vec2, &low2, &high2); CHKERRQ(info);
176: info = VecGetOwnershipRange(VMax, &low3, &high3); CHKERRQ(info);
178: if ( low1!= low2 || low1!= low3 || high1!= high2 || high1!= high3)
179: SETERRQ(1,"Vectors must be identically loaded over processors");
181: info = VecGetArray(Vec1,&v1); CHKERRQ(info);
182: info = VecGetArray(Vec2,&v2); CHKERRQ(info);
184: if ( VMax != Vec1 && VMax != Vec2){
185: info = VecGetArray(VMax,&vmax); CHKERRQ(info);
186: } else if ( VMax==Vec1 ){ vmax=v1;
187: } else { vmax =v2;
188: }
190: info=VecGetLocalSize(Vec1,&n); CHKERRQ(info);
192: for (i=0; i<n; i++){
193: vmax[i] = PetscMax(v1[i],v2[i]);
194: }
196: info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
197: info = VecRestoreArray(Vec2,&v2); CHKERRQ(info);
199: if ( VMax != Vec1 && VMax != Vec2){
200: info = VecRestoreArray(VMax,&vmax); CHKERRQ(info);
201: }
203: return(0);
204: }*/
208: inline static PetscScalar Fischer(PetscScalar a, PetscScalar b)
209: {
210: // Method suggested by Bob Vanderbei
211: if (a + b <= 0) {
212: return sqrt(a*a + b*b) - (a + b);
213: }
214: return -2.0*a*b / (sqrt(a*a + b*b) + (a + b));
215: }
219: int VecFischer(Vec X, Vec F, Vec L, Vec U, Vec FF)
220: {
221: PetscScalar *x, *f, *l, *u, *ff;
222: PetscScalar xval, fval, lval, uval;
223: int info, i;
224: PetscInt low[5], high[5], n;
233: info = VecGetOwnershipRange(X, low, high); CHKERRQ(info);
234: info = VecGetOwnershipRange(F, low + 1, high + 1); CHKERRQ(info);
235: info = VecGetOwnershipRange(L, low + 2, high + 2); CHKERRQ(info);
236: info = VecGetOwnershipRange(U, low + 3, high + 3); CHKERRQ(info);
237: info = VecGetOwnershipRange(FF, low + 4, high + 4); CHKERRQ(info);
239: for (i = 1; i < 4; ++i) {
240: if (low[0] != low[i] || high[0] != high[i])
241: SETERRQ(1,"Vectors must be identically loaded over processors");
242: }
244: info = VecGetArray(X, &x); CHKERRQ(info);
245: info = VecGetArray(F, &f); CHKERRQ(info);
246: info = VecGetArray(L, &l); CHKERRQ(info);
247: info = VecGetArray(U, &u); CHKERRQ(info);
248: info = VecGetArray(FF, &ff); CHKERRQ(info);
250: info = VecGetLocalSize(X, &n); CHKERRQ(info);
252: for (i = 0; i < n; ++i) {
253: xval = x[i]; fval = f[i];
254: lval = l[i]; uval = u[i];
256: if ((lval <= -TAO_INFINITY) && (uval >= TAO_INFINITY)) {
257: ff[i] = -fval;
258: }
259: else if (lval <= -TAO_INFINITY) {
260: ff[i] = -Fischer(uval - xval, -fval);
261: }
262: else if (uval >= TAO_INFINITY) {
263: ff[i] = Fischer(xval - lval, fval);
264: }
265: else if (lval == uval) {
266: ff[i] = lval - xval;
267: }
268: else {
269: fval = Fischer(uval - xval, -fval);
270: ff[i] = Fischer(xval - lval, fval);
271: }
272: }
273:
274: info = VecRestoreArray(X, &x); CHKERRQ(info);
275: info = VecRestoreArray(F, &f); CHKERRQ(info);
276: info = VecRestoreArray(L, &l); CHKERRQ(info);
277: info = VecRestoreArray(U, &u); CHKERRQ(info);
278: info = VecRestoreArray(FF, &ff); CHKERRQ(info);
280: return(0);
281: }
285: inline static PetscScalar SFischer(PetscScalar a, PetscScalar b, PetscScalar c)
286: {
287: // Method suggested by Bob Vanderbei
288: if (a + b <= 0) {
289: return sqrt(a*a + b*b + 2.0*c*c) - (a + b);
290: }
291: return 2.0*(c*c - a*b) / (sqrt(a*a + b*b + 2.0*c*c) + (a + b));
292: }
296: int VecSFischer(Vec X, Vec F, Vec L, Vec U, double mu, Vec FF)
297: {
298: PetscScalar *x, *f, *l, *u, *ff;
299: PetscScalar xval, fval, lval, uval;
300: int info, i;
301: PetscInt low[5], high[5], n;
310: info = VecGetOwnershipRange(X, low, high); CHKERRQ(info);
311: info = VecGetOwnershipRange(F, low + 1, high + 1); CHKERRQ(info);
312: info = VecGetOwnershipRange(L, low + 2, high + 2); CHKERRQ(info);
313: info = VecGetOwnershipRange(U, low + 3, high + 3); CHKERRQ(info);
314: info = VecGetOwnershipRange(FF, low + 4, high + 4); CHKERRQ(info);
316: for (i = 1; i < 4; ++i) {
317: if (low[0] != low[i] || high[0] != high[i])
318: SETERRQ(1,"Vectors must be identically loaded over processors");
319: }
321: info = VecGetArray(X, &x); CHKERRQ(info);
322: info = VecGetArray(F, &f); CHKERRQ(info);
323: info = VecGetArray(L, &l); CHKERRQ(info);
324: info = VecGetArray(U, &u); CHKERRQ(info);
325: info = VecGetArray(FF, &ff); CHKERRQ(info);
327: info = VecGetLocalSize(X, &n); CHKERRQ(info);
329: for (i = 0; i < n; ++i) {
330: xval = (*x++); fval = (*f++);
331: lval = (*l++); uval = (*u++);
333: if ((lval <= -TAO_INFINITY) && (uval >= TAO_INFINITY)) {
334: (*ff++) = -fval - mu*xval;
335: }
336: else if (lval <= -TAO_INFINITY) {
337: (*ff++) = -SFischer(uval - xval, -fval, mu);
338: }
339: else if (uval >= TAO_INFINITY) {
340: (*ff++) = SFischer(xval - lval, fval, mu);
341: }
342: else if (lval == uval) {
343: (*ff++) = lval - xval;
344: }
345: else {
346: fval = SFischer(uval - xval, -fval, mu);
347: (*ff++) = SFischer(xval - lval, fval, mu);
348: }
349: }
350: x -= n; f -= n; l -=n; u -= n; ff -= n;
352: info = VecRestoreArray(X, &x); CHKERRQ(info);
353: info = VecRestoreArray(F, &f); CHKERRQ(info);
354: info = VecRestoreArray(L, &l); CHKERRQ(info);
355: info = VecRestoreArray(U, &u); CHKERRQ(info);
356: info = VecRestoreArray(FF, &ff); CHKERRQ(info);
357: return(0);
358: }
362: int VecBoundProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP){
364: int i,info;
365: PetscInt n;
366: PetscScalar *xptr,*xlptr,*xuptr,*gptr,*gpptr;
367: PetscScalar xval,gpval;
369: /* Project variables at the lower and upper bound */
372: info = VecGetLocalSize(X,&n); CHKERRQ(info);
374: info=VecGetArray(X,&xptr); CHKERRQ(info);
375: info=VecGetArray(XL,&xlptr); CHKERRQ(info);
376: info=VecGetArray(XU,&xuptr); CHKERRQ(info);
377: info=VecGetArray(G,&gptr); CHKERRQ(info);
378: if (G!=GP){
379: info=VecGetArray(GP,&gpptr); CHKERRQ(info);
380: } else { gpptr=gptr; }
382: for (i=0; i<n; ++i){
383: gpval = gptr[i]; xval = xptr[i];
385: if (gpval>0 && xval<=xlptr[i]){
386: gpval = 0;
387: } else if (gpval<0 && xval>=xuptr[i]){
388: gpval = 0;
389: }
390: gpptr[i] = gpval;
392: /*
393: if (xlptr[i] >= xuptr[i]){
394: gpptr[i]=0.0;
395: } else if (xptr[i] <= xlptr[i]+eps){
396: gpptr[i] = PetscMin(gptr[i],zero);
397: } else if (xptr[i] >= xuptr[i]-eps){
398: gpptr[i] = PetscMax(gptr[i],zero);
399: } else {
400: gpptr[i] = gptr[i];
401: }
402: */
403: }
405: info=VecRestoreArray(X,&xptr); CHKERRQ(info);
406: info=VecRestoreArray(XL,&xlptr); CHKERRQ(info);
407: info=VecRestoreArray(XU,&xuptr); CHKERRQ(info);
408: info=VecRestoreArray(G,&gptr); CHKERRQ(info);
409: if (G!=GP){
410: info=VecRestoreArray(GP,&gpptr); CHKERRQ(info);
411: }
412: return(0);
413: }
417: int VecISAXPY(Vec Y, PetscScalar alpha, Vec X, IS YIS){
419: int i,info;
420: PetscInt in1,in2,xlow,xhigh,ylow,yhigh;
421: const PetscInt *yis;
422: PetscScalar *x,*y;
425: info=ISGetLocalSize(YIS,&in1); CHKERRQ(info);
426: info=VecGetLocalSize(X,&in2); CHKERRQ(info);
428: if ( in1 != in2 )
429: SETERRQ(1,"Index set and X vector must be identically loaded over processors");
431: info=VecGetOwnershipRange(X, &xlow, &xhigh); CHKERRQ(info);
432: info=VecGetOwnershipRange(Y, &ylow, &yhigh); CHKERRQ(info);
434: info=ISGetIndices(YIS, &yis); CHKERRQ(info);
436: info=VecGetArray(X,&x); CHKERRQ(info);
437: if (X != Y){
438: info=VecGetArray(Y,&y); CHKERRQ(info);
439: } else {
440: y=x;
441: }
443: for (i=0; i<in1; i++){
444: if (yis[i]>=ylow && yis[i]<yhigh && i+xlow < xhigh){
445: y[yis[i]-ylow] = y[yis[i]-ylow] + (alpha)*x[i];
446: } else {
447: SETERRQ(1,"IS index out of range");
448: }
449: }
450:
451: info=ISRestoreIndices(YIS, &yis); CHKERRQ(info);
453: info=VecRestoreArray(X,&x); CHKERRQ(info);
454: if (X != Y){
455: info=VecRestoreArray(Y,&y); CHKERRQ(info);
456: }
458: info = PetscLogFlops(2*in1); CHKERRQ(info);
460: return(0);
461: }
465: int VecCreateSubVec(Vec A,IS S,Vec *B)
466: {
467: int i,info;
468: PetscInt nloc,n,low,high,nnloc;
469: const PetscInt *ptrS;
470: PetscScalar zero=0.0;
471: PetscScalar *ptrB,*ptrA;
472: const VecType type_name;
473: MPI_Comm comm;
478: info=ISGetLocalSize(S,&nloc);CHKERRQ(info);
479: info=ISGetSize(S,&n);CHKERRQ(info);
480: info=PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(info);
481: info=VecCreate(comm,B);CHKERRQ(info);
482: info=VecSetSizes(*B,nloc,n);CHKERRQ(info);
483: info=VecGetType(A,&type_name);CHKERRQ(info);
484: info=VecSetType(*B,type_name);CHKERRQ(info);
485: info=VecSet(*B, zero);CHKERRQ(info);
486:
487: info=VecGetOwnershipRange(A,&low,&high);CHKERRQ(info);
488: info=VecGetLocalSize(A,&nnloc);CHKERRQ(info);
489: info=VecGetArray(A,&ptrA);CHKERRQ(info);
490: info=VecGetArray(*B,&ptrB);CHKERRQ(info);
491: info=ISGetIndices(S,&ptrS);CHKERRQ(info);
492: for (i=0;i<nloc;i++){ ptrB[i]=ptrA[ptrS[i]-low]; }
493: info=VecRestoreArray(A,&ptrA);CHKERRQ(info);
494: info=VecRestoreArray(*B,&ptrB);CHKERRQ(info);
495: info=ISRestoreIndices(S,&ptrS);CHKERRQ(info);
497: return(0);
498: }
502: int VecStepMax(Vec X, Vec DX, PetscReal*step){
503: int i,info;
504: PetscInt nn;
505: PetscScalar stepmax=1.0e300;
506: PetscScalar *xx,*dx;
507: double t1,t2;
510: info = VecGetLocalSize(X,&nn);CHKERRQ(info);
511: info = VecGetArray(X,&xx);CHKERRQ(info);
512: info = VecGetArray(DX,&dx);CHKERRQ(info);
513: for (i=0;i<nn;i++){
514: if (xx[i] < 0){
515: SETERRQ(1,"Vector must be positive");
516: } else if (dx[i]<0){ stepmax=PetscMin(stepmax,-xx[i]/dx[i]);
517: }
518: }
519: info = VecRestoreArray(X,&xx);CHKERRQ(info);
520: info = VecRestoreArray(DX,&dx);CHKERRQ(info);
521: t1=(double)stepmax;
522: info = MPI_Allreduce(&t1,&t2,1,MPI_DOUBLE,MPI_MIN,((PetscObject)X)->comm);
523: CHKERRQ(info);
524: *step=(PetscReal)t2;
525: return(0);
526: }
530: int MatCreateProjectionOperator(Vec A, IS S , Vec B, Mat *MM){
532: PetscInt m,n,M,N,low,high;
533: int info,i;
534: const PetscInt *iptr;
535: PetscScalar one=1.0;
536: Mat MMM;
544: info = VecGetSize(B,&M);CHKERRQ(info);
545: info = VecGetLocalSize(B,&m);CHKERRQ(info);
546: info = VecGetSize(A,&N);CHKERRQ(info);
547: info = VecGetLocalSize(A,&n);CHKERRQ(info);
549: info = MatCreateMPIAIJ(((PetscObject)A)->comm,m,n,M,N,1,PETSC_NULL,1,PETSC_NULL,&MMM);
550: CHKERRQ(info);
551: /*
552: info = MatCreateMPIBDiag(((PetscObject)pv)->comm,m,M,N,M,1,PETSC_NULL,PETSC_NULL,&A);
553: */
554:
555: info = ISGetSize(S,&n);CHKERRQ(info);
556: if (n!=M){
557: SETERRQ(1,"Size of index set does not match the second vector.");
558: }
559:
560: info = ISGetLocalSize(S,&n);CHKERRQ(info);
561: if (n!=m){
562: SETERRQ(1,"Local size of index set does not match the second vector.");
563: }
564: info=VecGetOwnershipRange(B,&low,&high);CHKERRQ(info);
565: info = ISGetIndices(S,&iptr);CHKERRQ(info);
566: for (i=0; i<n; i++){
567: info=MatSetValue(MMM,low+i,iptr[i],one,INSERT_VALUES);CHKERRQ(info);
568: }
569: info = ISRestoreIndices(S,&iptr);CHKERRQ(info);
571: info = MatAssemblyBegin(MMM,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
572: info = MatAssemblyEnd(MMM,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
576: *MM = MMM;
578: return(0);
579: }
583: int VecStepMax2(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *step2){
585: int i,info;
586: PetscInt nn1,nn2;
587: PetscScalar *xx,*dx,*xl,*xu;
588: PetscScalar stepmax2=0;
589: double t1,t2;
591: TaoFunctionBegin;
592: info = VecGetArray(X,&xx);CHKERRQ(info);
593: info = VecGetArray(XL,&xl);CHKERRQ(info);
594: info = VecGetArray(XU,&xu);CHKERRQ(info);
595: info = VecGetArray(DX,&dx);CHKERRQ(info);
596: info = VecGetLocalSize(X,&nn1);CHKERRQ(info);
597: info = VecGetLocalSize(DX,&nn2);CHKERRQ(info);
599: for (i=0;i<nn1;i++){
600: if (dx[i] > 0){
601: stepmax2=PetscMax(stepmax2,(xu[i]-xx[i])/dx[i]);
602: } else if (dx[i]<0){
603: stepmax2=PetscMax(stepmax2,(xl[i]-xx[i])/dx[i]);
604: }
605: }
606: info = VecRestoreArray(X,&xx);CHKERRQ(info);
607: info = VecRestoreArray(XL,&xl);CHKERRQ(info);
608: info = VecRestoreArray(XU,&xu);CHKERRQ(info);
609: info = VecRestoreArray(DX,&dx);CHKERRQ(info);
611: t1=(double)stepmax2;
612: info = MPI_Allreduce(&t1,&t2,1,MPI_DOUBLE,MPI_MAX,((PetscObject)X)->comm);
613: CHKERRQ(info);
614: *step2=(PetscReal)t2;
616: TaoFunctionReturn(0);
617: }
622: int VecBoundStepInfo(Vec X, Vec XL, Vec XU, Vec DX, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax){
623: int info,i;
624: PetscInt n;
625: PetscScalar *x,*xl,*xu,*dx;
626: double tt,t1=1.0e+20, t2=1.0e+20, t3=0;
627: double tt1,tt2,tt3;
628: MPI_Comm comm;
629:
631: info=VecGetArray(X,&x);CHKERRQ(info);
632: info=VecGetArray(XL,&xl);CHKERRQ(info);
633: info=VecGetArray(XU,&xu);CHKERRQ(info);
634: info=VecGetArray(DX,&dx);CHKERRQ(info);
635: info = VecGetLocalSize(X,&n);CHKERRQ(info);
636: for (i=0;i<n;i++){
637: if (dx[i]>0){
638: tt=(xu[i]-x[i])/dx[i];
639: t1=PetscMin(t1,tt);
640: if (t1>0){
641: t2=PetscMin(t2,tt);
642: }
643: t3=PetscMax(t3,tt);
644: } else if (dx[i]<0){
645: tt=(xl[i]-x[i])/dx[i];
646: t1=PetscMin(t1,tt);
647: if (t1>0){
648: t2=PetscMin(t2,tt);
649: }
650: t3=PetscMax(t3,tt);
651: }
652: }
653: info=VecRestoreArray(X,&x);CHKERRQ(info);
654: info=VecRestoreArray(XL,&xl);CHKERRQ(info);
655: info=VecRestoreArray(XU,&xu);CHKERRQ(info);
656: info=VecRestoreArray(DX,&dx);CHKERRQ(info);
657: info=PetscObjectGetComm((PetscObject)X,&comm);CHKERRQ(info);
658:
659: if (boundmin){ info = MPI_Allreduce(&t1,&tt1,1,MPI_DOUBLE,MPI_MIN,comm);CHKERRQ(info);}
660: if (wolfemin){ info = MPI_Allreduce(&t2,&tt2,1,MPI_DOUBLE,MPI_MIN,comm);CHKERRQ(info);}
661: if (boundmax) { info = MPI_Allreduce(&t3,&tt3,1,MPI_DOUBLE,MPI_MAX,comm);CHKERRQ(info);}
663: *boundmin=(PetscReal)tt1;
664: *wolfemin=(PetscReal)tt2;
665: *boundmax=(PetscReal)tt3;
666: info = PetscInfo3(X,"Step Bound Info: Closest Bound: %6.4e, Wolfe: %6.4e, Max: %6.4e \n",*boundmin,*wolfemin,*boundmax); CHKERRQ(info);
667: return(0);
668: }
672: /*@C
673: VecISSetToConstant - Sets the elements of a vector, specified by an index set, to a constant
675: Input Parameter:
676: + S - a PETSc IS
677: . c - the constant
678: - V - a Vec
680: .seealso VecSet()
682: Level: advanced
683: @*/
684: int VecISSetToConstant(IS S, PetscScalar c, Vec V){
685: int info,i;
686: PetscInt nloc,low,high;
687: const PetscInt *s;
688: PetscScalar *v;
695: info = VecGetOwnershipRange(V, &low, &high); CHKERRQ(info);
696: info = ISGetLocalSize(S,&nloc);CHKERRQ(info);
698: info = ISGetIndices(S, &s); CHKERRQ(info);
699: info = VecGetArray(V,&v); CHKERRQ(info);
700: for (i=0; i<nloc; i++){
701: v[s[i]-low] = c;
702: /*
703: if (s[i]>=low && s[i]<high){
704: v[s[i]-low] = c;
705: } else {
706: SETERRQ(1,"IS index out of range");
707: }
708: */
709: }
710:
711: info = ISRestoreIndices(S, &s); CHKERRQ(info);
712: info = VecRestoreArray(V,&v); CHKERRQ(info);
714: return(0);
716: }
720: int VecPow(Vec Vec1, double p)
721: {
722: int i,info;
723: PetscInt n;
724: PetscScalar *v1;
729: info = VecGetArray(Vec1, &v1); CHKERRQ(info);
730: info = VecGetLocalSize(Vec1, &n); CHKERRQ(info);
732: if (1.0 == p) {
733: }
734: else if (-1.0 == p) {
735: for (i = 0; i < n; ++i){
736: v1[i] = 1.0 / v1[i];
737: }
738: }
739: else if (0.0 == p) {
740: for (i = 0; i < n; ++i){
741: // Not-a-number left alone
742: // Infinity set to one
743: if (v1[i] == v1[i]) {
744: v1[i] = 1.0;
745: }
746: }
747: }
748: else if (0.5 == p) {
749: for (i = 0; i < n; ++i) {
750: if (v1[i] >= 0) {
751: v1[i] = sqrt(v1[i]);
752: }
753: else {
754: v1[i] = TAO_INFINITY;
755: }
756: }
757: }
758: else if (-0.5 == p) {
759: for (i = 0; i < n; ++i) {
760: if (v1[i] >= 0) {
761: v1[i] = 1.0 / sqrt(v1[i]);
762: }
763: else {
764: v1[i] = TAO_INFINITY;
765: }
766: }
767: }
768: else if (2.0 == p) {
769: for (i = 0; i < n; ++i){
770: v1[i] *= v1[i];
771: }
772: }
773: else if (-2.0 == p) {
774: for (i = 0; i < n; ++i){
775: v1[i] = 1.0 / (v1[i] * v1[i]);
776: }
777: }
778: else {
779: for (i = 0; i < n; ++i) {
780: if (v1[i] >= 0) {
781: v1[i] = pow(v1[i], p);
782: }
783: else {
784: v1[i] = TAO_INFINITY;
785: }
786: }
787: }
789: info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
790: return(0);
791: }