Changeset 589 for soft/giet_vm/applications/ocean/jacobcalc.C
- Timestamp:
- Jul 8, 2015, 3:57:15 PM (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
soft/giet_vm/applications/ocean/jacobcalc.C
r581 r589 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 /*************************************************************************/ 1 #line 115 "/Users/alain/soc/giet_vm/applications/ocean/null_macros/c.m4.null.GIET" 16 2 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 EXTERN_ENV21 22 #include <stdio.h>23 #include <math.h>24 25 #include "decs.h"26 27 void jacobcalc(double ***x, double ***y, double ***z, long pid, long firstrow, long lastrow, long firstcol, long lastcol)28 {29 double f1;30 double f2;31 double f3;32 double f4;33 double f5;34 double f6;35 double f7;36 double f8;37 long iindex;38 long indexp1;39 long indexm1;40 long im1;41 long ip1;42 long i;43 long j;44 long jj;45 double **t2a;46 double **t2b;47 double **t2c;48 double *t1a;49 double *t1b;50 double *t1c;51 double *t1d;52 double *t1e;53 double *t1f;54 double *t1g;55 56 t2a = (double **) z[pid];57 if ((gp[pid].neighbors[UP] == -1) && (gp[pid].neighbors[LEFT] == -1)) {58 t2a[0][0] = 0.0;59 }60 if ((gp[pid].neighbors[DOWN] == -1) && (gp[pid].neighbors[LEFT] == -1)) {61 t2a[im - 1][0] = 0.0;62 }63 if ((gp[pid].neighbors[UP] == -1) && (gp[pid].neighbors[RIGHT] == -1)) {64 t2a[0][jm - 1] = 0.0;65 }66 if ((gp[pid].neighbors[DOWN] == -1) && (gp[pid].neighbors[RIGHT] == -1)) {67 t2a[im - 1][jm - 1] = 0.0;68 }69 70 t2a = (double **) x[pid];71 jj = gp[pid].neighbors[UPLEFT];72 if (jj != -1) {73 t2a[0][0] = x[jj][im - 2][jm - 2];74 }75 jj = gp[pid].neighbors[UPRIGHT];76 if (jj != -1) {77 t2a[0][jm - 1] = x[jj][im - 2][1];78 }79 jj = gp[pid].neighbors[DOWNLEFT];80 if (jj != -1) {81 t2a[im - 1][0] = x[jj][1][jm - 2];82 }83 jj = gp[pid].neighbors[DOWNRIGHT];84 if (jj != -1) {85 t2a[im - 1][jm - 1] = x[jj][1][1];86 }87 88 t2a = (double **) y[pid];89 jj = gp[pid].neighbors[UPLEFT];90 if (jj != -1) {91 t2a[0][0] = y[jj][im - 2][jm - 2];92 }93 jj = gp[pid].neighbors[UPRIGHT];94 if (jj != -1) {95 t2a[0][jm - 1] = y[jj][im - 2][1];96 }97 jj = gp[pid].neighbors[DOWNLEFT];98 if (jj != -1) {99 t2a[im - 1][0] = y[jj][1][jm - 2];100 }101 jj = gp[pid].neighbors[DOWNRIGHT];102 if (jj != -1) {103 t2a[im - 1][jm - 1] = y[jj][1][1];104 }105 106 t2a = (double **) x[pid];107 if (gp[pid].neighbors[UP] == -1) {108 jj = gp[pid].neighbors[LEFT];109 if (jj != -1) {110 t2a[0][0] = x[jj][0][jm - 2];111 } else {112 jj = gp[pid].neighbors[DOWN];113 if (jj != -1) {114 t2a[im - 1][0] = x[jj][1][0];115 }116 }117 jj = gp[pid].neighbors[RIGHT];118 if (jj != -1) {119 t2a[0][jm - 1] = x[jj][0][1];120 } else {121 jj = gp[pid].neighbors[DOWN];122 if (jj != -1) {123 t2a[im - 1][jm - 1] = x[jj][1][jm - 1];124 }125 }126 } else if (gp[pid].neighbors[DOWN] == -1) {127 jj = gp[pid].neighbors[LEFT];128 if (jj != -1) {129 t2a[im - 1][0] = x[jj][im - 1][jm - 2];130 } else {131 jj = gp[pid].neighbors[UP];132 if (jj != -1) {133 t2a[0][0] = x[jj][im - 2][0];134 }135 }136 jj = gp[pid].neighbors[RIGHT];137 if (jj != -1) {138 t2a[im - 1][jm - 1] = x[jj][im - 1][1];139 } else {140 jj = gp[pid].neighbors[UP];141 if (jj != -1) {142 t2a[0][jm - 1] = x[jj][im - 2][jm - 1];143 }144 }145 } else if (gp[pid].neighbors[LEFT] == -1) {146 jj = gp[pid].neighbors[UP];147 if (jj != -1) {148 t2a[0][0] = x[jj][im - 2][0];149 }150 jj = gp[pid].neighbors[DOWN];151 if (jj != -1) {152 t2a[im - 1][0] = x[jj][1][0];153 }154 } else if (gp[pid].neighbors[RIGHT] == -1) {155 jj = gp[pid].neighbors[UP];156 if (jj != -1) {157 t2a[0][jm - 1] = x[jj][im - 2][jm - 1];158 }159 jj = gp[pid].neighbors[DOWN];160 if (jj != -1) {161 t2a[im - 1][jm - 1] = x[jj][1][jm - 1];162 }163 }164 165 t2a = (double **) y[pid];166 if (gp[pid].neighbors[UP] == -1) {167 jj = gp[pid].neighbors[LEFT];168 if (jj != -1) {169 t2a[0][0] = y[jj][0][jm - 2];170 } else {171 jj = gp[pid].neighbors[DOWN];172 if (jj != -1) {173 t2a[im - 1][0] = y[jj][1][0];174 }175 }176 jj = gp[pid].neighbors[RIGHT];177 if (jj != -1) {178 t2a[0][jm - 1] = y[jj][0][1];179 } else {180 jj = gp[pid].neighbors[DOWN];181 if (jj != -1) {182 t2a[im - 1][jm - 1] = y[jj][1][jm - 1];183 }184 }185 } else if (gp[pid].neighbors[DOWN] == -1) {186 jj = gp[pid].neighbors[LEFT];187 if (jj != -1) {188 t2a[im - 1][0] = y[jj][im - 1][jm - 2];189 } else {190 jj = gp[pid].neighbors[UP];191 if (jj != -1) {192 t2a[0][0] = y[jj][im - 2][0];193 }194 }195 jj = gp[pid].neighbors[RIGHT];196 if (jj != -1) {197 t2a[im - 1][jm - 1] = y[jj][im - 1][1];198 } else {199 jj = gp[pid].neighbors[UP];200 if (jj != -1) {201 t2a[0][jm - 1] = y[jj][im - 2][jm - 1];202 }203 }204 } else if (gp[pid].neighbors[LEFT] == -1) {205 jj = gp[pid].neighbors[UP];206 if (jj != -1) {207 t2a[0][0] = y[jj][im - 2][0];208 }209 jj = gp[pid].neighbors[DOWN];210 if (jj != -1) {211 t2a[im - 1][0] = y[jj][1][0];212 }213 } else if (gp[pid].neighbors[RIGHT] == -1) {214 jj = gp[pid].neighbors[UP];215 if (jj != -1) {216 t2a[0][jm - 1] = y[jj][im - 2][jm - 1];217 }218 jj = gp[pid].neighbors[DOWN];219 if (jj != -1) {220 t2a[im - 1][jm - 1] = y[jj][1][jm - 1];221 }222 }223 224 j = gp[pid].neighbors[UP];225 if (j != -1) {226 t1a = (double *) t2a[0];227 t1b = (double *) y[j][im - 2];228 for (i = 1; i <= lastcol; i++) {229 t1a[i] = t1b[i];230 }231 }232 j = gp[pid].neighbors[DOWN];233 if (j != -1) {234 t1a = (double *) t2a[im - 1];235 t1b = (double *) y[j][1];236 for (i = 1; i <= lastcol; i++) {237 t1a[i] = t1b[i];238 }239 }240 j = gp[pid].neighbors[LEFT];241 if (j != -1) {242 t2b = (double **) y[j];243 for (i = 1; i <= lastrow; i++) {244 t2a[i][0] = t2b[i][jm - 2];245 }246 }247 j = gp[pid].neighbors[RIGHT];248 if (j != -1) {249 t2b = (double **) y[j];250 for (i = 1; i <= lastrow; i++) {251 t2a[i][jm - 1] = t2b[i][1];252 }253 }254 255 t2a = (double **) x[pid];256 j = gp[pid].neighbors[UP];257 if (j != -1) {258 t1a = (double *) t2a[0];259 t1b = (double *) x[j][im - 2];260 for (i = 1; i <= lastcol; i++) {261 t1a[i] = t1b[i];262 }263 }264 j = gp[pid].neighbors[DOWN];265 if (j != -1) {266 t1a = (double *) t2a[im - 1];267 t1b = (double *) x[j][1];268 for (i = 1; i <= lastcol; i++) {269 t1a[i] = t1b[i];270 }271 }272 j = gp[pid].neighbors[LEFT];273 if (j != -1) {274 t2b = (double **) x[j];275 for (i = 1; i <= lastrow; i++) {276 t2a[i][0] = t2b[i][jm - 2];277 }278 }279 j = gp[pid].neighbors[RIGHT];280 if (j != -1) {281 t2b = (double **) x[j];282 for (i = 1; i <= lastrow; i++) {283 t2a[i][jm - 1] = t2b[i][1];284 }285 }286 287 t2a = (double **) x[pid];288 t2b = (double **) y[pid];289 t2c = (double **) z[pid];290 for (i = firstrow; i <= lastrow; i++) {291 ip1 = i + 1;292 im1 = i - 1;293 t1a = (double *) t2a[i];294 t1b = (double *) t2b[i];295 t1c = (double *) t2c[i];296 t1d = (double *) t2b[ip1];297 t1e = (double *) t2b[im1];298 t1f = (double *) t2a[ip1];299 t1g = (double *) t2a[im1];300 for (iindex = firstcol; iindex <= lastcol; iindex++) {301 indexp1 = iindex + 1;302 indexm1 = iindex - 1;303 f1 = (t1b[indexm1] + t1d[indexm1] - t1b[indexp1] - t1d[indexp1]) * (t1f[iindex] - t1a[iindex]);304 f2 = (t1e[indexm1] + t1b[indexm1] - t1e[indexp1] - t1b[indexp1]) * (t1a[iindex] - t1g[iindex]);305 f3 = (t1d[iindex] + t1d[indexp1] - t1e[iindex] - t1e[indexp1]) * (t1a[indexp1] - t1a[iindex]);306 f4 = (t1d[indexm1] + t1d[iindex] - t1e[indexm1] - t1e[iindex]) * (t1a[iindex] - t1a[indexm1]);307 f5 = (t1d[iindex] - t1b[indexp1]) * (t1f[indexp1] - t1a[iindex]);308 f6 = (t1b[indexm1] - t1e[iindex]) * (t1a[iindex] - t1g[indexm1]);309 f7 = (t1b[indexp1] - t1e[iindex]) * (t1g[indexp1] - t1a[iindex]);310 f8 = (t1d[iindex] - t1b[indexm1]) * (t1a[iindex] - t1f[indexm1]);311 312 t1c[iindex] = factjacob * (f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8);313 }314 }315 316 if (gp[pid].neighbors[UP] == -1) {317 t1c = (double *) t2c[0];318 for (j = firstcol; j <= lastcol; j++) {319 t1c[j] = 0.0;320 }321 }322 if (gp[pid].neighbors[DOWN] == -1) {323 t1c = (double *) t2c[im - 1];324 for (j = firstcol; j <= lastcol; j++) {325 t1c[j] = 0.0;326 }327 }328 if (gp[pid].neighbors[LEFT] == -1) {329 for (j = firstrow; j <= lastrow; j++) {330 t2c[j][0] = 0.0;331 }332 }333 if (gp[pid].neighbors[RIGHT] == -1) {334 for (j = firstrow; j <= lastrow; j++) {335 t2c[j][jm - 1] = 0.0;336 }337 }338 339 }
Note: See TracChangeset
for help on using the changeset viewer.