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