source: soft/giet_vm/applications/ocean/slave2.c @ 808

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

Introducing a port of ocean application without macros.
File size: 42.2 KB
RevLine 
[799]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/*    ****************
18      subroutine slave2
19      ****************  */
20
21#include <user_barrier.h>
22#include <user_lock.h>
23#include <math.h>
24#include <stdio.h>
25
26#include "decs.h"
27
28///////////////////////////
29void slave2(long procid, 
30            long firstrow, 
31            long lastrow, 
32            long numrows, 
33            long firstcol, 
34            long lastcol, 
35            long numcols)
36{
37    long       i;
38    long       j;
39    long       iindex;
40    double     hh1;
41    double     hh3;
42    double     hinv;
43    double     h1inv;
44    long       istart;
45    long       iend;
46    long       jstart;
47    long       jend;
48    long       ist;
49    long       ien;
50    long       jst;
51    long       jen;
52    double     ressqr;
53    double     psiaipriv;
54    double     f4;
55    double     timst;
56    long       psiindex;
57    long       i_off;
58    long       j_off;
59    long       multi_start = 0;
60    double **  t2a;
61    double **  t2b;
62    double **  t2c;
63    double **  t2d;
64    double **  t2e;
65    double **  t2f;
66    double **  t2g;
67    double **  t2h;
68    double *   t1a;
69    double *   t1b;
70    double *   t1c;
71    double *   t1d;
72    double *   t1e;
73    double *   t1f;
74    double *   t1g;
75    double *   t1h;
76    long       barrier_start;
77    long       step_start;
78
79    ressqr = lev_res[numlev - 1] * lev_res[numlev - 1];
80    i_off  = (gps[procid]->rownum) * numrows;
81    j_off  = (gps[procid]->colnum) * numcols;
82
83
84/*   ***************************************************************
85          f i r s t     p h a s e   (of timestep calculation)
86     ***************************************************************/
87
88    // BEGIN
89    step_start = giet_proctime();
90
91    t2a = (double **) ga[procid];
92    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) 
93    {
94        t2a[0][0] = 0.0;
95    }
96    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) 
97    {
98        t2a[im - 1][0] = 0.0;
99    }
100    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) 
101    {
102        t2a[0][jm - 1] = 0.0;
103    }
104    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) 
105    {
106        t2a[im - 1][jm - 1] = 0.0;
107    }
108    if (gps[procid]->neighbors[UP] == -1) 
109    {
110        t1a = (double *) t2a[0];
111        for (j = firstcol; j <= lastcol; j++) 
112        {
113            t1a[j] = 0.0;
114        }
115    }
116    if (gps[procid]->neighbors[DOWN] == -1) 
117    {
118        t1a = (double *) t2a[im - 1];
119        for (j = firstcol; j <= lastcol; j++) 
120        {
121            t1a[j] = 0.0;
122        }
123    }
124    if (gps[procid]->neighbors[LEFT] == -1) 
125    {
126        for (j = firstrow; j <= lastrow; j++) 
127        {
128            t2a[j][0] = 0.0;
129        }
130    }
131    if (gps[procid]->neighbors[RIGHT] == -1) 
132    {
133        for (j = firstrow; j <= lastrow; j++) 
134        {
135            t2a[j][jm - 1] = 0.0;
136        }
137    }
138    for (i = firstrow; i <= lastrow; i++) 
139    {
140        t1a = (double *) t2a[i];
141        for (iindex = firstcol; iindex <= lastcol; iindex++) 
142        {
143            t1a[iindex] = 0.0;
144        }
145    }
146
147    t2a = (double **) gb[procid];
148    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) 
149    {
150        t2a[0][0] = 0.0;
151    }
152    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) 
153    {
154        t2a[im - 1][0] = 0.0;
155    }
156    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) 
157    {
158        t2a[0][jm - 1] = 0.0;
159    }
160    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) 
161    {
162        t2a[im - 1][jm - 1] = 0.0;
163    }
164    if (gps[procid]->neighbors[UP] == -1) 
165    {
166        t1a = (double *) t2a[0];
167        for (j = firstcol; j <= lastcol; j++)
168        {
169            t1a[j] = 0.0;
170        }
171    }
172    if (gps[procid]->neighbors[DOWN] == -1) 
173    {
174        t1a = (double *) t2a[im - 1];
175        for (j = firstcol; j <= lastcol; j++) 
176        {
177            t1a[j] = 0.0;
178        }
179    }
180    if (gps[procid]->neighbors[LEFT] == -1) 
181    {
182        for (j = firstrow; j <= lastrow; j++) 
183        {
184            t2a[j][0] = 0.0;
185        }
186    }
187    if (gps[procid]->neighbors[RIGHT] == -1) 
188    {
189        for (j = firstrow; j <= lastrow; j++) 
190        {
191            t2a[j][jm - 1] = 0.0;
192        }
193    }
194    for (i = firstrow; i <= lastrow; i++) 
195    {
196        t1a = (double *) t2a[i];
197        for (iindex = firstcol; iindex <= lastcol; iindex++) 
198        {
199            t1a[iindex] = 0.0;
200        }
201    }
202
203/* put the laplacian of psi{1,3} in work1{1,2}
204   note that psi(i,j,2) represents the psi3 array in
205   the original equations  */
206
207    for (psiindex = 0; psiindex <= 1; psiindex++) 
208    {
209        t2a = (double **) work1[procid][psiindex];
210        if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) 
211        {
212            t2a[0][0] = 0;
213        }
214        if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) 
215        {
216            t2a[im - 1][0] = 0;
217        }
218        if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) 
219        {
220            t2a[0][jm - 1] = 0;
221        }
222        if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) 
223        {
224            t2a[im - 1][jm - 1] = 0;
225        }
226
227        laplacalc(procid, psi, work1, psiindex, firstrow, lastrow, firstcol, lastcol);
228    }
229
230/* set values of work2 array to psi1 - psi3   */
231
232    t2a = (double **) work2[procid];
233    t2b = (double **) psi[procid][0];
234    t2c = (double **) psi[procid][1];
235    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
236        t2a[0][0] = t2b[0][0] - t2c[0][0];
237    }
238    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
239        t2a[im - 1][0] = t2b[im - 1][0] - t2c[im - 1][0];
240    }
241    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
242        t2a[0][jm - 1] = t2b[0][jm - 1] - t2c[0][jm - 1];
243    }
244    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
245        t2a[im - 1][jm - 1] = t2b[im - 1][jm - 1] - t2c[im - 1][jm - 1];
246    }
247    if (gps[procid]->neighbors[UP] == -1) {
248        t1a = (double *) t2a[0];
249        t1b = (double *) t2b[0];
250        t1c = (double *) t2c[0];
251        for (j = firstcol; j <= lastcol; j++) {
252            t1a[j] = t1b[j] - t1c[j];
253        }
254    }
255    if (gps[procid]->neighbors[DOWN] == -1) {
256        t1a = (double *) t2a[im - 1];
257        t1b = (double *) t2b[im - 1];
258        t1c = (double *) t2c[im - 1];
259        for (j = firstcol; j <= lastcol; j++) {
260            t1a[j] = t1b[j] - t1c[j];
261        }
262    }
263    if (gps[procid]->neighbors[LEFT] == -1) {
264        for (j = firstrow; j <= lastrow; j++) {
265            t2a[j][0] = t2b[j][0] - t2c[j][0];
266        }
267    }
268    if (gps[procid]->neighbors[RIGHT] == -1) {
269        for (j = firstrow; j <= lastrow; j++) {
270            t2a[j][jm - 1] = t2b[j][jm - 1] - t2c[j][jm - 1];
271        }
272    }
273    for (i = firstrow; i <= lastrow; i++) {
274        t1a = (double *) t2a[i];
275        t1b = (double *) t2b[i];
276        t1c = (double *) t2c[i];
277        for (iindex = firstcol; iindex <= lastcol; iindex++) {
278            t1a[iindex] = t1b[iindex] - t1c[iindex];
279        }
280    }
281
282/* set values of work3 array to h3/h * psi1 + h1/h * psi3  */
283
284    t2a = (double **) work3[procid];
285    hh3 = h3 / h;
286    hh1 = h1 / h;
287    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
288        t2a[0][0] = hh3 * t2a[0][0] + hh1 * t2c[0][0];
289    }
290    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
291        t2a[im - 1][0] = hh3 * t2a[im - 1][0] + hh1 * t2c[im - 1][0];
292    }
293    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
294        t2a[0][jm - 1] = hh3 * t2a[0][jm - 1] + hh1 * t2c[0][jm - 1];
295    }
296    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
297        t2a[im - 1][jm - 1] = hh3 * t2a[im - 1][jm - 1] + hh1 * t2c[im - 1][jm - 1];
298    }
299    if (gps[procid]->neighbors[UP] == -1) {
300        for (j = firstcol; j <= lastcol; j++) {
301            t2a[0][j] = hh3 * t2a[0][j] + hh1 * t2c[0][j];
302        }
303    }
304    if (gps[procid]->neighbors[DOWN] == -1) {
305        for (j = firstcol; j <= lastcol; j++) {
306            t2a[im - 1][j] = hh3 * t2a[im - 1][j] + hh1 * t2c[im - 1][j];
307        }
308    }
309    if (gps[procid]->neighbors[LEFT] == -1) {
310        for (j = firstrow; j <= lastrow; j++) {
311            t2a[j][0] = hh3 * t2a[j][0] + hh1 * t2c[j][0];
312        }
313    }
314    if (gps[procid]->neighbors[RIGHT] == -1) {
315        for (j = firstrow; j <= lastrow; j++) {
316            t2a[j][jm - 1] = hh3 * t2a[j][jm - 1] + hh1 * t2c[j][jm - 1];
317        }
318    }
319    for (i = firstrow; i <= lastrow; i++) {
320        t1a = (double *) t2a[i];
321        t1c = (double *) t2c[i];
322        for (iindex = firstcol; iindex <= lastcol; iindex++) {
323            t1a[iindex] = hh3 * t1a[iindex] + hh1 * t1c[iindex];
324        }
325    }
326
327/* set values of temparray{1,3} to psi{1,3}  */
328
329    for (psiindex = 0; psiindex <= 1; psiindex++) {
330        t2a = (double **) temparray[procid][psiindex];
331        t2b = (double **) psi[procid][psiindex];
332        if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
333            t2a[0][0] = t2b[0][0];
334        }
335        if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
336            t2a[im - 1][0] = t2b[im - 1][0];
337        }
338        if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
339            t2a[0][jm - 1] = t2b[0][jm - 1];
340        }
341        if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
342            t2a[im - 1][jm - 1] = t2b[im - 1][jm - 1];
343        }
344        if (gps[procid]->neighbors[UP] == -1) {
345            for (j = firstcol; j <= lastcol; j++) {
346                t2a[0][j] = t2b[0][j];
347            }
348        }
349        if (gps[procid]->neighbors[DOWN] == -1) {
350            for (j = firstcol; j <= lastcol; j++) {
351                t2a[im - 1][j] = t2b[im - 1][j];
352            }
353        }
354        if (gps[procid]->neighbors[LEFT] == -1) {
355            for (j = firstrow; j <= lastrow; j++) {
356                t2a[j][0] = t2b[j][0];
357            }
358        }
359        if (gps[procid]->neighbors[RIGHT] == -1) {
360            for (j = firstrow; j <= lastrow; j++) {
361                t2a[j][jm - 1] = t2b[j][jm - 1];
362            }
363        }
364
365        for (i = firstrow; i <= lastrow; i++) {
366            t1a = (double *) t2a[i];
367            t1b = (double *) t2b[i];
368            for (iindex = firstcol; iindex <= lastcol; iindex++) {
369                t1a[iindex] = t1b[iindex];
370            }
371        }
372    }
373
374    // END
375    gps[procid]->steps_time[0] = giet_proctime() - step_start;
376
377    // BARRIER
378    barrier_start = giet_proctime();
379    sqt_barrier_wait( &barrier );
380    gps[procid]->sync_time += (giet_proctime() - barrier_start);
381   
382/*     *******************************************************
383              s e c o n d   p h a s e
384       *******************************************************
385
386   set values of psi{1,3} to psim{1,3}   */
387
388    // BEGIN
389    step_start = giet_proctime();
390
391    for (psiindex = 0; psiindex <= 1; psiindex++) {
392        t2a = (double **) psi[procid][psiindex];
393        t2b = (double **) psim[procid][psiindex];
394        if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
395            t2a[0][0] = t2b[0][0];
396        }
397        if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
398            t2a[0][jm - 1] = t2b[0][jm - 1];
399        }
400        if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
401            t2a[im - 1][0] = t2b[im - 1][0];
402        }
403        if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
404            t2a[im - 1][jm - 1] = t2b[im - 1][jm - 1];
405        }
406        if (gps[procid]->neighbors[UP] == -1) {
407            for (j = firstcol; j <= lastcol; j++) {
408                t2a[0][j] = t2b[0][j];
409            }
410        }
411        if (gps[procid]->neighbors[DOWN] == -1) {
412            for (j = firstcol; j <= lastcol; j++) {
413                t2a[im - 1][j] = t2b[im - 1][j];
414            }
415        }
416        if (gps[procid]->neighbors[LEFT] == -1) {
417            for (j = firstrow; j <= lastrow; j++) {
418                t2a[j][0] = t2b[j][0];
419            }
420        }
421        if (gps[procid]->neighbors[RIGHT] == -1) {
422            for (j = firstrow; j <= lastrow; j++) {
423                t2a[j][jm - 1] = t2b[j][jm - 1];
424            }
425        }
426
427        for (i = firstrow; i <= lastrow; i++) {
428            t1a = (double *) t2a[i];
429            t1b = (double *) t2b[i];
430            for (iindex = firstcol; iindex <= lastcol; iindex++) {
431                t1a[iindex] = t1b[iindex];
432            }
433        }
434    }
435
436/* put the laplacian of the psim array
437   into the work7 array; first part of a three-laplacian
438   calculation to compute the friction terms  */
439
440    for (psiindex = 0; psiindex <= 1; psiindex++) {
441        t2a = (double **) work7[procid][psiindex];
442        if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
443            t2a[0][0] = 0;
444        }
445        if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
446            t2a[im - 1][0] = 0;
447        }
448        if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
449            t2a[0][jm - 1] = 0;
450        }
451        if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
452            t2a[im - 1][jm - 1] = 0;
453        }
454        laplacalc(procid, psim, work7, psiindex, firstrow, lastrow, firstcol, lastcol);
455    }
456
457/* to the values of the work1{1,2} arrays obtained from the
458   laplacians of psi{1,2} in the previous phase, add to the
459   elements of every column the corresponding value in the
460   one-dimenional f array  */
461
462    for (psiindex = 0; psiindex <= 1; psiindex++) {
463        t2a = (double **) work1[procid][psiindex];
464        if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
465            t2a[0][0] = t2a[0][0] + f[0];
466        }
467        if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
468            t2a[im - 1][0] = t2a[im - 1][0] + f[0];
469        }
470        if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
471            t2a[0][jm - 1] = t2a[0][jm - 1] + f[jmx[numlev - 1] - 1];
472        }
473        if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
474            t2a[im - 1][jm - 1] = t2a[im - 1][jm - 1] + f[jmx[numlev - 1] - 1];
475        }
476        if (gps[procid]->neighbors[UP] == -1) {
477            for (j = firstcol; j <= lastcol; j++) {
478                t2a[0][j] = t2a[0][j] + f[j + j_off];
479            }
480        }
481        if (gps[procid]->neighbors[DOWN] == -1) {
482            for (j = firstcol; j <= lastcol; j++) {
483                t2a[im - 1][j] = t2a[im - 1][j] + f[j + j_off];
484            }
485        }
486        if (gps[procid]->neighbors[LEFT] == -1) {
487            for (j = firstrow; j <= lastrow; j++) {
488                t2a[j][0] = t2a[j][0] + f[j + i_off];
489            }
490        }
491        if (gps[procid]->neighbors[RIGHT] == -1) {
492            for (j = firstrow; j <= lastrow; j++) {
493                t2a[j][jm - 1] = t2a[j][jm - 1] + f[j + i_off];
494            }
495        }
496        for (i = firstrow; i <= lastrow; i++) {
497            t1a = (double *) t2a[i];
498            for (iindex = firstcol; iindex <= lastcol; iindex++) {
499                t1a[iindex] = t1a[iindex] + f[iindex + j_off];
500            }
501        }
502    }
503
504    // END
505    gps[procid]->steps_time[1] = giet_proctime() - step_start;
506   
507    barrier_start = giet_proctime();
508    sqt_barrier_wait( &barrier );
509    gps[procid]->sync_time += giet_proctime() - barrier_start;
510
511/*      *******************************************************
512                 t h i r d   p h a s e
513        *******************************************************
514
515   put the jacobian of the work1{1,2} and psi{1,3} arrays
516   (the latter currently in temparray) in the work5{1,2} arrays  */
517
518    // BEGIN
519    step_start = giet_proctime();
520
521    for (psiindex = 0; psiindex <= 1; psiindex++) 
522    {
523        jacobcalc2( work1, 
524                    temparray, 
525                    work5, 
526                    psiindex, 
527                    procid, 
528                    firstrow, 
529                    lastrow, 
530                    firstcol, 
531                    lastcol );
532    }
533
534/* set values of psim{1,3} to temparray{1,3}  */
535
536    for (psiindex = 0; psiindex <= 1; psiindex++) {
537        t2a = (double **) psim[procid][psiindex];
538        t2b = (double **) temparray[procid][psiindex];
539        if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
540            t2a[0][0] = t2b[0][0];
541        }
542        if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
543            t2a[im - 1][0] = t2b[im - 1][0];
544        }
545        if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
546            t2a[0][jm - 1] = t2b[0][jm - 1];
547        }
548        if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
549            t2a[im - 1][jm - 1] = t2b[im - 1][jm - 1];
550        }
551        if (gps[procid]->neighbors[UP] == -1) {
552            t1a = (double *) t2a[0];
553            t1b = (double *) t2b[0];
554            for (j = firstcol; j <= lastcol; j++) {
555                t1a[j] = t1b[j];
556            }
557        }
558        if (gps[procid]->neighbors[DOWN] == -1) {
559            t1a = (double *) t2a[im - 1];
560            t1b = (double *) t2b[im - 1];
561            for (j = firstcol; j <= lastcol; j++) {
562                t1a[j] = t1b[j];
563            }
564        }
565        if (gps[procid]->neighbors[LEFT] == -1) {
566            for (j = firstrow; j <= lastrow; j++) {
567                t2a[j][0] = t2b[j][0];
568            }
569        }
570        if (gps[procid]->neighbors[RIGHT] == -1) {
571            for (j = firstrow; j <= lastrow; j++) {
572                t2a[j][jm - 1] = t2b[j][jm - 1];
573            }
574        }
575        for (i = firstrow; i <= lastrow; i++) {
576            t1a = (double *) t2a[i];
577            t1b = (double *) t2b[i];
578            for (iindex = firstcol; iindex <= lastcol; iindex++) {
579                t1a[iindex] = t1b[iindex];
580            }
581        }
582    }
583
584/* put the laplacian of the work7{1,2} arrays in the work4{1,2}
585   arrays; second step in the three-laplacian friction calculation  */
586
587    for (psiindex = 0; psiindex <= 1; psiindex++) 
588    {
589        laplacalc(procid, work7, work4, psiindex, firstrow, lastrow, firstcol, lastcol);
590    }
591
592    // END
593    gps[procid]->steps_time[2] = giet_proctime() - step_start;
594
595    // BARRIER
596    barrier_start = giet_proctime();
597    sqt_barrier_wait( &barrier );
598    gps[procid]->sync_time += (giet_proctime() - barrier_start);
599
600/*     *******************************************************
601                f o u r t h   p h a s e
602       *******************************************************
603
604   put the jacobian of the work2 and work3 arrays in the work6
605   array  */
606
607    // BEGIN
608    step_start = giet_proctime();
609
610    jacobcalc(work2, work3, work6, procid, firstrow, lastrow, firstcol, lastcol);
611
612/* put the laplacian of the work4{1,2} arrays in the work7{1,2}
613   arrays; third step in the three-laplacian friction calculation  */
614
615    for (psiindex = 0; psiindex <= 1; psiindex++) {
616        laplacalc(procid, work4, work7, psiindex, firstrow, lastrow, firstcol, lastcol);
617    }
618   
619    // END
620    gps[procid]->steps_time[3] = giet_proctime() - step_start;
621
622    // BARRIER
623    barrier_start = giet_proctime();
624    sqt_barrier_wait( &barrier );
625    gps[procid]->sync_time += (giet_proctime() - barrier_start);
626
627/*     *******************************************************
628                f i f t h   p h a s e
629       *******************************************************
630
631   use the values of the work5, work6 and work7 arrays
632   computed in the previous time-steps to compute the
633   ga and gb arrays   */
634
635    // BEGIN
636    step_start = giet_proctime();
637
638    hinv = 1.0 / h;
639    h1inv = 1.0 / h1;
640
641    t2a = (double **) ga[procid];
642    t2b = (double **) gb[procid];
643    t2c = (double **) work5[procid][0];
644    t2d = (double **) work5[procid][1];
645    t2e = (double **) work7[procid][0];
646    t2f = (double **) work7[procid][1];
647    t2g = (double **) work6[procid];
648    t2h = (double **) tauz[procid];
649    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
650        t2a[0][0] = t2c[0][0] - t2d[0][0] + eig2 * t2g[0][0] + h1inv * t2h[0][0] + lf * t2e[0][0] - lf * t2f[0][0];
651        t2b[0][0] = hh1 * t2c[0][0] + hh3 * t2d[0][0] + hinv * t2h[0][0] + lf * hh1 * t2e[0][0] + lf * hh3 * t2f[0][0];
652    }
653    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
654        t2a[im - 1][0] = t2c[im - 1][0] - t2d[im - 1][0] + eig2 * t2g[im - 1][0] + h1inv * t2h[im - 1][0] + lf * t2e[im - 1][0] - lf * t2f[im - 1][0];
655        t2b[im - 1][0] = hh1 * t2c[im - 1][0] + hh3 * t2d[im - 1][0] + hinv * t2h[im - 1][0] + lf * hh1 * t2e[im - 1][0] + lf * hh3 * t2f[im - 1][0];
656    }
657    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
658        t2a[0][jm - 1] = t2c[0][jm - 1] - t2d[0][jm - 1] + eig2 * t2g[0][jm - 1] + h1inv * t2h[0][jm - 1] + lf * t2e[0][jm - 1] - lf * t2f[0][jm - 1];
659        t2b[0][jm - 1] = hh1 * t2c[0][jm - 1] + hh3 * t2d[0][jm - 1] + hinv * t2h[0][jm - 1] + lf * hh1 * t2e[0][jm - 1] + lf * hh3 * t2f[0][jm - 1];
660    }
661    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
662        t2a[im - 1][jm - 1] = t2c[im - 1][jm - 1] - t2d[im - 1][jm - 1] + eig2 * t2g[im - 1][jm - 1] + h1inv * t2h[im - 1][jm - 1] + lf * t2e[im - 1][jm - 1] - lf * t2f[im - 1][jm - 1];
663        t2b[im - 1][jm - 1] = hh1 * t2c[im - 1][jm - 1] + hh3 * t2d[im - 1][jm - 1] + hinv * t2h[im - 1][jm - 1] + lf * hh1 * t2e[im - 1][jm - 1] + lf * hh3 * t2f[im - 1][jm - 1];
664    }
665    if (gps[procid]->neighbors[UP] == -1) 
666    {
667        t1a = (double *) t2a[0];
668        t1b = (double *) t2b[0];
669        t1c = (double *) t2c[0];
670        t1d = (double *) t2d[0];
671        t1e = (double *) t2e[0];
672        t1f = (double *) t2f[0];
673        t1g = (double *) t2g[0];
674        t1h = (double *) t2h[0];
675        for (j = firstcol; j <= lastcol; j++) 
676        {
677            t1a[j] = t1c[j] - t1d[j] + eig2 * t1g[j] + h1inv * t1h[j] + lf * t1e[j] - lf * t1f[j];
678            t1b[j] = hh1 * t1c[j] + hh3 * t1d[j] + hinv * t1h[j] + lf * hh1 * t1e[j] + lf * hh3 * t1f[j];
679        }
680    }
681    if (gps[procid]->neighbors[DOWN] == -1) {
682        t1a = (double *) t2a[im - 1];
683        t1b = (double *) t2b[im - 1];
684        t1c = (double *) t2c[im - 1];
685        t1d = (double *) t2d[im - 1];
686        t1e = (double *) t2e[im - 1];
687        t1f = (double *) t2f[im - 1];
688        t1g = (double *) t2g[im - 1];
689        t1h = (double *) t2h[im - 1];
690        for (j = firstcol; j <= lastcol; j++) {
691            t1a[j] = t1c[j] - t1d[j] + eig2 * t1g[j] + h1inv * t1h[j] + lf * t1e[j] - lf * t1f[j];
692            t1b[j] = hh1 * t1c[j] + hh3 * t1d[j] + hinv * t1h[j] + lf * hh1 * t1e[j] + lf * hh3 * t1f[j];
693        }
694    }
695    if (gps[procid]->neighbors[LEFT] == -1) {
696        for (j = firstrow; j <= lastrow; j++) {
697            t2a[j][0] = t2c[j][0] - t2d[j][0] + eig2 * t2g[j][0] + h1inv * t2h[j][0] + lf * t2e[j][0] - lf * t2f[j][0];
698            t2b[j][0] = hh1 * t2c[j][0] + hh3 * t2d[j][0] + hinv * t2h[j][0] + lf * hh1 * t2e[j][0] + lf * hh3 * t2f[j][0];
699        }
700    }
701    if (gps[procid]->neighbors[RIGHT] == -1) {
702        for (j = firstrow; j <= lastrow; j++) {
703            t2a[j][jm - 1] = t2c[j][jm - 1] - t2d[j][jm - 1] + eig2 * t2g[j][jm - 1] + h1inv * t2h[j][jm - 1] + lf * t2e[j][jm - 1] - lf * t2f[j][jm - 1];
704            t2b[j][jm - 1] = hh1 * t2c[j][jm - 1] + hh3 * t2d[j][jm - 1] + hinv * t2h[j][jm - 1] + lf * hh1 * t2e[j][jm - 1] + lf * hh3 * t2f[j][jm - 1];
705        }
706    }
707
708    for (i = firstrow; i <= lastrow; i++) {
709        t1a = (double *) t2a[i];
710        t1b = (double *) t2b[i];
711        t1c = (double *) t2c[i];
712        t1d = (double *) t2d[i];
713        t1e = (double *) t2e[i];
714        t1f = (double *) t2f[i];
715        t1g = (double *) t2g[i];
716        t1h = (double *) t2h[i];
717        for (iindex = firstcol; iindex <= lastcol; iindex++) {
718            t1a[iindex] = t1c[iindex] - t1d[iindex] + eig2 * t1g[iindex] + h1inv * t1h[iindex] + lf * t1e[iindex] - lf * t1f[iindex];
719            t1b[iindex] = hh1 * t1c[iindex] + hh3 * t1d[iindex] + hinv * t1h[iindex] + lf * hh1 * t1e[iindex] + lf * hh3 * t1f[iindex];
720        }
721    }
722
723    // END
724    gps[procid]->steps_time[4] = giet_proctime() - step_start;
725
726    // BARRIER
727    barrier_start = giet_proctime();
728    sqt_barrier_wait( &barrier );
729    gps[procid]->sync_time += (giet_proctime() - barrier_start);
730
731/*     *******************************************************
732               s i x t h   p h a s e
733       *******************************************************  */
734
735    // BEGIN
736    step_start = giet_proctime();
737
738    istart = 1;
739    iend = istart + gps[procid]->rel_num_y[numlev - 1] - 1;
740    jstart = 1;
741    jend = jstart + gps[procid]->rel_num_x[numlev - 1] - 1;
742    ist = istart;
743    ien = iend;
744    jst = jstart;
745    jen = jend;
746
747    if (gps[procid]->neighbors[UP] == -1) {
748        istart = 0;
749    }
750    if (gps[procid]->neighbors[LEFT] == -1) {
751        jstart = 0;
752    }
753    if (gps[procid]->neighbors[DOWN] == -1) {
754        iend = im - 1;
755    }
756    if (gps[procid]->neighbors[RIGHT] == -1) {
757        jend = jm - 1;
758    }
759    t2a = (double **) rhs_multi[procid][numlev - 1];
760    t2b = (double **) ga[procid];
761    t2c = (double **) oldga[procid];
762    t2d = (double **) q_multi[procid][numlev - 1];
763    for (i = istart; i <= iend; i++) {
764        t1a = (double *) t2a[i];
765        t1b = (double *) t2b[i];
766        for (j = jstart; j <= jend; j++) {
767            t1a[j] = t1b[j] * ressqr;
768        }
769    }
770
771    if (gps[procid]->neighbors[UP] == -1) {
772        t1d = (double *) t2d[0];
773        t1b = (double *) t2b[0];
774        for (j = jstart; j <= jend; j++) {
775            t1d[j] = t1b[j];
776        }
777    }
778    if (gps[procid]->neighbors[DOWN] == -1) {
779        t1d = (double *) t2d[im - 1];
780        t1b = (double *) t2b[im - 1];
781        for (j = jstart; j <= jend; j++) {
782            t1d[j] = t1b[j];
783        }
784    }
785    if (gps[procid]->neighbors[LEFT] == -1) {
786        for (i = istart; i <= iend; i++) {
787            t2d[i][0] = t2b[i][0];
788        }
789    }
790    if (gps[procid]->neighbors[RIGHT] == -1) {
791        for (i = istart; i <= iend; i++) {
792            t2d[i][jm - 1] = t2b[i][jm - 1];
793        }
794    }
795    //fac = 1.0 / (4.0 - ressqr*eig2);
796    for (i = ist; i <= ien; i++) {
797        t1d = (double *) t2d[i];
798        t1c = (double *) t2c[i];
799        for (j = jst; j <= jen; j++) {
800            t1d[j] = t1c[j];
801        }
802    }
803
804    multi_start = giet_proctime();
805    multig(procid);
806    gps[procid]->multi_time += (giet_proctime() - multi_start); 
807
808/* the shared sum variable psiai is initialized to 0 at
809   every time-step  */
810
811    if (procid == MASTER) 
812    {
813        global->psiai = 0.0;
814    }
815
816/*  copy the solution for use as initial guess in next time-step  */
817
818    for (i = istart; i <= iend; i++) 
819    {
820        t1b = (double *) t2b[i];
821        t1c = (double *) t2c[i];
822        t1d = (double *) t2d[i];
823        for (j = jstart; j <= jend; j++) 
824        {
825            t1b[j] = t1d[j];
826            t1c[j] = t1d[j];
827        }
828    }
829
830    // END
831    gps[procid]->steps_time[5] = giet_proctime() - step_start;
832
833    // BARRIER
834    barrier_start = giet_proctime();
835    sqt_barrier_wait( &barrier );
836    gps[procid]->sync_time += (giet_proctime() - barrier_start);
837
838/*     *******************************************************
839                s e v e n t h   p h a s e
840       *******************************************************
841
842   every process computes the running sum for its assigned portion
843   in a private variable psiaipriv   */
844
845    // BEGIN
846    step_start = giet_proctime();
847
848    psiaipriv = 0.0;
849    t2a = (double **) ga[procid];
850    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
851        psiaipriv = psiaipriv + 0.25 * (t2a[0][0]);
852    }
853    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
854        psiaipriv = psiaipriv + 0.25 * (t2a[0][jm - 1]);
855    }
856    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
857        psiaipriv = psiaipriv + 0.25 * (t2a[im - 1][0]);
858    }
859    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
860        psiaipriv = psiaipriv + 0.25 * (t2a[im - 1][jm - 1]);
861    }
862    if (gps[procid]->neighbors[UP] == -1) {
863        t1a = (double *) t2a[0];
864        for (j = firstcol; j <= lastcol; j++) {
865            psiaipriv = psiaipriv + 0.5 * t1a[j];
866        }
867    }
868    if (gps[procid]->neighbors[DOWN] == -1) {
869        t1a = (double *) t2a[im - 1];
870        for (j = firstcol; j <= lastcol; j++) {
871            psiaipriv = psiaipriv + 0.5 * t1a[j];
872        }
873    }
874    if (gps[procid]->neighbors[LEFT] == -1) {
875        for (j = firstrow; j <= lastrow; j++) {
876            psiaipriv = psiaipriv + 0.5 * t2a[j][0];
877        }
878    }
879    if (gps[procid]->neighbors[RIGHT] == -1) {
880        for (j = firstrow; j <= lastrow; j++) {
881            psiaipriv = psiaipriv + 0.5 * t2a[j][jm - 1];
882        }
883    }
884    for (i = firstrow; i <= lastrow; i++) {
885        t1a = (double *) t2a[i];
886        for (iindex = firstcol; iindex <= lastcol; iindex++) {
887            psiaipriv = psiaipriv + t1a[iindex];
888        }
889    }
890
891/* after computing its private sum, every process adds that to the
892   shared running sum psiai  */
893
894    sqt_lock_acquire( &psiai_lock );
895    global->psiai = global->psiai + psiaipriv;
896    sqt_lock_release( &psiai_lock );
897
898    // END
899    gps[procid]->steps_time[6] = giet_proctime() - step_start;
900
901    // BARRIER
902    barrier_start = giet_proctime();
903    sqt_barrier_wait( &barrier );
904    gps[procid]->sync_time += (giet_proctime() - barrier_start);
905   
906/*      *******************************************************
907                e i g h t h   p h a s e
908        *******************************************************
909
910   augment ga(i,j) with [-psiai/psibi]*psib(i,j) */
911
912    // BEGIN
913    step_start = giet_proctime();
914
915    f4 = (-global->psiai) /(global->psibi);
916
917    t2a = (double **) ga[procid];
918    t2b = (double **) psib[procid];
919    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
920        t2a[0][0] = t2a[0][0] + f4 * t2b[0][0];
921    }
922    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
923        t2a[im - 1][0] = t2a[im - 1][0] + f4 * t2b[im - 1][0];
924    }
925    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
926        t2a[0][jm - 1] = t2a[0][jm - 1] + f4 * t2b[0][jm - 1];
927    }
928    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
929        t2a[im - 1][jm - 1] = t2a[im - 1][jm - 1] + f4 * t2b[im - 1][jm - 1];
930    }
931    if (gps[procid]->neighbors[UP] == -1) {
932        t1a = (double *) t2a[0];
933        t1b = (double *) t2b[0];
934        for (j = firstcol; j <= lastcol; j++) {
935            t1a[j] = t1a[j] + f4 * t1b[j];
936        }
937    }
938    if (gps[procid]->neighbors[DOWN] == -1) {
939        t1a = (double *) t2a[im - 1];
940        t1b = (double *) t2b[im - 1];
941        for (j = firstcol; j <= lastcol; j++) {
942            t1a[j] = t1a[j] + f4 * t1b[j];
943        }
944    }
945    if (gps[procid]->neighbors[LEFT] == -1) {
946        for (j = firstrow; j <= lastrow; j++) {
947            t2a[j][0] = t2a[j][0] + f4 * t2b[j][0];
948        }
949    }
950    if (gps[procid]->neighbors[RIGHT] == -1) {
951        for (j = firstrow; j <= lastrow; j++) {
952            t2a[j][jm - 1] = t2a[j][jm - 1] + f4 * t2b[j][jm - 1];
953        }
954    }
955    for (i = firstrow; i <= lastrow; i++) {
956        t1a = (double *) t2a[i];
957        t1b = (double *) t2b[i];
958        for (iindex = firstcol; iindex <= lastcol; iindex++) {
959            t1a[iindex] = t1a[iindex] + f4 * t1b[iindex];
960        }
961    }
962
963    t2a = (double **) rhs_multi[procid][numlev - 1];
964    t2b = (double **) gb[procid];
965    t2c = (double **) oldgb[procid];
966    t2d = (double **) q_multi[procid][numlev - 1];
967    for (i = istart; i <= iend; i++) {
968        t1a = (double *) t2a[i];
969        t1b = (double *) t2b[i];
970        for (j = jstart; j <= jend; j++) {
971            t1a[j] = t1b[j] * ressqr;
972        }
973    }
974    if (gps[procid]->neighbors[UP] == -1) {
975        t1d = (double *) t2d[0];
976        t1b = (double *) t2b[0];
977        for (j = jstart; j <= jend; j++) {
978            t1d[j] = t1b[j];
979        }
980    }
981    if (gps[procid]->neighbors[DOWN] == -1) {
982        t1d = (double *) t2d[im - 1];
983        t1b = (double *) t2b[im - 1];
984        for (j = jstart; j <= jend; j++) {
985            t1d[j] = t1b[j];
986        }
987    }
988    if (gps[procid]->neighbors[LEFT] == -1) {
989        for (i = istart; i <= iend; i++) {
990            t2d[i][0] = t2b[i][0];
991        }
992    }
993    if (gps[procid]->neighbors[RIGHT] == -1) {
994        for (i = istart; i <= iend; i++) {
995            t2d[i][jm - 1] = t2b[i][jm - 1];
996        }
997    }
998    //fac = 1.0 / (4.0 - ressqr*eig2);
999    for (i = ist; i <= ien; i++) {
1000        t1d = (double *) t2d[i];
1001        t1c = (double *) t2c[i];
1002        for (j = jst; j <= jen; j++) {
1003            t1d[j] = t1c[j];
1004        }
1005    }
1006
1007    multi_start = giet_proctime();
1008    multig(procid);
1009    gps[procid]->multi_time += (giet_proctime() - multi_start);
1010
1011    for (i = istart; i <= iend; i++) 
1012    {
1013        t1b = (double *) t2b[i];
1014        t1c = (double *) t2c[i];
1015        t1d = (double *) t2d[i];
1016        for (j = jstart; j <= jend; j++) 
1017        {
1018            t1b[j] = t1d[j];
1019            t1c[j] = t1d[j];
1020        }
1021    }
1022
1023    // END
1024    gps[procid]->steps_time[7] = giet_proctime() - step_start;
1025
1026    // BARRIER
1027    barrier_start = giet_proctime();
1028    sqt_barrier_wait( &barrier );
1029    gps[procid]->sync_time += (giet_proctime() - barrier_start);
1030
1031/*      *******************************************************
1032                n i n t h   p h a s e
1033        *******************************************************
1034
1035   put appropriate linear combinations of ga and gb in work2 and work3;
1036   note that here (as in most cases) the constant multipliers are made
1037   private variables; the specific order in which things are done is
1038   chosen in order to hopefully reuse things brought into the cache
1039
1040   note that here again we choose to have all processes share the work
1041   on both matrices despite the fact that the work done per element
1042   is the same, because the operand matrices are the same in both cases */
1043
1044    // BEGIN
1045    step_start = giet_proctime();
1046
1047    t2a = (double **) ga[procid];
1048    t2b = (double **) gb[procid];
1049    t2c = (double **) work2[procid];
1050    t2d = (double **) work3[procid];
1051    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
1052        t2c[0][0] = t2b[0][0] - hh1 * t2a[0][0];
1053        t2d[0][0] = t2b[0][0] + hh3 * t2a[0][0];
1054    }
1055    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) {
1056        t2c[im - 1][0] = t2b[im - 1][0] - hh1 * t2a[im - 1][0];
1057        t2d[im - 1][0] = t2b[im - 1][0] + hh3 * t2a[im - 1][0];
1058    }
1059    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
1060        t2c[0][jm - 1] = t2b[0][jm - 1] - hh1 * t2a[0][jm - 1];
1061        t2d[0][jm - 1] = t2b[0][jm - 1] + hh3 * t2a[0][jm - 1];
1062    }
1063    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) {
1064        t2c[im - 1][jm - 1] = t2b[im - 1][jm - 1] - hh1 * t2a[im - 1][jm - 1];
1065        t2d[im - 1][jm - 1] = t2b[im - 1][jm - 1] + hh3 * t2a[im - 1][jm - 1];
1066    }
1067    if (gps[procid]->neighbors[UP] == -1) {
1068        t1a = (double *) t2a[0];
1069        t1b = (double *) t2b[0];
1070        t1c = (double *) t2c[0];
1071        t1d = (double *) t2d[0];
1072        for (j = firstcol; j <= lastcol; j++) {
1073            t1d[j] = t1b[j] + hh3 * t1a[j];
1074            t1c[j] = t1b[j] - hh1 * t1a[j];
1075        }
1076    }
1077    if (gps[procid]->neighbors[DOWN] == -1) {
1078        t1a = (double *) t2a[im - 1];
1079        t1b = (double *) t2b[im - 1];
1080        t1c = (double *) t2c[im - 1];
1081        t1d = (double *) t2d[im - 1];
1082        for (j = firstcol; j <= lastcol; j++) {
1083            t1d[j] = t1b[j] + hh3 * t1a[j];
1084            t1c[j] = t1b[j] - hh1 * t1a[j];
1085        }
1086    }
1087    if (gps[procid]->neighbors[LEFT] == -1) {
1088        for (j = firstrow; j <= lastrow; j++) {
1089            t2d[j][0] = t2b[j][0] + hh3 * t2a[j][0];
1090            t2c[j][0] = t2b[j][0] - hh1 * t2a[j][0];
1091        }
1092    }
1093    if (gps[procid]->neighbors[RIGHT] == -1) {
1094        for (j = firstrow; j <= lastrow; j++) {
1095            t2d[j][jm - 1] = t2b[j][jm - 1] + hh3 * t2a[j][jm - 1];
1096            t2c[j][jm - 1] = t2b[j][jm - 1] - hh1 * t2a[j][jm - 1];
1097        }
1098    }
1099
1100    for (i = firstrow; i <= lastrow; i++) {
1101        t1a = (double *) t2a[i];
1102        t1b = (double *) t2b[i];
1103        t1c = (double *) t2c[i];
1104        t1d = (double *) t2d[i];
1105        for (iindex = firstcol; iindex <= lastcol; iindex++) {
1106            t1d[iindex] = t1b[iindex] + hh3 * t1a[iindex];
1107            t1c[iindex] = t1b[iindex] - hh1 * t1a[iindex];
1108        }
1109    }
1110
1111    // END
1112    gps[procid]->steps_time[8] = giet_proctime() - step_start;
1113
1114    // BARRIER
1115    barrier_start = giet_proctime();
1116    sqt_barrier_wait( &barrier );
1117    gps[procid]->sync_time += (giet_proctime() - barrier_start);
1118
1119/*      *******************************************************
1120                t e n t h    p h a s e
1121        *******************************************************/
1122
1123    // BEGIN
1124    step_start = giet_proctime();
1125
1126    timst = 2 * dtau;
1127
1128/* update the psi{1,3} matrices by adding 2*dtau*work3 to each */
1129
1130    t2a = (double **) psi[procid][0];
1131    t2b = (double **) work3[procid];
1132    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) 
1133    {
1134        t2a[0][0] = t2a[0][0] + timst * t2b[0][0];
1135    }
1136    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) 
1137    {
1138        t2a[im - 1][0] = t2a[im - 1][0] + timst * t2b[im - 1][0];
1139    }
1140    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) 
1141    {
1142        t2a[0][jm - 1] = t2a[0][jm - 1] + timst * t2b[0][jm - 1];
1143    }
1144    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) 
1145    {
1146        t2a[im - 1][jm - 1] = t2a[im - 1][jm - 1] + timst * t2b[im - 1][jm - 1];
1147    }
1148    if (gps[procid]->neighbors[UP] == -1) 
1149    {
1150        t1a = (double *) t2a[0];
1151        t1b = (double *) t2b[0];
1152        for (j = firstcol; j <= lastcol; j++) 
1153        {
1154            t1a[j] = t1a[j] + timst * t1b[j];
1155        }
1156    }
1157    if (gps[procid]->neighbors[DOWN] == -1) 
1158    {
1159        t1a = (double *) t2a[im - 1];
1160        t1b = (double *) t2b[im - 1];
1161        for (j = firstcol; j <= lastcol; j++) 
1162        {
1163            t1a[j] = t1a[j] + timst * t1b[j];
1164        }
1165    }
1166    if (gps[procid]->neighbors[LEFT] == -1) 
1167    {
1168        for (j = firstrow; j <= lastrow; j++) 
1169        {
1170            t2a[j][0] = t2a[j][0] + timst * t2b[j][0];
1171        }
1172    }
1173    if (gps[procid]->neighbors[RIGHT] == -1) 
1174    {
1175        for (j = firstrow; j <= lastrow; j++) 
1176        {
1177            t2a[j][jm - 1] = t2a[j][jm - 1] + timst * t2b[j][jm - 1];
1178        }
1179    }
1180    for (i = firstrow; i <= lastrow; i++) 
1181    {
1182        t1a = (double *) t2a[i];
1183        t1b = (double *) t2b[i];
1184        for (iindex = firstcol; iindex <= lastcol; iindex++) 
1185        {
1186            t1a[iindex] = t1a[iindex] + timst * t1b[iindex];
1187        }
1188    }
1189
1190    t2a = (double **) psi[procid][1];
1191    t2b = (double **) work2[procid];
1192    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[LEFT] == -1)) 
1193    {
1194        t2a[0][0] = t2a[0][0] + timst * t2b[0][0];
1195    }
1196    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[LEFT] == -1)) 
1197    {
1198        t2a[im - 1][0] = t2a[im - 1][0] + timst * t2b[im - 1][0];
1199    }
1200    if ((gps[procid]->neighbors[UP] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) 
1201    {
1202        t2a[0][jm - 1] = t2a[0][jm - 1] + timst * t2b[0][jm - 1];
1203    }
1204    if ((gps[procid]->neighbors[DOWN] == -1) && (gps[procid]->neighbors[RIGHT] == -1)) 
1205    {
1206        t2a[im - 1][jm - 1] = t2a[im - 1][jm - 1] + timst * t2b[im - 1][jm - 1];
1207    }
1208    if (gps[procid]->neighbors[UP] == -1) 
1209    {
1210        t1a = (double *) t2a[0];
1211        t1b = (double *) t2b[0];
1212        for (j = firstcol; j <= lastcol; j++) 
1213        {
1214            t1a[j] = t1a[j] + timst * t1b[j];
1215        }
1216    }
1217    if (gps[procid]->neighbors[DOWN] == -1) 
1218    {
1219        t1a = (double *) t2a[im - 1];
1220        t1b = (double *) t2b[im - 1];
1221        for (j = firstcol; j <= lastcol; j++) 
1222        {
1223            t1a[j] = t1a[j] + timst * t1b[j];
1224        }
1225    }
1226    if (gps[procid]->neighbors[LEFT] == -1) 
1227    {
1228        for (j = firstrow; j <= lastrow; j++) 
1229        {
1230            t2a[j][0] = t2a[j][0] + timst * t2b[j][0];
1231        }
1232    }
1233    if (gps[procid]->neighbors[RIGHT] == -1) 
1234    {
1235        for (j = firstrow; j <= lastrow; j++) 
1236        {
1237            t2a[j][jm - 1] = t2a[j][jm - 1] + timst * t2b[j][jm - 1];
1238        }
1239    }
1240
1241    for (i = firstrow; i <= lastrow; i++) 
1242    {
1243        t1a = (double *) t2a[i];
1244        t1b = (double *) t2b[i];
1245        for (iindex = firstcol; iindex <= lastcol; iindex++) 
1246        {
1247            t1a[iindex] = t1a[iindex] + timst * t1b[iindex];
1248        }
1249    }
1250
1251    // END
1252    gps[procid]->steps_time[9] = giet_proctime() - step_start;
1253
1254    // BARRIER
1255    barrier_start = giet_proctime();
1256    sqt_barrier_wait( &barrier );
1257    gps[procid]->sync_time += (giet_proctime() - barrier_start);
1258
1259} // end slave2()
Note: See TracBrowser for help on using the repository browser.