Actual source code: pmat.c

  1: #include "tao_general.h"

  3: #ifdef TAO_USE_PETSC
  4: #include "petscmat.h"

  6: //#include "src/mat/matimpl.h"

  8: int MatD_Fischer(Mat, Vec, Vec, Vec, Vec, Vec, Vec, Vec, Vec);
  9: int MatD_SFischer(Mat, Vec, Vec, Vec, Vec, double, Vec, Vec, Vec, Vec, Vec);

 13: inline static PetscScalar Fischer(PetscScalar a, PetscScalar b)
 14: {

 16: #ifdef TODD

 18:   if (PetscAbsScalar(a) > PetscAbsScalar(b)) {
 19:     return sqrt(a*a + b*b) - a - b;
 20:   }
 21:   return sqrt(a*a + b*b) - b - a;

 23: #else

 25:    // Method suggested by Bob Vanderbei

 27:    if (a + b <= 0) {
 28:      return sqrt(a*a + b*b) - (a + b);
 29:    }
 30:    return -2.0*a*b / (sqrt(a*a + b*b) + (a + b));

 32: #endif

 34: }

 38: inline static PetscScalar Norm(PetscScalar a, PetscScalar b)
 39: {
 40:   return sqrt(a*a + b*b);
 41: }

 45: int MatD_Fischer(Mat m, Vec X, Vec F, Vec L, Vec U,
 46:                  Vec T1, Vec T2, Vec DA, Vec DB)
 47: {
 48:   PetscScalar *x, *f, *l, *u, *da, *db, *t1, *t2;
 49:   PetscScalar ai, bi, ci, di, ei;
 50:   int info, i, n;
 51:   int low[8], high[8];


 63:   info = VecGetOwnershipRange(X, low, high); CHKERRQ(info);
 64:   info = VecGetOwnershipRange(F, low + 1, high + 1); CHKERRQ(info);
 65:   info = VecGetOwnershipRange(L, low + 2, high + 2); CHKERRQ(info);
 66:   info = VecGetOwnershipRange(U, low + 3, high + 3); CHKERRQ(info);
 67:   info = VecGetOwnershipRange(DA, low + 4, high + 4); CHKERRQ(info);
 68:   info = VecGetOwnershipRange(DB, low + 5, high + 5); CHKERRQ(info);
 69:   info = VecGetOwnershipRange(T1, low + 6, high + 6); CHKERRQ(info);
 70:   info = VecGetOwnershipRange(T2, low + 7, high + 7); CHKERRQ(info);

 72:   for (i = 1; i < 8; i++) {
 73:     if (low[0] != low[i] || high[0] != high[i])
 74:       SETERRQ(1,"Vectors must be identically loaded over processors");
 75:   }

 77:   info = VecGetArray(X, &x); CHKERRQ(info);
 78:   info = VecGetArray(F, &f); CHKERRQ(info);
 79:   info = VecGetArray(L, &l); CHKERRQ(info);
 80:   info = VecGetArray(U, &u); CHKERRQ(info);
 81:   info = VecGetArray(DA, &da); CHKERRQ(info);
 82:   info = VecGetArray(DB, &db); CHKERRQ(info);
 83:   info = VecGetArray(T1, &t1); CHKERRQ(info);

 85:   info = VecGetLocalSize(X, &n); CHKERRQ(info);

 87:   for (i = 0; i < n; i++) {
 88:     da[i] = 0;
 89:     db[i] = 0;
 90:     t1[i] = 0;

 92:     if (PetscAbsScalar(f[i]) <= TAO_EPSILON) {
 93:       if (l[i] > -TAO_INFINITY && PetscAbsScalar(x[i] - l[i]) <= TAO_EPSILON) {
 94:         t1[i] = 1;
 95:         da[i] = 1;
 96:       }

 98:       if (u[i] <  TAO_INFINITY && PetscAbsScalar(u[i] - x[i]) <= TAO_EPSILON) {
 99:         t1[i] = 1;
100:         db[i] = 1;
101:       }
102:     }
103:   }

105:   info = VecRestoreArray(T1, &t1); CHKERRQ(info);
106:   info = MatMult(m, T1, T2); CHKERRQ(info);
107:   info = VecGetArray(T2, &t2); CHKERRQ(info);
108: 
109:   for (i = 0; i < n; i++) {
110:     if ((l[i] <= -TAO_INFINITY) && (u[i] >= TAO_INFINITY)) {
111:       da[i] = 0;
112:       db[i] = -1;
113:     }
114:     else if (l[i] <= -TAO_INFINITY) {
115:       if (db[i] >= 1) {
116:         ai = Norm(1, t2[i]);

118:         da[i] = -1/ai - 1;
119:         db[i] = -t2[i]/ai - 1;
120:       }
121:       else {
122:         bi = u[i] - x[i];
123:         ai = Norm(bi, f[i]);
124:         ai = PetscMax(TAO_EPSILON, ai);

126:         da[i] = bi / ai - 1;
127:         db[i] = -f[i] / ai - 1;
128:       }
129:     }
130:     else if (u[i] >=  TAO_INFINITY) {
131:       if (da[i] >= 1) {
132:         ai = Norm(1, t2[i]);

134:         da[i] = 1 / ai - 1;
135:         db[i] = t2[i] / ai - 1;
136:       }
137:       else {
138:         bi = x[i] - l[i];
139:         ai = Norm(bi, f[i]);
140:         ai = PetscMax(TAO_EPSILON, ai);

142:         da[i] = bi / ai - 1;
143:         db[i] = f[i] / ai - 1;
144:       }
145:     }
146:     else if (l[i] == u[i]) {
147:       da[i] = -1;
148:       db[i] = 0;
149:     }
150:     else {
151:       if (db[i] >= 1) {
152:         ai = Norm(1, t2[i]);

154:         ci = 1 / ai + 1;
155:         di = t2[i] / ai + 1;
156:       }
157:       else {
158:         bi = x[i] - u[i];
159:         ai = Norm(bi, f[i]);
160:         ai = PetscMax(TAO_EPSILON, ai);

162:         ci = bi / ai + 1;
163:         di = f[i] / ai + 1;
164:       }

166:       if (da[i] >= 1) {
167:         bi = ci + di*t2[i];
168:         ai = Norm(1, bi);
169: 
170:         bi = bi / ai - 1;
171:         ai = 1 / ai - 1;
172:       }
173:       else {
174:         ei = Fischer(u[i] - x[i], -f[i]);
175:         ai = Norm(x[i] - l[i], ei);
176:         ai = PetscMax(TAO_EPSILON, ai);

178:         bi = ei / ai - 1;
179:         ai = (x[i] - l[i]) / ai - 1;
180:       }

182:       da[i] = ai + bi*ci;
183:       db[i] = bi*di;
184:     }
185:   }

187:   info = VecRestoreArray(DA, &da); CHKERRQ(info);
188:   info = VecRestoreArray(DB, &db); CHKERRQ(info);

190:   info = VecRestoreArray(X, &x); CHKERRQ(info);
191:   info = VecRestoreArray(F, &f); CHKERRQ(info);
192:   info = VecRestoreArray(L, &l); CHKERRQ(info);
193:   info = VecRestoreArray(U, &u); CHKERRQ(info);
194:   info = VecRestoreArray(T2, &t2); CHKERRQ(info);
195:   return(0);
196: }

200: inline static PetscScalar SFischer(PetscScalar a, PetscScalar b, PetscScalar c)
201: {

203: #ifdef TODD

205:   if (PetscAbsScalar(a) > PetscAbsScalar(b)) {
206:     return sqrt(a*a + b*b + 2.0*c*c) - a - b;
207:   }
208:   return sqrt(a*a + b*b + 2.0*c*c) - b - a;

210: #else

212:    // Method suggested by Bob Vanderbei

214:    if (a + b <= 0) {
215:      return sqrt(a*a + b*b + 2.0*c*c) - (a + b);
216:    }
217:    return 2.0*(c*c - a*b) / (sqrt(a*a + b*b + 2.0*c*c) + (a + b));

219: #endif

221: }

225: inline static PetscScalar SNorm(PetscScalar a, PetscScalar b, PetscScalar c)
226: {
227:   return sqrt(a*a + b*b + 2.0*c*c);
228: }

232: int MatD_SFischer(Mat m, Vec X, Vec F, Vec L, Vec U, double mu,
233:                   Vec T1, Vec T2, Vec DA, Vec DB, Vec DM)
234: {
235:   PetscScalar *x, *f, *l, *u, *da, *db, *dm;
236:   PetscScalar ai, bi, ci, di, ei, fi;
237:   int info, i, n;
238:   int low[7], high[7];


249:   info = VecGetOwnershipRange(X, low, high); CHKERRQ(info);
250:   info = VecGetOwnershipRange(F, low + 1, high + 1); CHKERRQ(info);
251:   info = VecGetOwnershipRange(L, low + 2, high + 2); CHKERRQ(info);
252:   info = VecGetOwnershipRange(U, low + 3, high + 3); CHKERRQ(info);
253:   info = VecGetOwnershipRange(DA, low + 4, high + 4); CHKERRQ(info);
254:   info = VecGetOwnershipRange(DB, low + 5, high + 5); CHKERRQ(info);
255:   info = VecGetOwnershipRange(DM, low + 6, high + 6); CHKERRQ(info);

257:   for (i = 1; i < 7; i++) {
258:     if (low[0] != low[i] || high[0] != high[i])
259:       SETERRQ(1,"Vectors must be identically loaded over processors");
260:   }

262:   info = VecGetArray(X, &x); CHKERRQ(info);
263:   info = VecGetArray(F, &f); CHKERRQ(info);
264:   info = VecGetArray(L, &l); CHKERRQ(info);
265:   info = VecGetArray(U, &u); CHKERRQ(info);
266:   info = VecGetArray(DA, &da); CHKERRQ(info);
267:   info = VecGetArray(DB, &db); CHKERRQ(info);
268:   info = VecGetArray(DM, &dm); CHKERRQ(info);

270:   info = VecGetLocalSize(X, &n); CHKERRQ(info);

272:   for (i = 0; i < n; i++) {
273:     if ((l[i] <= -TAO_INFINITY) && (u[i] >= TAO_INFINITY)) {
274:       da[i] = -mu;
275:       db[i] = -1;
276:       dm[i] = -x[i];
277:     }
278:     else if (l[i] <= -TAO_INFINITY) {
279:       bi = u[i] - x[i];
280:       ai = SNorm(bi, f[i], mu);
281:       ai = PetscMax(TAO_EPSILON, ai);

283:       da[i] = bi / ai - 1;
284:       db[i] = -f[i] / ai - 1;
285:       dm[i] = 2.0 * mu / ai;
286:     }
287:     else if (u[i] >=  TAO_INFINITY) {
288:       bi = x[i] - l[i];
289:       ai = SNorm(bi, f[i], mu);
290:       ai = PetscMax(TAO_EPSILON, ai);

292:       da[i] = bi / ai - 1;
293:       db[i] = f[i] / ai - 1;
294:       dm[i] = 2.0 * mu / ai;
295:     }
296:     else if (l[i] == u[i]) {
297:       da[i] = -1;
298:       db[i] = 0;
299:       dm[i] = 0;
300:     }
301:     else {
302:       bi = x[i] - u[i];
303:       ai = SNorm(bi, f[i], mu);
304:       ai = PetscMax(TAO_EPSILON, ai);

306:       ci = bi / ai + 1;
307:       di = f[i] / ai + 1;
308:       fi = 2.0 * mu / ai;

310:       ei = SFischer(u[i] - x[i], -f[i], mu);
311:       ai = SNorm(x[i] - l[i], ei, mu);
312:       ai = PetscMax(TAO_EPSILON, ai);

314:       bi = ei / ai - 1;
315:       ai = (x[i] - l[i]) / ai - 1;

317:       da[i] = ai + bi*ci;
318:       db[i] = bi*di;
319:       dm[i] = ei + bi*fi;
320:     }
321:   }

323:   info = VecRestoreArray(X, &x); CHKERRQ(info);
324:   info = VecRestoreArray(F, &f); CHKERRQ(info);
325:   info = VecRestoreArray(L, &l); CHKERRQ(info);
326:   info = VecRestoreArray(U, &u); CHKERRQ(info);
327:   info = VecRestoreArray(DA, &da); CHKERRQ(info);
328:   info = VecRestoreArray(DB, &db); CHKERRQ(info);
329:   info = VecRestoreArray(DM, &dm); CHKERRQ(info);
330:   return(0);
331: }

333: #endif