1c--------------------------------------------------------------------- 2c 3c drvftools.f : Bundle of Fortran driver routines for ftools.f 4c 5c Each routine is to be called by les**.c 6c 7c--------------------------------------------------------------------- 8c 9c---------------- 10c drvLesPrepDiag 11c---------------- 12c 13 subroutine drvlesPrepDiag ( flowDiag, ilwork, 14 & iBC, BC, iper, 15 & rowp, colm, 16 & lhsK, lhsP) 17c 18 use pointer_data 19 use pvsQbi 20 use convolImpFlow !brings in the current part of convol coef for imp BC 21 use convolRCRFlow !brings in the current part of convol coef for RCR BC 22 23 include "common.h" 24 include "mpif.h" 25c 26 dimension flowDiag(nshg,4), ilwork(nlwork) 27 dimension iBC(nshg), iper(nshg), BC(nshg,ndofBC) 28 real*8 lhsK(9,nnz_tot), lhsP(4,nnz_tot) 29 integer rowp(nnz_tot), colm(nshg+1) 30 integer n, k 31c 32 integer sparseloc 33c 34c 35c.... Clear the flowdiag 36c 37 if((flmpl.eq.1).or.(ipord.gt.1)) then 38 do n = 1, nshg 39 k = sparseloc( rowp(colm(n)), colm(n+1)-colm(n), n ) 40 & + colm(n)-1 41c 42 flowdiag(n,1) = lhsK(1,k) 43 flowdiag(n,2) = lhsK(5,k) 44 flowdiag(n,3) = lhsK(9,k) 45c 46 flowdiag(n,4) = lhsP(4,k) 47 enddo 48 else 49 flowDiag = zero 50 do n = 1, nshg ! rowsum put on the diagonal instead of diag entry 51 do k=colm(n),colm(n+1)-1 52 53c 54 flowdiag(n,1) = flowdiag(n,1) + abs(lhsK(1,k)) 55c & + lhsK(2,k) + lhsK(3,k) 56 flowdiag(n,2) = flowdiag(n,2) + abs(lhsK(5,k)) 57c & + lhsK(4,k) + lhsK(6,k) 58 flowdiag(n,3) = flowdiag(n,3) + abs(lhsK(9,k)) 59c & + lhsK(7,k) + lhsK(8,k) 60c 61 flowdiag(n,4) = flowdiag(n,4) + abs(lhsP(4,k)) 62 enddo 63 flowdiag(n,:)=flowdiag(n,:)*pt33 64 enddo 65 endif 66 if(ipvsq.ge.3) then ! for first cut only do diagonal extraction 67 ! this is not yet correct for multi procs I suspect if partition 68 ! boundary cuts a p=QR face 69 tfact=alfi * gami * Delt(1) 70 do n=1,nshg 71 if(numResistSrfs.gt.zero) then 72 do k = 1,numResistSrfs 73 if (nsrflistResist(k).eq.ndsurf(n)) then 74 irankCoupled=k 75 flowdiag(n,1:3) = flowdiag(n,1:3) 76 & + tfact*ValueListResist(irankCoupled)* 77 & NABI(n,:)*NABI(n,:) 78 exit 79 endif 80 enddo 81 elseif(numImpSrfs.gt.zero) then 82 do k = 1,numImpSrfs 83 if (nsrflistImp(k).eq.ndsurf(n)) then 84 irankCoupled=k 85 flowdiag(n,1:3) = flowdiag(n,1:3) 86 & + tfact*ImpConvCoef(ntimeptpT+2,irankCoupled)* 87 & NABI(n,:)*NABI(n,:) 88 exit 89 endif 90 enddo 91 elseif(numRCRSrfs.gt.zero) then 92 do k = 1,numRCRSrfs 93 if (nsrflistRCR(k).eq.ndsurf(n)) then 94 irankCoupled=k 95 flowdiag(n,1:3) = flowdiag(n,1:3) 96 & + tfact*RCRConvCoef(lstep+2,irankCoupled)* !check lstep+2 if restart from t.ne.0 97 & NABI(n,:)*NABI(n,:) 98 exit 99 endif 100 enddo 101 endif 102 enddo 103 endif 104c 105 106c 107 if(iabc==1) !are there any axisym bc's 108 & call rotabc(flowdiag, iBC, 'in ') 109c 110 111c 112c.... communicate : add the slaves part to the master's part of flowDiag 113c 114 if (numpe > 1) then 115 call commu (flowDiag, ilwork, nflow, 'in ') 116 endif 117c 118c.... satisfy the boundary conditions on the diagonal 119c 120 call bc3diag(iBC, BC, flowDiag) 121c 122c 123c.... on processor periodicity was not taken care of in the setting of the 124c boundary conditions on the matrix. Take care of it now. 125c 126 call bc3per(iBC, flowDiag, iper, ilwork, 4) 127c 128c... slaves and masters have the correct values 129c 130c 131c.... Calculate square root 132c 133 do i = 1, nshg 134 do j = 1, nflow 135 if (flowDiag(i,j).ne.0) 136 & flowDiag(i,j) = 1. / sqrt(abs(flowDiag(i,j))) 137 enddo 138 enddo 139c 140 return 141 end 142 143c 144c------------- 145c drvsclrDiag 146c------------- 147c 148 subroutine drvsclrDiag ( sclrDiag, ilwork, iBC, BC, iper, 149 & rowp, colm, lhsS ) 150c 151 use pointer_data 152 include "common.h" 153 include "mpif.h" 154c 155 integer ilwork(nlwork), iBC(nshg), iper(nshg), 156 & rowp(nnz_tot), colm(nshg+1) 157 158 real*8 sclrDiag(nshg), lhsS(nnz_tot), BC(nshg,ndofBC) 159 integer sparseloc 160 161 sclrDiag = zero 162 do n = 1, nshg 163 k = sparseloc( rowp(colm(n)), colm(n+1)-colm(n), n ) 164 & + colm(n)-1 165c 166 sclrDiag(n) = lhsS(k) 167 enddo 168c 169c.... communicate : add the slaves part to the master's part of sclrDiag 170c 171 if (numpe > 1) then 172 call commu (sclrDiag, ilwork, 1, 'in ') 173 endif 174c 175c.... satisfy the boundary conditions on the diagonal 176c 177 call bc3SclrDiag(iBC, sclrDiag) 178c 179c 180c.... on processor periodicity was not taken care of in the setting of the 181c boundary conditions on the matrix. Take care of it now. 182c 183 call bc3per(iBC, sclrDiag, iper, ilwork, 1) 184c 185c... slaves and masters have the correct values 186c 187c 188c.... Calculate square root 189c 190 do i = 1, nshg 191 if (sclrDiag(i).ne.0) then 192 sclrDiag(i) = 1. / sqrt(abs(sclrDiag(i))) 193 endif 194 enddo 195c 196 return 197 end 198 199C============================================================================ 200C 201C "fLesSparseApG": 202C 203C============================================================================ 204 subroutine fLesSparseApG( col, row, pLhs, 205 & p, q, nNodes, 206 & nnz_tot ) 207c 208c.... Data declaration 209c 210 implicit none 211 integer nNodes, nnz_tot 212 integer col(nNodes+1), row(nnz_tot) 213 real*8 pLhs(4,nnz_tot), p(nNodes), q(nNodes,3) 214c 215 real*8 pisave 216 integer i, j, k 217c 218c.... clear the vector 219c 220 do i = 1, nNodes 221 q(i,1) = 0 222 q(i,2) = 0 223 q(i,3) = 0 224 enddo 225c 226c.... Do an AP product 227c 228 do i = 1, nNodes 229c 230 pisave = p(i) 231cdir$ ivdep 232 do k = col(i), col(i+1)-1 233 j = row(k) 234c 235 q(j,1) = q(j,1) - pLhs(1,k) * pisave 236 q(j,2) = q(j,2) - pLhs(2,k) * pisave 237 q(j,3) = q(j,3) - pLhs(3,k) * pisave 238 enddo 239 enddo 240c 241c.... end 242c 243 return 244 end 245 246C============================================================================ 247C 248C "fLesSparseApKG": 249C 250C============================================================================ 251 252 subroutine fLesSparseApKG( col, row, kLhs, pLhs, 253 1 p, q, nNodes, 254 2 nnz_tot_hide ) 255c 256c.... Data declaration 257c 258c implicit none 259 use pvsQbi 260 include "common.h" 261 integer nNodes 262 integer col(nNodes+1), row(nnz_tot) 263 real*8 kLhs(9,nnz_tot), pLhs(4,nnz_tot) 264 real*8 p(nNodes,4), q(nNodes,3) 265c 266 real*8 tmp1, tmp2, tmp3, pisave 267 integer i, j, k 268c 269c.... clear the vector 270c 271 do i = 1, nNodes 272 q(i,1) = 0 273 q(i,2) = 0 274 q(i,3) = 0 275 enddo 276c 277c.... Do an AP product 278c 279 do i = 1, nNodes 280c 281 tmp1 = 0 282 tmp2 = 0 283 tmp3 = 0 284 pisave = p(i,4) 285cdir$ ivdep 286 do k = col(i), col(i+1)-1 287 j = row(k) 288 tmp1 = tmp1 289 1 + kLhs(1,k) * p(j,1) 290 2 + kLhs(4,k) * p(j,2) 291 3 + kLhs(7,k) * p(j,3) 292 tmp2 = tmp2 293 1 + kLhs(2,k) * p(j,1) 294 2 + kLhs(5,k) * p(j,2) 295 3 + kLhs(8,k) * p(j,3) 296 tmp3 = tmp3 297 1 + kLhs(3,k) * p(j,1) 298 2 + kLhs(6,k) * p(j,2) 299 3 + kLhs(9,k) * p(j,3) 300c 301 q(j,1) = q(j,1) - pLhs(1,k) * pisave 302 q(j,2) = q(j,2) - pLhs(2,k) * pisave 303 q(j,3) = q(j,3) - pLhs(3,k) * pisave 304 enddo 305 q(i,1) = q(i,1) + tmp1 306 q(i,2) = q(i,2) + tmp2 307 q(i,3) = q(i,3) + tmp3 308 enddo 309 310 if(ipvsq.ge.2) then 311 tfact=alfi * gami * Delt(1) 312 call ElmpvsQ(q,p,tfact) 313 endif 314c 315c.... end 316c 317 return 318 end 319 320 321C============================================================================ 322C 323C "fLesSparseApNGt": 324C 325C============================================================================ 326 327 subroutine fLesSparseApNGt( col, row, pLhs, 328 1 p, q, nNodes, 329 2 nnz_tot ) 330c 331c.... Data declaration 332c 333 implicit none 334 integer nNodes, nnz_tot 335 integer col(nNodes+1), row(nnz_tot) 336 real*8 pLhs(4,nnz_tot), p(nNodes,3), q(nNodes) 337c 338 real*8 tmp 339 integer i, j, k 340c 341c.... Do an AP product 342c 343 do i = nNodes, 1, -1 344c 345 tmp = 0 346 do k = col(i), col(i+1)-1 347 j = row(k) 348c 349 tmp = tmp 350 1 + pLhs(1,k) * p(j,1) 351 2 + pLhs(2,k) * p(j,2) 352 3 + pLhs(3,k) * p(j,3) 353 enddo 354 q(i) = tmp 355 enddo 356c 357c.... end 358c 359 return 360 end 361 362C============================================================================ 363C 364C "fLesSparseApNGtC": 365C 366C============================================================================ 367 368 subroutine fLesSparseApNGtC( col, row, pLhs, 369 1 p, q, nNodes, 370 2 nnz_tot ) 371c 372c.... Data declaration 373c 374 implicit none 375 integer nNodes, nnz_tot 376 integer col(nNodes+1), row(nnz_tot) 377 real*8 pLhs(4,nnz_tot), p(nNodes,4), q(nNodes) 378c 379 real*8 tmp 380 integer i, j, k 381c 382c.... Do an AP product 383c 384 do i = nNodes, 1, -1 385c 386 tmp = 0 387 do k = col(i), col(i+1)-1 388 j = row(k) 389c 390 tmp = tmp 391 1 + pLhs(1,k) * p(j,1) 392 2 + pLhs(2,k) * p(j,2) 393 3 + pLhs(3,k) * p(j,3) 394 4 + pLhs(4,k) * p(j,4) 395 enddo 396 q(i) = tmp 397 enddo 398c 399c.... end 400c 401 return 402 end 403 404C============================================================================ 405C 406C "fLesSparseApFull": 407C 408C============================================================================ 409 410 subroutine fLesSparseApFull( col, row, kLhs, pLhs, 411 1 p, q, nNodes, 412 2 nnz_tot_hide ) 413c 414c.... Data declaration 415c 416c implicit none 417 use pvsQbi 418 include "common.h" 419 420 integer nNodes, nnz_tot 421 integer col(nNodes+1), row(nnz_tot) 422 real*8 kLhs(9,nnz_tot), pLhs(4,nnz_tot) 423 real*8 p(nNodes,4), q(nNodes,4) 424c 425 real*8 tmp1, tmp2, tmp3, tmp4, pisave 426 integer i, j, k 427c 428c.... clear the vector 429c 430 do i = 1, nNodes 431 q(i,1) = 0 432 q(i,2) = 0 433 q(i,3) = 0 434 enddo 435c 436c.... Do an AP product 437c 438 do i = 1, nNodes 439c 440 tmp1 = 0 441 tmp2 = 0 442 tmp3 = 0 443 tmp4 = 0 444 pisave = p(i,4) 445cdir$ ivdep 446 do k = col(i), col(i+1)-1 447 j = row(k) 448c 449 tmp1 = tmp1 450 1 + kLhs(1,k) * p(j,1) 451 2 + kLhs(4,k) * p(j,2) 452 3 + kLhs(7,k) * p(j,3) 453 tmp2 = tmp2 454 1 + kLhs(2,k) * p(j,1) 455 2 + kLhs(5,k) * p(j,2) 456 3 + kLhs(8,k) * p(j,3) 457 tmp3 = tmp3 458 1 + kLhs(3,k) * p(j,1) 459 2 + kLhs(6,k) * p(j,2) 460 3 + kLhs(9,k) * p(j,3) 461c 462 tmp4 = tmp4 463 1 + pLhs(1,k) * p(j,1) 464 2 + pLhs(2,k) * p(j,2) 465 3 + pLhs(3,k) * p(j,3) 466 4 + pLhs(4,k) * p(j,4) 467c 468 q(j,1) = q(j,1) - pLhs(1,k) * pisave 469 q(j,2) = q(j,2) - pLhs(2,k) * pisave 470 q(j,3) = q(j,3) - pLhs(3,k) * pisave 471 enddo 472 q(i,1) = q(i,1) + tmp1 473 q(i,2) = q(i,2) + tmp2 474 q(i,3) = q(i,3) + tmp3 475 q(i,4) = tmp4 476 enddo 477 if(ipvsq.ge.2) then 478 tfact=alfi * gami * Delt(1) 479 call ElmpvsQ(q,p,tfact) 480 endif 481c 482c.... end 483c 484 return 485 end 486 487C============================================================================ 488C 489C "fLesSparseApSclr": 490C 491C============================================================================ 492 493 subroutine fLesSparseApSclr( col, row, lhs, 494 1 p, q, nNodes, 495 & nnz_tot) 496c 497c.... Data declaration 498c 499 implicit none 500 integer nNodes, nnz_tot 501 integer col(nNodes+1), row(nnz_tot) 502 real*8 lhs(nnz_tot), p(nNodes), q(nNodes) 503c 504 real*8 tmp 505 integer i, j, k 506c 507c.... Do an AP product 508c 509 do i = nNodes, 1, -1 510c 511 tmp = 0 512 do k = col(i), col(i+1)-1 513 tmp = tmp + lhs(k) * p(row(k)) 514 enddo 515 q(i) = tmp 516 enddo 517c 518c.... end 519c 520 return 521 end 522C============================================================================ 523 subroutine commOut( global, ilwork, n, 524 & iper, iBC, BC ) 525 526 include "common.h" 527 528 real*8 global(nshg,n), BC(nshg,ndofBC) 529 integer ilwork(nlwork), iper(nshg), iBC(nshg) 530c 531 if ( numpe .gt. 1) then 532 call commu ( global, ilwork, n, 'out') 533 endif 534c 535c before doing AP product P must be made periodic 536c on processor slaves did not get updated with the 537c commu (out) so do it here 538c 539 do i=1,n 540 global(:,i) = global(iper(:),i) ! iper(i)=i if non-slave so no danger 541 enddo 542c 543c slave has masters value, for abc we need to rotate it 544c (if this is a vector only no SCALARS) 545 if((iabc==1) .and. (n.gt.1)) !are there any axisym bc's 546 & call rotabc(global, iBC, 'out') 547 548 549c$$$ do j = 1,nshg 550c$$$ if (btest(iBC(j),10)) then 551c$$$ i = iper(j) 552c$$$ res(j,:) = res(i,:) 553c$$$ endif 554c$$$ enddo 555 556 return 557 end 558 559C============================================================================ 560 subroutine commIn( global, ilwork, n, 561 & iper, iBC, BC ) 562 563 include "common.h" 564 565 real*8 global(nshg,n), BC(nshg,ndofBC) 566 integer ilwork(nlwork), iper(nshg), iBC(nshg) 567c 568 if((iabc==1) .and. (n.gt.1)) !are there any axisym bc's 569 & call rotabc(global, iBC, 'in ') 570c 571 572 if ( numpe .gt. 1 ) then 573 call commu ( global, ilwork, n, 'in ') 574 endif 575 576 call bc3per ( iBC, global, iper, ilwork, n) 577 578 return 579 end 580 581