Actual source code: pvec.c
1: #include "tao_general.h"
2: #include "vecimpl.h" /*I "tao_solver.h" I*/
3: #include "petscksp.h"
6: int PetscMapSame(PetscMap m1, PetscMap m2, PetscTruth *flag){
7: int i,info;
8: int *grange1,*grange2;
9: int size1,size2;
10: MPI_Comm comm1,comm2;
13: *flag=PETSC_TRUE;
14: info = PetscMapGetGlobalRange(m1,&grange1);CHKERRQ(info);
15: info = PetscMapGetGlobalRange(m2,&grange2);CHKERRQ(info);
16: info = PetscObjectGetComm((PetscObject)m1,&comm1);CHKERRQ(info);
17: info = PetscObjectGetComm((PetscObject)m2,&comm2);CHKERRQ(info);
18: info = MPI_Comm_size(comm1,&size1); CHKERRQ(info);
19: info = MPI_Comm_size(comm2,&size2); CHKERRQ(info);
20: if (size1!=size2){
21: *flag=PETSC_FALSE;
22: } else {
23: for (i=0;i<size1+1;i++){
24: if (grange1[i]!=grange2[i] ){
25: *flag=PETSC_FALSE;
26: break;
27: }
28: }
29: }
30: return(0);
31: }
35: int VecCompare(Vec V1,Vec V2, PetscTruth *flg){
36: int info;
37: PetscMap m1,m2;
41: info=VecGetPetscMap(V1,&m1);CHKERRQ(info);
42: info=VecGetPetscMap(V2,&m2);CHKERRQ(info);
43: info=PetscMapSame(m1,m2,flg);CHKERRQ(info);
44: return(0);
45: }
47: /* ---------------------------------------------------------- */
50: int VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
51: {
52: int i,n,low1,low2,low3,low4,high1,high2,high3,high4,info;
53: PetscScalar *v1,*v2,*v3,*vmed;
61: if (Vec1==Vec2 || Vec1==Vec3){
62: info=VecCopy(Vec1,VMedian); CHKERRQ(info);
63: return(0);
64: }
65: if (Vec2==Vec3){
66: info=VecCopy(Vec2,VMedian); CHKERRQ(info);
67: return(0);
68: }
76: info = VecGetOwnershipRange(Vec1, &low1, &high1); CHKERRQ(info);
77: info = VecGetOwnershipRange(Vec2, &low2, &high2); CHKERRQ(info);
78: info = VecGetOwnershipRange(Vec3, &low3, &high3); CHKERRQ(info);
79: info = VecGetOwnershipRange(VMedian, &low4, &high4); CHKERRQ(info);
80: if ( low1!= low2 || low1!= low3 || low1!= low4 ||
81: high1!= high2 || high1!= high3 || high1!= high4)
82: SETERRQ(1,"InCompatible vector local lengths");
84: info = VecGetArray(Vec1,&v1); CHKERRQ(info);
85: info = VecGetArray(Vec2,&v2); CHKERRQ(info);
86: info = VecGetArray(Vec3,&v3); CHKERRQ(info);
88: if ( VMedian != Vec1 && VMedian != Vec2 && VMedian != Vec3){
89: info = VecGetArray(VMedian,&vmed); CHKERRQ(info);
90: } else if ( VMedian==Vec1 ){
91: vmed=v1;
92: } else if ( VMedian==Vec2 ){
93: vmed=v2;
94: } else {
95: vmed=v3;
96: }
98: info=VecGetLocalSize(Vec1,&n); CHKERRQ(info);
100: for (i=0;i<n;i++){
101: vmed[i]=PetscMax(PetscMax(PetscMin(v1[i],v2[i]),PetscMin(v1[i],v3[i])),PetscMin(v2[i],v3[i]));
102: }
104: info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
105: info = VecRestoreArray(Vec2,&v2); CHKERRQ(info);
106: info = VecRestoreArray(Vec3,&v2); CHKERRQ(info);
107:
108: if (VMedian!=Vec1 && VMedian != Vec2 && VMedian != Vec3){
109: info = VecRestoreArray(VMedian,&vmed); CHKERRQ(info);
110: }
112: return(0);
113: }
115: /* ----------------------------------------------------------
118: int VecPointwiseMin(Vec Vec1, Vec Vec2, Vec VMin)
119: {
120: int i,n,low1,low2,low3,high1,high2,high3,info;
121: PetscScalar *v1,*v2,*vmin;
128: if (Vec1==Vec2){
129: info=VecCopy(Vec1,VMin); CHKERRQ(info);
130: return(0);
131: }
136: info = VecGetOwnershipRange(Vec1, &low1, &high1); CHKERRQ(info);
137: info = VecGetOwnershipRange(Vec2, &low2, &high2); CHKERRQ(info);
138: info = VecGetOwnershipRange(VMin, &low3, &high3); CHKERRQ(info);
140: if ( low1!= low2 || low1!= low3 || high1!= high2 || high1!= high3)
141: SETERRQ(1,"InCompatible vector local lengths");
143: info = VecGetArray(Vec1,&v1); CHKERRQ(info);
144: info = VecGetArray(Vec2,&v2); CHKERRQ(info);
146: if ( VMin != Vec1 && VMin != Vec2){
147: info = VecGetArray(VMin,&vmin); CHKERRQ(info);
148: } else if ( VMin==Vec1 ){
149: vmin=v1;
150: } else {
151: vmin =v2;
152: }
154: info=VecGetLocalSize(Vec1,&n); CHKERRQ(info);
156: for (i=0; i<n; i++){
157: vmin[i] = PetscMin(v1[i],v2[i]);
158: }
160: info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
161: info = VecRestoreArray(Vec2,&v2); CHKERRQ(info);
163: if (VMin!=Vec1 && VMin != Vec2){
164: info = VecRestoreArray(VMin,&vmin); CHKERRQ(info);
165: }
167: return(0);
168: }
169: */
170: /* ----------------------------------------------------------
173: int VecPointwiseMax(Vec Vec1, Vec Vec2, Vec VMax)
174: {
175: int i,n,low1,low2,low3,high1,high2,high3,info;
176: PetscScalar *v1,*v2,*vmax;
180: if (Vec1==Vec2){
181: info=VecCopy(Vec1,VMax); CHKERRQ(info);
182: return(0);
183: }
192: info = VecGetOwnershipRange(Vec1, &low1, &high1); CHKERRQ(info);
193: info = VecGetOwnershipRange(Vec2, &low2, &high2); CHKERRQ(info);
194: info = VecGetOwnershipRange(VMax, &low3, &high3); CHKERRQ(info);
196: if ( low1!= low2 || low1!= low3 || high1!= high2 || high1!= high3)
197: SETERRQ(1,"Vectors must be identically loaded over processors");
199: info = VecGetArray(Vec1,&v1); CHKERRQ(info);
200: info = VecGetArray(Vec2,&v2); CHKERRQ(info);
202: if ( VMax != Vec1 && VMax != Vec2){
203: info = VecGetArray(VMax,&vmax); CHKERRQ(info);
204: } else if ( VMax==Vec1 ){ vmax=v1;
205: } else { vmax =v2;
206: }
208: info=VecGetLocalSize(Vec1,&n); CHKERRQ(info);
210: for (i=0; i<n; i++){
211: vmax[i] = PetscMax(v1[i],v2[i]);
212: }
214: info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
215: info = VecRestoreArray(Vec2,&v2); CHKERRQ(info);
217: if ( VMax != Vec1 && VMax != Vec2){
218: info = VecRestoreArray(VMax,&vmax); CHKERRQ(info);
219: }
221: return(0);
222: }*/
226: inline static PetscScalar Fischer(PetscScalar a, PetscScalar b)
227: {
229: #ifdef TODD
231: if (PetscAbsScalar(a) > PetscAbsScalar(b)) {
232: return sqrt(a*a + b*b) - a - b;
233: }
234: return sqrt(a*a + b*b) - b - a;
236: #else
238: // Method suggested by Bob Vanderbei
240: if (a + b <= 0) {
241: return sqrt(a*a + b*b) - (a + b);
242: }
243: return -2.0*a*b / (sqrt(a*a + b*b) + (a + b));
245: #endif
247: }
251: int VecFischer(Vec X, Vec F, Vec L, Vec U, Vec FF)
252: {
253: PetscScalar *x, *f, *l, *u, *ff;
254: PetscScalar xval, fval, lval, uval;
255: int info, i, n;
256: int low[5], high[5];
265: info = VecGetOwnershipRange(X, low, high); CHKERRQ(info);
266: info = VecGetOwnershipRange(F, low + 1, high + 1); CHKERRQ(info);
267: info = VecGetOwnershipRange(L, low + 2, high + 2); CHKERRQ(info);
268: info = VecGetOwnershipRange(U, low + 3, high + 3); CHKERRQ(info);
269: info = VecGetOwnershipRange(FF, low + 4, high + 4); CHKERRQ(info);
271: for (i = 1; i < 4; ++i) {
272: if (low[0] != low[i] || high[0] != high[i])
273: SETERRQ(1,"Vectors must be identically loaded over processors");
274: }
276: info = VecGetArray(X, &x); CHKERRQ(info);
277: info = VecGetArray(F, &f); CHKERRQ(info);
278: info = VecGetArray(L, &l); CHKERRQ(info);
279: info = VecGetArray(U, &u); CHKERRQ(info);
280: info = VecGetArray(FF, &ff); CHKERRQ(info);
282: info = VecGetLocalSize(X, &n); CHKERRQ(info);
284: for (i = 0; i < n; ++i) {
285: xval = (*x++); fval = (*f++);
286: lval = (*l++); uval = (*u++);
288: if ((lval <= -TAO_INFINITY) && (uval >= TAO_INFINITY)) {
289: (*ff++) = -fval;
290: }
291: else if (lval <= -TAO_INFINITY) {
292: (*ff++) = -Fischer(uval - xval, -fval);
293: }
294: else if (uval >= TAO_INFINITY) {
295: (*ff++) = Fischer(xval - lval, fval);
296: }
297: else if (lval == uval) {
298: (*ff++) = lval - xval;
299: }
300: else {
301: fval = Fischer(uval - xval, -fval);
302: (*ff++) = Fischer(xval - lval, fval);
303: }
305: }
306: x -= n; f -= n; l -=n; u -= n; ff -= n;
307:
308: info = VecRestoreArray(X, &x); CHKERRQ(info);
309: info = VecRestoreArray(F, &f); CHKERRQ(info);
310: info = VecRestoreArray(L, &l); CHKERRQ(info);
311: info = VecRestoreArray(U, &u); CHKERRQ(info);
312: info = VecRestoreArray(FF, &ff); CHKERRQ(info);
314: return(0);
315: }
319: inline static PetscScalar SFischer(PetscScalar a, PetscScalar b, PetscScalar c)
320: {
322: #ifdef TODD
324: if (PetscAbsScalar(a) > PetscAbsScalar(b)) {
325: return sqrt(a*a + b*b + 2.0*c*c) - a - b;
326: }
327: return sqrt(a*a + b*b + 2.0*c*c) - b - a;
329: #else
331: // Method suggested by Bob Vanderbei
333: if (a + b <= 0) {
334: return sqrt(a*a + b*b + 2.0*c*c) - (a + b);
335: }
336: return 2.0*(c*c - a*b) / (sqrt(a*a + b*b + 2.0*c*c) + (a + b));
338: #endif
340: }
344: int VecSFischer(Vec X, Vec F, Vec L, Vec U, double mu, Vec FF)
345: {
346: PetscScalar *x, *f, *l, *u, *ff;
347: PetscScalar xval, fval, lval, uval;
348: int info, i, n;
349: int low[5], high[5];
358: info = VecGetOwnershipRange(X, low, high); CHKERRQ(info);
359: info = VecGetOwnershipRange(F, low + 1, high + 1); CHKERRQ(info);
360: info = VecGetOwnershipRange(L, low + 2, high + 2); CHKERRQ(info);
361: info = VecGetOwnershipRange(U, low + 3, high + 3); CHKERRQ(info);
362: info = VecGetOwnershipRange(FF, low + 4, high + 4); CHKERRQ(info);
364: for (i = 1; i < 4; ++i) {
365: if (low[0] != low[i] || high[0] != high[i])
366: SETERRQ(1,"Vectors must be identically loaded over processors");
367: }
369: info = VecGetArray(X, &x); CHKERRQ(info);
370: info = VecGetArray(F, &f); CHKERRQ(info);
371: info = VecGetArray(L, &l); CHKERRQ(info);
372: info = VecGetArray(U, &u); CHKERRQ(info);
373: info = VecGetArray(FF, &ff); CHKERRQ(info);
375: info = VecGetLocalSize(X, &n); CHKERRQ(info);
377: for (i = 0; i < n; ++i) {
378: xval = (*x++); fval = (*f++);
379: lval = (*l++); uval = (*u++);
381: if ((lval <= -TAO_INFINITY) && (uval >= TAO_INFINITY)) {
382: (*ff++) = -fval - mu*xval;
383: }
384: else if (lval <= -TAO_INFINITY) {
385: (*ff++) = -SFischer(uval - xval, -fval, mu);
386: }
387: else if (uval >= TAO_INFINITY) {
388: (*ff++) = SFischer(xval - lval, fval, mu);
389: }
390: else if (lval == uval) {
391: (*ff++) = lval - xval;
392: }
393: else {
394: fval = SFischer(uval - xval, -fval, mu);
395: (*ff++) = SFischer(xval - lval, fval, mu);
396: }
397: }
398: x -= n; f -= n; l -=n; u -= n; ff -= n;
400: info = VecRestoreArray(X, &x); CHKERRQ(info);
401: info = VecRestoreArray(F, &f); CHKERRQ(info);
402: info = VecRestoreArray(L, &l); CHKERRQ(info);
403: info = VecRestoreArray(U, &u); CHKERRQ(info);
404: info = VecRestoreArray(FF, &ff); CHKERRQ(info);
405: return(0);
406: }
410: int VecBoundProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP){
412: int i,n,info;
413: PetscScalar *xptr,*xlptr,*xuptr,*gptr,*gpptr;
414: PetscScalar xval,gpval;
416: /* Project variables at the lower and upper bound */
419: info = VecGetLocalSize(X,&n); CHKERRQ(info);
421: info=VecGetArray(X,&xptr); CHKERRQ(info);
422: info=VecGetArray(XL,&xlptr); CHKERRQ(info);
423: info=VecGetArray(XU,&xuptr); CHKERRQ(info);
424: info=VecGetArray(G,&gptr); CHKERRQ(info);
425: if (G!=GP){
426: info=VecGetArray(GP,&gpptr); CHKERRQ(info);
427: } else { gpptr=gptr; }
429: for (i=0; i<n; ++i){
430: gpval = gptr[i]; xval = xptr[i];
432: if (gpval>0 && xval<=xlptr[i]){
433: gpval = 0;
434: } else if (gpval<0 && xval>=xuptr[i]){
435: gpval = 0;
436: }
437: gpptr[i] = gpval;
439: /*
440: if (xlptr[i] >= xuptr[i]){
441: gpptr[i]=0.0;
442: } else if (xptr[i] <= xlptr[i]+eps){
443: gpptr[i] = PetscMin(gptr[i],zero);
444: } else if (xptr[i] >= xuptr[i]-eps){
445: gpptr[i] = PetscMax(gptr[i],zero);
446: } else {
447: gpptr[i] = gptr[i];
448: }
449: */
450: }
452: info=VecRestoreArray(X,&xptr); CHKERRQ(info);
453: info=VecRestoreArray(XL,&xlptr); CHKERRQ(info);
454: info=VecRestoreArray(XU,&xuptr); CHKERRQ(info);
455: info=VecRestoreArray(G,&gptr); CHKERRQ(info);
456: if (G!=GP){
457: info=VecRestoreArray(GP,&gpptr); CHKERRQ(info);
458: }
459: return(0);
460: }
464: int VecISAXPY(Vec Y, PetscScalar alpha, Vec X, IS YIS){
466: int i,info,in1,in2,xlow,xhigh,ylow,yhigh,*yis;
467: PetscScalar *x,*y;
470: info=ISGetLocalSize(YIS,&in1); CHKERRQ(info);
471: info=VecGetLocalSize(X,&in2); CHKERRQ(info);
473: if ( in1 != in2 )
474: SETERRQ(1,"Index set and X vector must be identically loaded over processors");
476: info=VecGetOwnershipRange(X, &xlow, &xhigh); CHKERRQ(info);
477: info=VecGetOwnershipRange(Y, &ylow, &yhigh); CHKERRQ(info);
479: info=ISGetIndices(YIS, &yis); CHKERRQ(info);
481: info=VecGetArray(X,&x); CHKERRQ(info);
482: if (X != Y){
483: info=VecGetArray(Y,&y); CHKERRQ(info);
484: } else {
485: y=x;
486: }
488: for (i=0; i<in1; i++){
489: if (yis[i]>=ylow && yis[i]<yhigh && i+xlow < xhigh){
490: y[yis[i]-ylow] = y[yis[i]-ylow] + (alpha)*x[i];
491: } else {
492: SETERRQ(1,"IS index out of range");
493: }
494: }
495:
496: info=ISRestoreIndices(YIS, &yis); CHKERRQ(info);
498: info=VecRestoreArray(X,&x); CHKERRQ(info);
499: if (X != Y){
500: info=VecRestoreArray(Y,&y); CHKERRQ(info);
501: }
503: info = PetscLogFlops(2*in1); CHKERRQ(info);
505: return(0);
506: }
510: int VecCreateSubVec(Vec A,IS S,Vec *B)
511: {
512: int info,nloc,n;
513: int i,low,high,nnloc;
514: int *ptrS;
515: PetscScalar zero=0.0;
516: PetscScalar *ptrB,*ptrA;
517: VecType type_name;
518: MPI_Comm comm;
523: info=ISGetLocalSize(S,&nloc);CHKERRQ(info);
524: info=ISGetSize(S,&n);CHKERRQ(info);
525: info=PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(info);
526: info=VecCreate(comm,B);CHKERRQ(info);
527: info=VecSetSizes(*B,nloc,n);CHKERRQ(info);
528: info=VecGetType(A,&type_name);CHKERRQ(info);
529: info=VecSetType(*B,type_name);CHKERRQ(info);
530: info=VecSet(*B, zero);CHKERRQ(info);
531:
532: info=VecGetOwnershipRange(A,&low,&high);CHKERRQ(info);
533: info=VecGetLocalSize(A,&nnloc);CHKERRQ(info);
534: info=VecGetArray(A,&ptrA);CHKERRQ(info);
535: info=VecGetArray(*B,&ptrB);CHKERRQ(info);
536: info=ISGetIndices(S,&ptrS);CHKERRQ(info);
537: for (i=0;i<nloc;i++){ ptrB[i]=ptrA[ptrS[i]-low]; }
538: info=VecRestoreArray(A,&ptrA);CHKERRQ(info);
539: info=VecRestoreArray(*B,&ptrB);CHKERRQ(info);
540: info=ISRestoreIndices(S,&ptrS);CHKERRQ(info);
542: return(0);
543: }
547: int VecStepMax(Vec X, Vec DX, PetscReal*step){
548: int i,info,nn;
549: PetscScalar stepmax=1.0e300;
550: PetscScalar *xx,*dx;
551: double t1,t2;
554: info = VecGetLocalSize(X,&nn);CHKERRQ(info);
555: info = VecGetArray(X,&xx);CHKERRQ(info);
556: info = VecGetArray(DX,&dx);CHKERRQ(info);
557: for (i=0;i<nn;i++){
558: if (xx[i] < 0){
559: SETERRQ(1,"Vector must be positive");
560: } else if (dx[i]<0){ stepmax=PetscMin(stepmax,-xx[i]/dx[i]);
561: }
562: }
563: info = VecRestoreArray(X,&xx);CHKERRQ(info);
564: info = VecRestoreArray(DX,&dx);CHKERRQ(info);
565: t1=(double)stepmax;
566: info = MPI_Allreduce(&t1,&t2,1,MPI_DOUBLE,MPI_MIN,X->comm);
567: CHKERRQ(info);
568: *step=(PetscReal)t2;
569: return(0);
570: }
574: int MatCreateProjectionOperator(Vec A, IS S , Vec B, Mat *MM){
576: int m,n,M,N,low,high;
577: int info,i,*iptr;
578: PetscScalar one=1.0;
579: Mat MMM;
587: info = VecGetSize(B,&M);CHKERRQ(info);
588: info = VecGetLocalSize(B,&m);CHKERRQ(info);
589: info = VecGetSize(A,&N);CHKERRQ(info);
590: info = VecGetLocalSize(A,&n);CHKERRQ(info);
592: info = MatCreateMPIAIJ(A->comm,m,n,M,N,1,PETSC_NULL,1,PETSC_NULL,&MMM);
593: CHKERRQ(info);
594: /*
595: info = MatCreateMPIBDiag(pv->comm,m,M,N,M,1,PETSC_NULL,PETSC_NULL,&A);
596: */
597:
598: info = ISGetSize(S,&n);CHKERRQ(info);
599: if (n!=M){
600: SETERRQ(1,"Size of index set does not match the second vector.");
601: }
602:
603: info = ISGetLocalSize(S,&n);CHKERRQ(info);
604: if (n!=m){
605: SETERRQ(1,"Local size of index set does not match the second vector.");
606: }
607: info=VecGetOwnershipRange(B,&low,&high);CHKERRQ(info);
608: info = ISGetIndices(S,&iptr);CHKERRQ(info);
609: for (i=0; i<n; i++){
610: info=MatSetValue(MMM,low+i,iptr[i],one,INSERT_VALUES);CHKERRQ(info);
611: }
612: info = ISRestoreIndices(S,&iptr);CHKERRQ(info);
614: info = MatAssemblyBegin(MMM,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
615: info = MatAssemblyEnd(MMM,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
617: info = MatCompress(MMM);CHKERRQ(info);
619: *MM = MMM;
621: return(0);
622: }
626: int VecStepMax2(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *step2){
628: int i,nn1,nn2,info;
629: PetscScalar *xx,*dx,*xl,*xu;
630: PetscScalar stepmax2=0;
631: double t1,t2;
633: TaoFunctionBegin;
634: info = VecGetArray(X,&xx);CHKERRQ(info);
635: info = VecGetArray(XL,&xl);CHKERRQ(info);
636: info = VecGetArray(XU,&xu);CHKERRQ(info);
637: info = VecGetArray(DX,&dx);CHKERRQ(info);
638: info = VecGetLocalSize(X,&nn1);CHKERRQ(info);
639: info = VecGetLocalSize(DX,&nn2);CHKERRQ(info);
641: for (i=0;i<nn1;i++){
642: if (dx[i] > 0){
643: stepmax2=PetscMax(stepmax2,(xu[i]-xx[i])/dx[i]);
644: } else if (dx[i]<0){
645: stepmax2=PetscMax(stepmax2,(xl[i]-xx[i])/dx[i]);
646: }
647: }
648: info = VecRestoreArray(X,&xx);CHKERRQ(info);
649: info = VecRestoreArray(XL,&xl);CHKERRQ(info);
650: info = VecRestoreArray(XU,&xu);CHKERRQ(info);
651: info = VecRestoreArray(DX,&dx);CHKERRQ(info);
653: t1=(double)stepmax2;
654: info = MPI_Allreduce(&t1,&t2,1,MPI_DOUBLE,MPI_MAX,X->comm);
655: CHKERRQ(info);
656: *step2=(PetscReal)t2;
658: TaoFunctionReturn(0);
659: }
664: int VecBoundStepInfo(Vec X, Vec XL, Vec XU, Vec DX, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax){
665: int info,i,n;
666: PetscScalar *x,*xl,*xu,*dx;
667: double tt,t1=1.0e+20, t2=1.0e+20, t3=0;
668: double tt1,tt2,tt3;
669: MPI_Comm comm;
670:
672: info=VecGetArray(X,&x);CHKERRQ(info);
673: info=VecGetArray(XL,&xl);CHKERRQ(info);
674: info=VecGetArray(XU,&xu);CHKERRQ(info);
675: info=VecGetArray(DX,&dx);CHKERRQ(info);
676: info = VecGetLocalSize(X,&n);CHKERRQ(info);
677: for (i=0;i<n;i++){
678: if (dx[i]>0){
679: tt=(xu[i]-x[i])/dx[i];
680: t1=PetscMin(t1,tt);
681: if (t1>0){
682: t2=PetscMin(t2,tt);
683: }
684: t3=PetscMax(t3,tt);
685: } else if (dx[i]<0){
686: tt=(xl[i]-x[i])/dx[i];
687: t1=PetscMin(t1,tt);
688: if (t1>0){
689: t2=PetscMin(t2,tt);
690: }
691: t3=PetscMax(t3,tt);
692: }
693: }
694: info=VecRestoreArray(X,&x);CHKERRQ(info);
695: info=VecRestoreArray(XL,&xl);CHKERRQ(info);
696: info=VecRestoreArray(XU,&xu);CHKERRQ(info);
697: info=VecRestoreArray(DX,&dx);CHKERRQ(info);
698: info=PetscObjectGetComm((PetscObject)X,&comm);CHKERRQ(info);
699:
700: if (boundmin){ info = MPI_Allreduce(&t1,&tt1,1,MPI_DOUBLE,MPI_MIN,comm);CHKERRQ(info);}
701: if (wolfemin){ info = MPI_Allreduce(&t2,&tt2,1,MPI_DOUBLE,MPI_MIN,comm);CHKERRQ(info);}
702: if (boundmax) { info = MPI_Allreduce(&t3,&tt3,1,MPI_DOUBLE,MPI_MAX,comm);CHKERRQ(info);}
704: *boundmin=(PetscReal)tt1;
705: *wolfemin=(PetscReal)tt2;
706: *boundmax=(PetscReal)tt3;
707: info = PetscLogInfo((X,"Step Bound Info: Closest Bound: %6.4e, Wolfe: %6.4e, Max: %6.4e \n",*boundmin,*wolfemin,*boundmax)); CHKERRQ(info);
708: return(0);
709: }
713: /*@C
714: VecISSetToConstant - Sets the elements of a vector, specified by an index set, to a constant
716: Input Parameter:
717: + S - a PETSc IS
718: . c - the constant
719: - V - a Vec
721: .seealso VecSet()
723: Level: advanced
724: @*/
725: int VecISSetToConstant(IS S, PetscScalar c, Vec V){
726: int info,i,nloc,low,high,*s;
727: PetscScalar *v;
734: info = VecGetOwnershipRange(V, &low, &high); CHKERRQ(info);
735: info = ISGetLocalSize(S,&nloc);CHKERRQ(info);
737: info = ISGetIndices(S, &s); CHKERRQ(info);
738: info = VecGetArray(V,&v); CHKERRQ(info);
739: for (i=0; i<nloc; i++){
740: v[s[i]-low] = c;
741: /*
742: if (s[i]>=low && s[i]<high){
743: v[s[i]-low] = c;
744: } else {
745: SETERRQ(1,"IS index out of range");
746: }
747: */
748: }
749:
750: info = ISRestoreIndices(S, &s); CHKERRQ(info);
751: info = VecRestoreArray(V,&v); CHKERRQ(info);
753: return(0);
755: }