Actual source code: projection.c
1: #include "taosolver.h"
5: /*@C
6: VecBoundGradientProjection - Projects vector according to this definition.
7: If XL[i] < X[i] < XU[i], then GP[i] = G[i];
8: If X[i]<=XL[i], then GP[i] = min(G[i],0);
9: If X[i]>=XU[i], then GP[i] = max(G[i],0);
11: Input Parameters:
12: + G - current gradient vector
13: . X - current solution vector
14: . XL - lower bounds
15: - XU - upper bounds
17: Output Parameter:
18: . GP - gradient projection vector
20: Level: advanced
21: @*/
22: PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP){
25: PetscInt n,i;
26: PetscReal *xptr,*xlptr,*xuptr,*gptr,*gpptr;
27: PetscReal xval,gpval;
29: /* Project variables at the lower and upper bound */
38: VecGetLocalSize(X,&n);
40: ierr=VecGetArray(X,&xptr);
41: ierr=VecGetArray(XL,&xlptr);
42: ierr=VecGetArray(XU,&xuptr);
43: ierr=VecGetArray(G,&gptr);
44: if (G!=GP){
45: ierr=VecGetArray(GP,&gpptr);
46: } else { gpptr=gptr; }
48: for (i=0; i<n; ++i){
49: gpval = gptr[i]; xval = xptr[i];
51: if (gpval>0 && xval<=xlptr[i]){
52: gpval = 0;
53: } else if (gpval<0 && xval>=xuptr[i]){
54: gpval = 0;
55: }
56: gpptr[i] = gpval;
57: }
59: ierr=VecRestoreArray(X,&xptr);
60: ierr=VecRestoreArray(XL,&xlptr);
61: ierr=VecRestoreArray(XU,&xuptr);
62: ierr=VecRestoreArray(G,&gptr);
63: if (G!=GP){
64: ierr=VecRestoreArray(GP,&gpptr);
65: }
66: return(0);
67: }
71: PetscErrorCode VecStepMaxBounded(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *stepmax){
74: PetscInt i,nn;
75: PetscReal *xx,*dx,*xl,*xu;
76: PetscReal localmax=0;
84: VecGetArray(X,&xx);
85: VecGetArray(XL,&xl);
86: VecGetArray(XU,&xu);
87: VecGetArray(DX,&dx);
88: VecGetLocalSize(X,&nn);
90: for (i=0;i<nn;i++){
91: if (dx[i] > 0){
92: localmax=PetscMax(localmax,(xu[i]-xx[i])/dx[i]);
93: } else if (dx[i]<0){
94: localmax=PetscMax(localmax,(xl[i]-xx[i])/dx[i]);
95: }
96: }
97: VecRestoreArray(X,&xx);
98: VecRestoreArray(XL,&xl);
99: VecRestoreArray(XU,&xu);
100: VecRestoreArray(DX,&dx);
102: MPI_Allreduce(&localmax,stepmax,1,MPIU_REAL,MPIU_MAX,((PetscObject)X)->comm);
103:
106: return(0);
107: }
110: PetscErrorCode VecStepBoundInfo(Vec X, Vec XL, Vec XU, Vec DX, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax){
113: PetscInt n,i;
114: PetscReal *x,*xl,*xu,*dx;
115: PetscReal t;
116: PetscReal localmin=1.0e300,localwolfemin=1.0e300,localmax=0;
117: MPI_Comm comm;
118:
125: ierr=VecGetArray(X,&x);
126: ierr=VecGetArray(XL,&xl);
127: ierr=VecGetArray(XU,&xu);
128: ierr=VecGetArray(DX,&dx);
129: VecGetLocalSize(X,&n);
130: for (i=0;i<n;i++){
131: if (dx[i]>0){
132: t=(xu[i]-x[i])/dx[i];
133: localmin=PetscMin(t,localmin);
134: if (localmin>0){
135: localwolfemin = PetscMin(t,localwolfemin);
136: }
137: localmax = PetscMax(t,localmax);
138: } else if (dx[i]<0){
139: t=(xl[i]-x[i])/dx[i];
140: localmin = PetscMin(t,localmin);
141: if (localmin>0){
142: localwolfemin = PetscMin(t,localwolfemin);
143: }
144: localmax = PetscMax(t,localmax);
145: }
146: }
147: ierr=VecRestoreArray(X,&x);
148: ierr=VecRestoreArray(XL,&xl);
149: ierr=VecRestoreArray(XU,&xu);
150: ierr=VecRestoreArray(DX,&dx);
151: ierr=PetscObjectGetComm((PetscObject)X,&comm);
152:
153: if (boundmin){ MPI_Allreduce(&localmin,boundmin,1,MPIU_REAL,MPIU_MIN,comm);}
154: if (wolfemin){ MPI_Allreduce(&localwolfemin,wolfemin,1,MPIU_REAL,MPIU_MIN,comm);}
155: if (boundmax) { MPI_Allreduce(&localmax,boundmax,1,MPIU_REAL,MPIU_MAX,comm);}
157: PetscInfo3(X,"Step Bound Info: Closest Bound: %G, Wolfe: %G, Max: %G \n",*boundmin,*wolfemin,*boundmax);
158: return(0);
159: }
164: PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step){
166: PetscInt i, nn;
167: PetscReal stepmax=TAO_INFINITY;
168: PetscReal *xx, *dx;
169:
174: VecGetLocalSize(X,&nn);
175: VecGetArray(X,&xx);
176: VecGetArray(DX,&dx);
177: for (i=0;i<nn;i++){
178: if (xx[i] < 0){
179: SETERRQ(PETSC_COMM_SELF,1,"Vector must be positive");
180: } else if (dx[i]<0){ stepmax=PetscMin(stepmax,-xx[i]/dx[i]);
181: }
182: }
183: VecRestoreArray(X,&xx);
184: VecRestoreArray(DX,&dx);
185: MPI_Allreduce(&stepmax,step,1,MPIU_REAL,MPIU_MIN,((PetscObject)X)->comm);
186:
187: return(0);
188: }
189:
190: