1//--------------------------------------------------------------------- 2//--------------------------------------------------------------------- 3#include "header.h" 4#include "mpinpb.h" 5#define G_MAIN 6#include "work_lhs.h" 7#undef G_MAIN 8 9extern void x_sendrecv_solve(int c, int cprev); 10extern void x_sendrecv_back(int c, int cprev); 11extern void x_backsubstitute(int first, int last, int c); 12extern void x_solve_cell(int first, int last, int c); 13 14void x_solve() { 15// printf("entered x_solve()...\n");fflush(0); 16 17//--------------------------------------------------------------------- 18//--------------------------------------------------------------------- 19 20//--------------------------------------------------------------------- 21// Performs line solves in X direction by first factoring 22// the block-tridiagonal matrix into an upper triangular matrix, 23// and then performing back substitution to solve for the unknown 24// vectors of each line. 25// 26// Make sure we treat elements zero to cell_size in the direction 27// of the sweep. 28//--------------------------------------------------------------------- 29 30 int c, cprev, stage, first, last, error; 31 32//--------------------------------------------------------------------- 33// in our terminology stage is the number of the cell in the x-direction 34// i.e. stage = 1 means the start of the line stage=ncells means end 35//--------------------------------------------------------------------- 36 for (stage = 1; stage <= ncells; stage++) { 37 c = slice(1,stage); 38// printf("passed slice.1(stage %0d)...\n", stage);fflush(0); 39//--------------------------------------------------------------------- 40// set first/last-cell flags 41//--------------------------------------------------------------------- 42 first = (stage == 1); 43 last = (stage == ncells); 44 45 if (stage >1) { 46 cprev = slice(1,stage-1); 47// printf("passed slice.2(stage %0d)...\n", stage);fflush(0); 48 x_sendrecv_solve(c, cprev); 49// printf("passed x_sendrecv_solve(stage %0d)...\n", stage);fflush(0); 50 } 51 x_solve_cell(first,last,c); 52// printf("passed x_solve_cell(stage %0d)...\n", stage);fflush(0); 53 } 54// printf("passed \"set first/last-cell flags\"...\n");fflush(0); 55 56//--------------------------------------------------------------------- 57// now perform backsubstitution in reverse direction 58//--------------------------------------------------------------------- 59 for (stage = ncells; stage >= 1; stage--) { 60 c = slice(1,stage); 61 first = (stage == 1); 62 last = (stage == ncells); 63 64 if (stage <ncells) { 65 cprev = slice(1,stage+1); 66 x_sendrecv_back(c, cprev); 67 } 68 69 x_backsubstitute(first,last,c); 70 } 71 72// printf("passed \"perform backsubstitution in reverse direction\"...\n");fflush(0); 73 return; 74} 75 76 77void x_sendrecv_solve(int c, int cprev) { 78 79//--------------------------------------------------------------------- 80// pack up and send C'(iend) and rhs'(iend) for 81// all j and k of previous cell 82//--------------------------------------------------------------------- 83 84 int j,k,m,n,isize,ptr, istart; 85 int phase; 86 int error, buffer_size; 87// printf("Entered x_sendrecv_solve()...\n");fflush(0); 88 89 isize = cell_size(1,cprev)-1; 90 buffer_size=MAX_CELL_DIM*MAX_CELL_DIM* 91 (BLOCK_SIZE*BLOCK_SIZE + BLOCK_SIZE); 92 93//--------------------------------------------------------------------- 94// pack up buffer 95//--------------------------------------------------------------------- 96 ptr = 0; 97 for (k = 0; k <= KMAX-1; k++) { 98 for (j = 0; j <= JMAX-1; j++) { 99 for (m = 1; m <= BLOCK_SIZE; m++) { 100 for (n = 1; n <= BLOCK_SIZE; n++) { 101 in_buffer(ptr+n) = lhsc(m,n,isize,j,k,cprev); 102 } 103 ptr = ptr+BLOCK_SIZE; 104 } 105 for (n = 1; n <= BLOCK_SIZE; n++) { 106 in_buffer(ptr+n) = rhs(n,isize,j,k,cprev); 107 } 108 ptr = ptr+BLOCK_SIZE; 109 } 110 } 111// printf("Passed \"pack up buffer\"...\n");fflush(0); 112 113//--------------------------------------------------------------------- 114// send and receive buffer 115//--------------------------------------------------------------------- 116 117 for (phase = 0; phase < 3; phase++) { 118 119 if (send_color[EASTDIR]==phase) 120// printf("before RCCE_send (in_buffer=%08x, size=%08d)\n", in_buffer, buffer_size*sizeof(double));fflush(0); 121 RCCE_send((char*)in_buffer, buffer_size*sizeof(double), successor(1)); 122// printf("after RCCE_send");fflush(0); 123 if (recv_color[EASTDIR]==phase) 124// printf("before RCCE_recv");fflush(0); 125 RCCE_recv((char*)out_buffer, buffer_size*sizeof(double), predecessor(1)); 126// printf("after RCCE_recv");fflush(0); 127 } 128// printf("Passed \"send and receive buffer\"...\n");fflush(0); 129 130//--------------------------------------------------------------------- 131// unpack buffer 132//--------------------------------------------------------------------- 133 istart = 0; 134 ptr = 0; 135 for (k = 0; k <= KMAX-1; k++) { 136 for (j = 0; j <= JMAX-1; j++) { 137 for (m = 1; m <= BLOCK_SIZE; m++) { 138 for (n = 1; n <= BLOCK_SIZE; n++) { 139 lhsc(m,n,istart-1,j,k,c) = out_buffer(ptr+n); 140 } 141 ptr = ptr+BLOCK_SIZE; 142 } 143 for (n = 1; n <= BLOCK_SIZE; n++) { 144 rhs(n,istart-1,j,k,c) = out_buffer(ptr+n); 145 } 146 ptr = ptr+BLOCK_SIZE; 147 } 148 } 149// printf("Passed \"unpack buffer\"...\n");fflush(0); 150 151 return; 152} 153 154//--------------------------------------------------------------------- 155//--------------------------------------------------------------------- 156 157void x_sendrecv_back(int c, int cprev) { 158 159//--------------------------------------------------------------------- 160// pack up and send U(istart) for all j and k 161//--------------------------------------------------------------------- 162 163 int j,k,n,ptr,istart,jp,kp; 164 int phase; 165 int error, buffer_size; 166 167//--------------------------------------------------------------------- 168// Send element 0 to previous processor 169//--------------------------------------------------------------------- 170 istart = 0; 171 buffer_size=MAX_CELL_DIM*MAX_CELL_DIM*BLOCK_SIZE; 172 ptr = 0; 173 for (k = 0; k <= KMAX-1; k++) { 174 for (j = 0; j <= JMAX-1; j++) { 175 for (n = 1; n <= BLOCK_SIZE; n++) { 176 in_buffer(ptr+n) = rhs(n,istart,j,k,cprev); 177 } 178 ptr = ptr+BLOCK_SIZE; 179 } 180 } 181 182//--------------------------------------------------------------------- 183// send and receive buffer 184//--------------------------------------------------------------------- 185 186 for (phase = 0; phase < 3; phase++) { 187 188 if (send_color[WESTDIR]==phase) 189 RCCE_send((char*)in_buffer, buffer_size*sizeof(double), predecessor(1)); 190 if (recv_color[WESTDIR]==phase) 191 RCCE_recv((char*)out_buffer, buffer_size*sizeof(double), successor(1)); 192 } 193 194//--------------------------------------------------------------------- 195// unpack U(isize) for all j and k 196//--------------------------------------------------------------------- 197 198 ptr = 0; 199 for (k = 0; k <= KMAX-1; k++) { 200 for (j = 0; j <= JMAX-1; j++) { 201 for (n = 1; n <= BLOCK_SIZE; n++) { 202 backsub_info(n,j,k,c) = out_buffer(ptr+n); 203 } 204 ptr = ptr+BLOCK_SIZE; 205 } 206 } 207 208 return; 209} 210 211void x_backsubstitute(int first, int last, int c) { 212 213//--------------------------------------------------------------------- 214//--------------------------------------------------------------------- 215 216//--------------------------------------------------------------------- 217// back solve: if last cell, then generate U(isize)=rhs(isize) 218// else assume U(isize) is loaded in un pack backsub_info 219// so just use it 220// after call u(istart) will be sent to next cell 221//--------------------------------------------------------------------- 222 223 int i, j, k; 224 int m,n,isize,jsize,ksize,istart; 225 226 istart = 0; 227 isize = cell_size(1,c)-1; 228 jsize = cell_size(2,c)-end(2,c)-1 ; 229 ksize = cell_size(3,c)-end(3,c)-1; 230 if (last == 0) { 231 for (k = start(3,c); k <= ksize; k++) { 232 for (j = start(2,c); j <= jsize; j++) { 233//--------------------------------------------------------------------- 234// U(isize) uses info from previous cell if not last cell 235//--------------------------------------------------------------------- 236 for (m = 1; m <= BLOCK_SIZE; m++) { 237 for (n = 1; n <= BLOCK_SIZE; n++) { 238 rhs(m,isize,j,k,c) = rhs(m,isize,j,k,c) 239 - lhsc(m,n,isize,j,k,c)* 240 backsub_info(n,j,k,c); 241 } 242 } 243 } 244 } 245 } 246 for (k = start(3,c); k <= ksize; k++) { 247 for (j = start(2,c); j <= jsize; j++) { 248 for (i = isize-1; i >= istart; i--) { 249 for (m = 1; m <= BLOCK_SIZE; m++) { 250 for (n = 1; n <= BLOCK_SIZE; n++) { 251 rhs(m,i,j,k,c) = rhs(m,i,j,k,c) 252 - lhsc(m,n,i,j,k,c)*rhs(n,i+1,j,k,c); 253 } 254 } 255 } 256 } 257 } 258 259 return; 260} 261 262 263//--------------------------------------------------------------------- 264//--------------------------------------------------------------------- 265 266void x_solve_cell(int first, int last, int c) { 267 268//--------------------------------------------------------------------- 269//--------------------------------------------------------------------- 270 271//--------------------------------------------------------------------- 272// performs guaussian elimination on this cell. 273// 274// assumes that unpacking routines for non-first cells 275// preload C' and rhs' from previous cell. 276// 277// assumed send happens outside this routine, but that 278// c'(IMAX) and rhs'(IMAX) will be sent to next cell 279//--------------------------------------------------------------------- 280 281 int i,j,k,isize,ksize,jsize,istart; 282 283 istart = 0; 284 isize = cell_size(1,c)-1; 285 jsize = cell_size(2,c)-end(2,c)-1; 286 ksize = cell_size(3,c)-end(3,c)-1; 287 288 lhsabinit(lhsa, lhsb, isize); 289 290 for (k = start(3,c); k <= ksize; k++) { 291 for (j = start(2,c); j <= jsize; j++) { 292 293//--------------------------------------------------------------------- 294// This function computes the left hand side in the xi-direction 295//--------------------------------------------------------------------- 296 297//--------------------------------------------------------------------- 298// determine a (labeled f) and n jacobians for cell c 299//--------------------------------------------------------------------- 300 for (i = start(1,c)-1; i <= cell_size(1,c) - end(1,c); i++) { 301 302 tmp1 = rho_i(i,j,k,c); 303 tmp2 = tmp1 * tmp1; 304 tmp3 = tmp1 * tmp2; 305//--------------------------------------------------------------------- 306// 307//--------------------------------------------------------------------- 308 fjac(1,1,i) = 0.0e+00; 309 fjac(1,2,i) = 1.0e+00; 310 fjac(1,3,i) = 0.0e+00; 311 fjac(1,4,i) = 0.0e+00; 312 fjac(1,5,i) = 0.0e+00; 313 314 fjac(2,1,i) = -(u(2,i,j,k,c) * tmp2 * 315 u(2,i,j,k,c)) 316 + c2 * qs(i,j,k,c); 317 fjac(2,2,i) = ( 2.0e+00 - c2 ) 318 * ( u(2,i,j,k,c) * tmp1 ); 319 fjac(2,3,i) = - c2 * ( u(3,i,j,k,c) * tmp1 ); 320 fjac(2,4,i) = - c2 * ( u(4,i,j,k,c) * tmp1 ); 321 fjac(2,5,i) = c2; 322 323 fjac(3,1,i) = - ( u(2,i,j,k,c)*u(3,i,j,k,c) ) * tmp2; 324 fjac(3,2,i) = u(3,i,j,k,c) * tmp1; 325 fjac(3,3,i) = u(2,i,j,k,c) * tmp1; 326 fjac(3,4,i) = 0.0e+00; 327 fjac(3,5,i) = 0.0e+00; 328 329 fjac(4,1,i) = - ( u(2,i,j,k,c)*u(4,i,j,k,c) ) * tmp2; 330 fjac(4,2,i) = u(4,i,j,k,c) * tmp1; 331 fjac(4,3,i) = 0.0e+00; 332 fjac(4,4,i) = u(2,i,j,k,c) * tmp1; 333 fjac(4,5,i) = 0.0e+00; 334 335 fjac(5,1,i) = ( c2 * 2.0e0 * qs(i,j,k,c) 336 - c1 * ( u(5,i,j,k,c) * tmp1 ) ) 337 * ( u(2,i,j,k,c) * tmp1 ); 338 fjac(5,2,i) = c1 * u(5,i,j,k,c) * tmp1 339 - c2 340 * ( u(2,i,j,k,c)*u(2,i,j,k,c) * tmp2 341 + qs(i,j,k,c) ); 342 fjac(5,3,i) = - c2 * ( u(3,i,j,k,c)*u(2,i,j,k,c) ) 343 * tmp2; 344 fjac(5,4,i) = - c2 * ( u(4,i,j,k,c)*u(2,i,j,k,c) ) 345 * tmp2; 346 fjac(5,5,i) = c1 * ( u(2,i,j,k,c) * tmp1 ); 347 348 njac(1,1,i) = 0.0e+00; 349 njac(1,2,i) = 0.0e+00; 350 njac(1,3,i) = 0.0e+00; 351 njac(1,4,i) = 0.0e+00; 352 njac(1,5,i) = 0.0e+00; 353 354 njac(2,1,i) = - con43 * c3c4 * tmp2 * u(2,i,j,k,c); 355 njac(2,2,i) = con43 * c3c4 * tmp1; 356 njac(2,3,i) = 0.0e+00; 357 njac(2,4,i) = 0.0e+00; 358 njac(2,5,i) = 0.0e+00; 359 360 njac(3,1,i) = - c3c4 * tmp2 * u(3,i,j,k,c); 361 njac(3,2,i) = 0.0e+00; 362 njac(3,3,i) = c3c4 * tmp1; 363 njac(3,4,i) = 0.0e+00; 364 njac(3,5,i) = 0.0e+00; 365 366 njac(4,1,i) = - c3c4 * tmp2 * u(4,i,j,k,c); 367 njac(4,2,i) = 0.0e+00 ; 368 njac(4,3,i) = 0.0e+00; 369 njac(4,4,i) = c3c4 * tmp1; 370 njac(4,5,i) = 0.0e+00; 371 372 njac(5,1,i) = - ( con43 * c3c4 373 - c1345 ) * tmp3 * SQR(u(2,i,j,k,c)) 374 - ( c3c4 - c1345 ) * tmp3 * SQR(u(3,i,j,k,c)) 375 - ( c3c4 - c1345 ) * tmp3 * SQR(u(4,i,j,k,c)) 376 - c1345 * tmp2 * u(5,i,j,k,c); 377 378 njac(5,2,i) = ( con43 * c3c4 379 - c1345 ) * tmp2 * u(2,i,j,k,c); 380 njac(5,3,i) = ( c3c4 - c1345 ) * tmp2 * u(3,i,j,k,c); 381 njac(5,4,i) = ( c3c4 - c1345 ) * tmp2 * u(4,i,j,k,c); 382 njac(5,5,i) = ( c1345 ) * tmp1; 383 384 } 385//--------------------------------------------------------------------- 386// now jacobians set, so form left hand side in x direction 387//--------------------------------------------------------------------- 388 for (i = start(1,c); i <= isize - end(1,c); i++) { 389 390 tmp1 = dt * tx1; 391 tmp2 = dt * tx2; 392 393 lhsa(1,1,i) = - tmp2 * fjac(1,1,i-1) 394 - tmp1 * njac(1,1,i-1) 395 - tmp1 * dx1 ; 396 lhsa(1,2,i) = - tmp2 * fjac(1,2,i-1) 397 - tmp1 * njac(1,2,i-1); 398 lhsa(1,3,i) = - tmp2 * fjac(1,3,i-1) 399 - tmp1 * njac(1,3,i-1); 400 lhsa(1,4,i) = - tmp2 * fjac(1,4,i-1) 401 - tmp1 * njac(1,4,i-1); 402 lhsa(1,5,i) = - tmp2 * fjac(1,5,i-1) 403 - tmp1 * njac(1,5,i-1); 404 405 lhsa(2,1,i) = - tmp2 * fjac(2,1,i-1) 406 - tmp1 * njac(2,1,i-1); 407 lhsa(2,2,i) = - tmp2 * fjac(2,2,i-1) 408 - tmp1 * njac(2,2,i-1) 409 - tmp1 * dx2; 410 lhsa(2,3,i) = - tmp2 * fjac(2,3,i-1) 411 - tmp1 * njac(2,3,i-1); 412 lhsa(2,4,i) = - tmp2 * fjac(2,4,i-1) 413 - tmp1 * njac(2,4,i-1); 414 lhsa(2,5,i) = - tmp2 * fjac(2,5,i-1) 415 - tmp1 * njac(2,5,i-1); 416 417 lhsa(3,1,i) = - tmp2 * fjac(3,1,i-1) 418 - tmp1 * njac(3,1,i-1); 419 lhsa(3,2,i) = - tmp2 * fjac(3,2,i-1) 420 - tmp1 * njac(3,2,i-1); 421 lhsa(3,3,i) = - tmp2 * fjac(3,3,i-1) 422 - tmp1 * njac(3,3,i-1) 423 - tmp1 * dx3 ; 424 lhsa(3,4,i) = - tmp2 * fjac(3,4,i-1) 425 - tmp1 * njac(3,4,i-1); 426 lhsa(3,5,i) = - tmp2 * fjac(3,5,i-1) 427 - tmp1 * njac(3,5,i-1); 428 429 lhsa(4,1,i) = - tmp2 * fjac(4,1,i-1) 430 - tmp1 * njac(4,1,i-1); 431 lhsa(4,2,i) = - tmp2 * fjac(4,2,i-1) 432 - tmp1 * njac(4,2,i-1); 433 lhsa(4,3,i) = - tmp2 * fjac(4,3,i-1) 434 - tmp1 * njac(4,3,i-1); 435 lhsa(4,4,i) = - tmp2 * fjac(4,4,i-1) 436 - tmp1 * njac(4,4,i-1) 437 - tmp1 * dx4; 438 lhsa(4,5,i) = - tmp2 * fjac(4,5,i-1) 439 - tmp1 * njac(4,5,i-1); 440 441 lhsa(5,1,i) = - tmp2 * fjac(5,1,i-1) 442 - tmp1 * njac(5,1,i-1); 443 lhsa(5,2,i) = - tmp2 * fjac(5,2,i-1) 444 - tmp1 * njac(5,2,i-1); 445 lhsa(5,3,i) = - tmp2 * fjac(5,3,i-1) 446 - tmp1 * njac(5,3,i-1); 447 lhsa(5,4,i) = - tmp2 * fjac(5,4,i-1) 448 - tmp1 * njac(5,4,i-1); 449 lhsa(5,5,i) = - tmp2 * fjac(5,5,i-1) 450 - tmp1 * njac(5,5,i-1) 451 - tmp1 * dx5; 452 453 lhsb(1,1,i) = 1.0e+00 454 + tmp1 * 2.0e+00 * njac(1,1,i) 455 + tmp1 * 2.0e+00 * dx1; 456 lhsb(1,2,i) = tmp1 * 2.0e+00 * njac(1,2,i); 457 lhsb(1,3,i) = tmp1 * 2.0e+00 * njac(1,3,i); 458 lhsb(1,4,i) = tmp1 * 2.0e+00 * njac(1,4,i); 459 lhsb(1,5,i) = tmp1 * 2.0e+00 * njac(1,5,i); 460 461 lhsb(2,1,i) = tmp1 * 2.0e+00 * njac(2,1,i); 462 lhsb(2,2,i) = 1.0e+00 463 + tmp1 * 2.0e+00 * njac(2,2,i) 464 + tmp1 * 2.0e+00 * dx2; 465 lhsb(2,3,i) = tmp1 * 2.0e+00 * njac(2,3,i); 466 lhsb(2,4,i) = tmp1 * 2.0e+00 * njac(2,4,i); 467 lhsb(2,5,i) = tmp1 * 2.0e+00 * njac(2,5,i); 468 469 lhsb(3,1,i) = tmp1 * 2.0e+00 * njac(3,1,i); 470 lhsb(3,2,i) = tmp1 * 2.0e+00 * njac(3,2,i); 471 lhsb(3,3,i) = 1.0e+00 472 + tmp1 * 2.0e+00 * njac(3,3,i) 473 + tmp1 * 2.0e+00 * dx3; 474 lhsb(3,4,i) = tmp1 * 2.0e+00 * njac(3,4,i); 475 lhsb(3,5,i) = tmp1 * 2.0e+00 * njac(3,5,i); 476 477 lhsb(4,1,i) = tmp1 * 2.0e+00 * njac(4,1,i); 478 lhsb(4,2,i) = tmp1 * 2.0e+00 * njac(4,2,i); 479 lhsb(4,3,i) = tmp1 * 2.0e+00 * njac(4,3,i); 480 lhsb(4,4,i) = 1.0e+00 481 + tmp1 * 2.0e+00 * njac(4,4,i) 482 + tmp1 * 2.0e+00 * dx4; 483 lhsb(4,5,i) = tmp1 * 2.0e+00 * njac(4,5,i); 484 485 lhsb(5,1,i) = tmp1 * 2.0e+00 * njac(5,1,i); 486 lhsb(5,2,i) = tmp1 * 2.0e+00 * njac(5,2,i); 487 lhsb(5,3,i) = tmp1 * 2.0e+00 * njac(5,3,i); 488 lhsb(5,4,i) = tmp1 * 2.0e+00 * njac(5,4,i); 489 lhsb(5,5,i) = 1.0e+00 490 + tmp1 * 2.0e+00 * njac(5,5,i) 491 + tmp1 * 2.0e+00 * dx5; 492 493 lhsc(1,1,i,j,k,c) = tmp2 * fjac(1,1,i+1) 494 - tmp1 * njac(1,1,i+1) 495 - tmp1 * dx1; 496 lhsc(1,2,i,j,k,c) = tmp2 * fjac(1,2,i+1) 497 - tmp1 * njac(1,2,i+1); 498 lhsc(1,3,i,j,k,c) = tmp2 * fjac(1,3,i+1) 499 - tmp1 * njac(1,3,i+1); 500 lhsc(1,4,i,j,k,c) = tmp2 * fjac(1,4,i+1) 501 - tmp1 * njac(1,4,i+1); 502 lhsc(1,5,i,j,k,c) = tmp2 * fjac(1,5,i+1) 503 - tmp1 * njac(1,5,i+1); 504 505 lhsc(2,1,i,j,k,c) = tmp2 * fjac(2,1,i+1) 506 - tmp1 * njac(2,1,i+1); 507 lhsc(2,2,i,j,k,c) = tmp2 * fjac(2,2,i+1) 508 - tmp1 * njac(2,2,i+1) 509 - tmp1 * dx2; 510 lhsc(2,3,i,j,k,c) = tmp2 * fjac(2,3,i+1) 511 - tmp1 * njac(2,3,i+1); 512 lhsc(2,4,i,j,k,c) = tmp2 * fjac(2,4,i+1) 513 - tmp1 * njac(2,4,i+1); 514 lhsc(2,5,i,j,k,c) = tmp2 * fjac(2,5,i+1) 515 - tmp1 * njac(2,5,i+1); 516 517 lhsc(3,1,i,j,k,c) = tmp2 * fjac(3,1,i+1) 518 - tmp1 * njac(3,1,i+1); 519 lhsc(3,2,i,j,k,c) = tmp2 * fjac(3,2,i+1) 520 - tmp1 * njac(3,2,i+1); 521 lhsc(3,3,i,j,k,c) = tmp2 * fjac(3,3,i+1) 522 - tmp1 * njac(3,3,i+1) 523 - tmp1 * dx3; 524 lhsc(3,4,i,j,k,c) = tmp2 * fjac(3,4,i+1) 525 - tmp1 * njac(3,4,i+1); 526 lhsc(3,5,i,j,k,c) = tmp2 * fjac(3,5,i+1) 527 - tmp1 * njac(3,5,i+1); 528 529 lhsc(4,1,i,j,k,c) = tmp2 * fjac(4,1,i+1) 530 - tmp1 * njac(4,1,i+1); 531 lhsc(4,2,i,j,k,c) = tmp2 * fjac(4,2,i+1) 532 - tmp1 * njac(4,2,i+1); 533 lhsc(4,3,i,j,k,c) = tmp2 * fjac(4,3,i+1) 534 - tmp1 * njac(4,3,i+1); 535 lhsc(4,4,i,j,k,c) = tmp2 * fjac(4,4,i+1) 536 - tmp1 * njac(4,4,i+1) 537 - tmp1 * dx4; 538 lhsc(4,5,i,j,k,c) = tmp2 * fjac(4,5,i+1) 539 - tmp1 * njac(4,5,i+1); 540 541 lhsc(5,1,i,j,k,c) = tmp2 * fjac(5,1,i+1) 542 - tmp1 * njac(5,1,i+1); 543 lhsc(5,2,i,j,k,c) = tmp2 * fjac(5,2,i+1) 544 - tmp1 * njac(5,2,i+1); 545 lhsc(5,3,i,j,k,c) = tmp2 * fjac(5,3,i+1) 546 - tmp1 * njac(5,3,i+1); 547 lhsc(5,4,i,j,k,c) = tmp2 * fjac(5,4,i+1) 548 - tmp1 * njac(5,4,i+1); 549 lhsc(5,5,i,j,k,c) = tmp2 * fjac(5,5,i+1) 550 - tmp1 * njac(5,5,i+1) 551 - tmp1 * dx5; 552 553 } 554 555 556//--------------------------------------------------------------------- 557// outer most do loops - sweeping in i direction 558//--------------------------------------------------------------------- 559 if (first == 1) { 560 561//--------------------------------------------------------------------- 562// multiply c(istart,j,k) by b_inverse and copy back to c 563// multiply rhs(istart) by b_inverse(istart) and copy to rhs 564//--------------------------------------------------------------------- 565 binvcrhs( &lhsb(1,1,istart), 566 &lhsc(1,1,istart,j,k,c), 567 &rhs(1,istart,j,k,c) ); 568 569 } 570 571//--------------------------------------------------------------------- 572// begin inner most do loop 573// do all the elements of the cell unless last 574//--------------------------------------------------------------------- 575 for (i = istart+first; i <= isize-last; i++) { 576 577//--------------------------------------------------------------------- 578// rhs(i) = rhs(i) - A*rhs(i-1) 579//--------------------------------------------------------------------- 580 matvec_sub(&lhsa(1,1,i), 581 &rhs(1,i-1,j,k,c),&rhs(1,i,j,k,c)); 582 583//--------------------------------------------------------------------- 584// B(i) = B(i) - C(i-1)*A(i) 585//--------------------------------------------------------------------- 586 matmul_sub(&lhsa(1,1,i), 587 &lhsc(1,1,i-1,j,k,c), 588 &lhsb(1,1,i)); 589 590 591//--------------------------------------------------------------------- 592// multiply c(i,j,k) by b_inverse and copy back to c 593// multiply rhs(1,j,k) by b_inverse(1,j,k) and copy to rhs 594//--------------------------------------------------------------------- 595 binvcrhs( &lhsb(1,1,i), 596 &lhsc(1,1,i,j,k,c), 597 &rhs(1,i,j,k,c) ); 598 599 } 600 601//--------------------------------------------------------------------- 602// Now finish up special cases for last cell 603//--------------------------------------------------------------------- 604 if (last == 1) { 605 606//--------------------------------------------------------------------- 607// rhs(isize) = rhs(isize) - A*rhs(isize-1) 608//--------------------------------------------------------------------- 609 matvec_sub(&lhsa(1,1,isize), 610 &rhs(1,isize-1,j,k,c),&rhs(1,isize,j,k,c)); 611 612//--------------------------------------------------------------------- 613// B(isize) = B(isize) - C(isize-1)*A(isize) 614//--------------------------------------------------------------------- 615 matmul_sub(&lhsa(1,1,isize), 616 &lhsc(1,1,isize-1,j,k,c), 617 &lhsb(1,1,isize)); 618 619//--------------------------------------------------------------------- 620// multiply rhs() by b_inverse() and copy to rhs 621//--------------------------------------------------------------------- 622 binvrhs( &lhsb(1,1,isize), 623 &rhs(1,isize,j,k,c) ); 624 625 } 626 } 627 } 628 629 630 return; 631} 632 633