source: soft/giet_vm/applications/ocean/jacobcalc.c @ 806

Last change on this file since 806 was 799, checked in by alain, 9 years ago

Introducing a port of ocean application without macros.
File size: 10.8 KB
Line 
1/*************************************************************************/
2/*                                                                       */
3/*  Copyright (c) 1994 Stanford University                               */
4/*                                                                       */
5/*  All rights reserved.                                                 */
6/*                                                                       */
7/*  Permission is given to use, copy, and modify this software for any   */
8/*  non-commercial purpose as long as this copyright notice is not       */
9/*  removed.  All other uses, including redistribution in whole or in    */
10/*  part, are forbidden without prior written permission.                */
11/*                                                                       */
12/*  This software is provided with absolutely no warranty and no         */
13/*  support.                                                             */
14/*                                                                       */
15/*************************************************************************/
16
17/* Does the arakawa jacobian calculation (of the x and y matrices,
18   putting the results in the z matrix) for a subblock.  */
19
20#include "decs.h"
21
22////////////////////////////
23void jacobcalc( double ***x, 
24                double ***y, 
25                double ***z, 
26                long procid, 
27                long firstrow, 
28                long lastrow, 
29                long firstcol, 
30                long lastcol)
31{
32    double f1;
33    double f2;
34    double f3;
35    double f4;
36    double f5;
37    double f6;
38    double f7;
39    double f8;
40    long iindex;
41    long indexp1;
42    long indexm1;
43    long im1;
44    long ip1;
45    long i;
46    long j;
47    long jj;
48    double **t2a;
49    double **t2b;
50    double **t2c;
51    double *t1a;
52    double *t1b;
53    double *t1c;
54    double *t1d;
55    double *t1e;
56    double *t1f;
57    double *t1g;
58
59    t2a = (double **) z[procid];
60    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) 
61    {
62        t2a[0][0] = 0.0;
63    }
64    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) 
65    {
66        t2a[im - 1][0] = 0.0;
67    }
68    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) 
69    {
70        t2a[0][jm - 1] = 0.0;
71    }
72    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) 
73    {
74        t2a[im - 1][jm - 1] = 0.0;
75    }
76
77    t2a = (double **) x[procid];
78    jj = gps[procid]->neighbors[UPLEFT];
79    if (jj != -1) 
80    {
81        t2a[0][0] = x[jj][im - 2][jm - 2];
82    }
83    jj = gps[procid]->neighbors[UPRIGHT];
84    if (jj != -1) 
85    {
86        t2a[0][jm - 1] = x[jj][im - 2][1];
87    }
88    jj = gps[procid]->neighbors[DOWNLEFT];
89    if (jj != -1) 
90    {
91        t2a[im - 1][0] = x[jj][1][jm - 2];
92    }
93    jj = gps[procid]->neighbors[DOWNRIGHT];
94    if (jj != -1) 
95    {
96        t2a[im - 1][jm - 1] = x[jj][1][1];
97    }
98
99    t2a = (double **) y[procid];
100    jj = gps[procid]->neighbors[UPLEFT];
101    if (jj != -1) 
102    {
103        t2a[0][0] = y[jj][im - 2][jm - 2];
104    }
105    jj = gps[procid]->neighbors[UPRIGHT];
106    if (jj != -1) 
107    {
108        t2a[0][jm - 1] = y[jj][im - 2][1];
109    }
110    jj = gps[procid]->neighbors[DOWNLEFT];
111    if (jj != -1) 
112    {
113        t2a[im - 1][0] = y[jj][1][jm - 2];
114    }
115    jj = gps[procid]->neighbors[DOWNRIGHT];
116    if (jj != -1) 
117    {
118        t2a[im - 1][jm - 1] = y[jj][1][1];
119    }
120
121    t2a = (double **) x[procid];
122    if (gps[procid]->neighbors[UP] == -1) 
123    {
124        jj = gps[procid]->neighbors[LEFT];
125        if (jj != -1) 
126        {
127            t2a[0][0] = x[jj][0][jm - 2];
128        } 
129        else 
130        {
131            jj = gps[procid]->neighbors[DOWN];
132            if (jj != -1) 
133            {
134                t2a[im - 1][0] = x[jj][1][0];
135            }
136        }
137        jj = gps[procid]->neighbors[RIGHT];
138        if (jj != -1) 
139        {
140            t2a[0][jm - 1] = x[jj][0][1];
141        } 
142        else 
143        {
144            jj = gps[procid]->neighbors[DOWN];
145            if (jj != -1) 
146            {
147                t2a[im - 1][jm - 1] = x[jj][1][jm - 1];
148            }
149        }
150    } 
151    else if (gps[procid]->neighbors[DOWN] == -1) 
152    {
153        jj = gps[procid]->neighbors[LEFT];
154        if (jj != -1) {
155            t2a[im - 1][0] = x[jj][im - 1][jm - 2];
156        } else {
157            jj = gps[procid]->neighbors[UP];
158            if (jj != -1) {
159                t2a[0][0] = x[jj][im - 2][0];
160            }
161        }
162        jj = gps[procid]->neighbors[RIGHT];
163        if (jj != -1) {
164            t2a[im - 1][jm - 1] = x[jj][im - 1][1];
165        } else {
166            jj = gps[procid]->neighbors[UP];
167            if (jj != -1) {
168                t2a[0][jm - 1] = x[jj][im - 2][jm - 1];
169            }
170        }
171    } else if (gps[procid]->neighbors[LEFT] == -1) {
172        jj = gps[procid]->neighbors[UP];
173        if (jj != -1) {
174            t2a[0][0] = x[jj][im - 2][0];
175        }
176        jj = gps[procid]->neighbors[DOWN];
177        if (jj != -1) {
178            t2a[im - 1][0] = x[jj][1][0];
179        }
180    } 
181    else if (gps[procid]->neighbors[RIGHT] == -1) 
182    {
183        jj = gps[procid]->neighbors[UP];
184        if (jj != -1) {
185            t2a[0][jm - 1] = x[jj][im - 2][jm - 1];
186        }
187        jj = gps[procid]->neighbors[DOWN];
188        if (jj != -1) {
189            t2a[im - 1][jm - 1] = x[jj][1][jm - 1];
190        }
191    }
192
193    t2a = (double **) y[procid];
194    if (gps[procid]->neighbors[UP] == -1) {
195        jj = gps[procid]->neighbors[LEFT];
196        if (jj != -1) {
197            t2a[0][0] = y[jj][0][jm - 2];
198        } else {
199            jj = gps[procid]->neighbors[DOWN];
200            if (jj != -1) {
201                t2a[im - 1][0] = y[jj][1][0];
202            }
203        }
204        jj = gps[procid]->neighbors[RIGHT];
205        if (jj != -1) {
206            t2a[0][jm - 1] = y[jj][0][1];
207        } else {
208            jj = gps[procid]->neighbors[DOWN];
209            if (jj != -1) {
210                t2a[im - 1][jm - 1] = y[jj][1][jm - 1];
211            }
212        }
213    } else if (gps[procid]->neighbors[DOWN] == -1) {
214        jj = gps[procid]->neighbors[LEFT];
215        if (jj != -1) {
216            t2a[im - 1][0] = y[jj][im - 1][jm - 2];
217        } else {
218            jj = gps[procid]->neighbors[UP];
219            if (jj != -1) {
220                t2a[0][0] = y[jj][im - 2][0];
221            }
222        }
223        jj = gps[procid]->neighbors[RIGHT];
224        if (jj != -1) {
225            t2a[im - 1][jm - 1] = y[jj][im - 1][1];
226        } else {
227            jj = gps[procid]->neighbors[UP];
228            if (jj != -1) {
229                t2a[0][jm - 1] = y[jj][im - 2][jm - 1];
230            }
231        }
232    } else if (gps[procid]->neighbors[LEFT] == -1) {
233        jj = gps[procid]->neighbors[UP];
234        if (jj != -1) {
235            t2a[0][0] = y[jj][im - 2][0];
236        }
237        jj = gps[procid]->neighbors[DOWN];
238        if (jj != -1) {
239            t2a[im - 1][0] = y[jj][1][0];
240        }
241    } else if (gps[procid]->neighbors[RIGHT] == -1) {
242        jj = gps[procid]->neighbors[UP];
243        if (jj != -1) {
244            t2a[0][jm - 1] = y[jj][im - 2][jm - 1];
245        }
246        jj = gps[procid]->neighbors[DOWN];
247        if (jj != -1) {
248            t2a[im - 1][jm - 1] = y[jj][1][jm - 1];
249        }
250    }
251
252    j = gps[procid]->neighbors[UP];
253    if (j != -1) {
254        t1a = (double *) t2a[0];
255        t1b = (double *) y[j][im - 2];
256        for (i = 1; i <= lastcol; i++) {
257            t1a[i] = t1b[i];
258        }
259    }
260    j = gps[procid]->neighbors[DOWN];
261    if (j != -1) {
262        t1a = (double *) t2a[im - 1];
263        t1b = (double *) y[j][1];
264        for (i = 1; i <= lastcol; i++) {
265            t1a[i] = t1b[i];
266        }
267    }
268    j = gps[procid]->neighbors[LEFT];
269    if (j != -1) {
270        t2b = (double **) y[j];
271        for (i = 1; i <= lastrow; i++) {
272            t2a[i][0] = t2b[i][jm - 2];
273        }
274    }
275    j = gps[procid]->neighbors[RIGHT];
276    if (j != -1) {
277        t2b = (double **) y[j];
278        for (i = 1; i <= lastrow; i++) {
279            t2a[i][jm - 1] = t2b[i][1];
280        }
281    }
282
283    t2a = (double **) x[procid];
284    j = gps[procid]->neighbors[UP];
285    if (j != -1) {
286        t1a = (double *) t2a[0];
287        t1b = (double *) x[j][im - 2];
288        for (i = 1; i <= lastcol; i++) {
289            t1a[i] = t1b[i];
290        }
291    }
292    j = gps[procid]->neighbors[DOWN];
293    if (j != -1) {
294        t1a = (double *) t2a[im - 1];
295        t1b = (double *) x[j][1];
296        for (i = 1; i <= lastcol; i++) {
297            t1a[i] = t1b[i];
298        }
299    }
300    j = gps[procid]->neighbors[LEFT];
301    if (j != -1) {
302        t2b = (double **) x[j];
303        for (i = 1; i <= lastrow; i++) {
304            t2a[i][0] = t2b[i][jm - 2];
305        }
306    }
307    j = gps[procid]->neighbors[RIGHT];
308    if (j != -1) {
309        t2b = (double **) x[j];
310        for (i = 1; i <= lastrow; i++) {
311            t2a[i][jm - 1] = t2b[i][1];
312        }
313    }
314
315    t2a = (double **) x[procid];
316    t2b = (double **) y[procid];
317    t2c = (double **) z[procid];
318    for (i = firstrow; i <= lastrow; i++) {
319        ip1 = i + 1;
320        im1 = i - 1;
321        t1a = (double *) t2a[i];
322        t1b = (double *) t2b[i];
323        t1c = (double *) t2c[i];
324        t1d = (double *) t2b[ip1];
325        t1e = (double *) t2b[im1];
326        t1f = (double *) t2a[ip1];
327        t1g = (double *) t2a[im1];
328        for (iindex = firstcol; iindex <= lastcol; iindex++) {
329            indexp1 = iindex + 1;
330            indexm1 = iindex - 1;
331            f1 = (t1b[indexm1] + t1d[indexm1] - t1b[indexp1] - t1d[indexp1]) * (t1f[iindex] - t1a[iindex]);
332            f2 = (t1e[indexm1] + t1b[indexm1] - t1e[indexp1] - t1b[indexp1]) * (t1a[iindex] - t1g[iindex]);
333            f3 = (t1d[iindex] + t1d[indexp1] - t1e[iindex] - t1e[indexp1]) * (t1a[indexp1] - t1a[iindex]);
334            f4 = (t1d[indexm1] + t1d[iindex] - t1e[indexm1] - t1e[iindex]) * (t1a[iindex] - t1a[indexm1]);
335            f5 = (t1d[iindex] - t1b[indexp1]) * (t1f[indexp1] - t1a[iindex]);
336            f6 = (t1b[indexm1] - t1e[iindex]) * (t1a[iindex] - t1g[indexm1]);
337            f7 = (t1b[indexp1] - t1e[iindex]) * (t1g[indexp1] - t1a[iindex]);
338            f8 = (t1d[iindex] - t1b[indexm1]) * (t1a[iindex] - t1f[indexm1]);
339
340            t1c[iindex] = factjacob * (f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8);
341        }
342    }
343
344    if (gps[procid]->neighbors[UP] == -1) 
345    {
346        t1c = (double *) t2c[0];
347        for (j = firstcol; j <= lastcol; j++) 
348        {
349            t1c[j] = 0.0;
350        }
351    }
352    if (gps[procid]->neighbors[DOWN] == -1) 
353    {
354        t1c = (double *) t2c[im - 1];
355        for (j = firstcol; j <= lastcol; j++) 
356        {
357            t1c[j] = 0.0;
358        }
359    }
360    if (gps[procid]->neighbors[LEFT] == -1) 
361    {
362        for (j = firstrow; j <= lastrow; j++) 
363        {
364            t2c[j][0] = 0.0;
365        }
366    }
367    if (gps[procid]->neighbors[RIGHT] == -1) 
368    {
369        for (j = firstrow; j <= lastrow; j++) 
370        {
371            t2c[j][jm - 1] = 0.0;
372        }
373    }
374
375}
Note: See TracBrowser for help on using the repository browser.