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