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: