1 /*====================================================================== 2 * 3 * lestools.c : Linear Algebra Solver Tools 4 * 5 * small single character : vector or matrix 6 * 7 *====================================================================== 8 */ 9 #include "les.h" 10 #include "usr.h" 11 #ifdef intel 12 void DRVSCLRDIAG( double *sclrDiag, int *ilwork, int *iBC, 13 double *BC, int *iper, int *rowp, 14 int *colm, double *lhsS); 15 16 void FMTXBLKDAXPY(double *srcpnt, double *dstpnt, double *coef, 17 int *mDims, int *dim); 18 19 void FMTXBLKDYEAX(double *srcpnt, double *dstpnt, double *coef, 20 int *mDims, int *dim); 21 22 void FMTXBLKDMAXPY(double *srcpnt, double *dstpnt, double *coef, 23 int *mDims, int *dim); 24 25 void FMTXVDIMVECCP(double *srcpnt, double *dstpnt, int *nSrcDims, 26 int *nDstDims, int *nDims, int *nNodes ); 27 28 void DRVLESPREPDIAG( double *flowDiag, int *ilwork, int *iBC, 29 double *BC, int *iper, int *rowp, 30 int *colm, double *lhsK, double *lhsP) ; 31 32 void FMTXVDIMVECMULT( double* , double*, 33 double *dstpnt, int *nSrcDims, 34 int *nDofs, int *nDstDims, 35 int *nDims, int *nNodes); 36 37 void FLESZERO( double *dstpnt, int *nDims, int *nNodes); 38 39 void FLESCP( double *srcpnt, double *dstpnt, int *nDims, 40 int *dim ) ; 41 42 void FLESSCALE( double *dstpnt, 43 double *coef, 44 int *nDims, 45 int *dim ) ; 46 47 void FLESSCALECP( double *srcpnt, 48 double *dstpnt, 49 double *coef, 50 int *nDims, 51 int *dim ) ; 52 53 void FLESADD( double *srcpnt, 54 double *dstpnt, 55 int *nDims, 56 int *dim ) ; 57 58 void FLESSUB( double *srcpnt, 59 double *dstpnt, 60 int *nDims, 61 int *dim ) ; 62 63 void DRVALLREDUCESCLR(double *tmpp, double *tmp); 64 65 void DRVALLREDUCE(double *valuesp, double *values, int* mDims); 66 67 double FLESDOT1( double *srcpnt, int* nDims, int* dim ); 68 69 double FLESDOT2(double *src1pnt, double *src2pnt, int* nDims, int* dim ); 70 71 void FLESDAXPY( double *srcpnt, double *lesP, double *sclrRegFct, 72 int *nDstDims, int *nNodes); 73 74 void FLESDXPAY( double *srcpnt, double *dstpnt, double *coef, 75 int *nDims, int *dim); 76 77 void FLESINV( double *dstpnt, int *nDims, int *dim ) ; 78 79 80 81 void FMTXBLKDOT2 (double *src1pnt, double *src2pnt, double *valuesp, 82 int* mDims, int* dim); 83 84 void COMMIN( double *lesQ, int *ilwork, int *nPs, 85 int *iper, int *iBC, double *BC); 86 87 void COMMOUT( double *lesP, int *ilwork, int *nQs, 88 int *iper, int *iBC, double *BC); 89 90 void FLESSPARSEAPFULL(int *colm, int *rowp, double *lhsK, 91 double *lhsP, double *lesP, double *lesQ, 92 int *nNodes, int *nnz); 93 94 void FLESSPARSEAPSCLR(int *colm, int *rowp, 95 double *lhsS, double *lesP, double *lesQ, 96 int *nNodes, int *nnz); 97 98 void FMTXVDIMVECDOT2 (double *src1pnt, double *src2pnt, double *coefp, 99 int *nSrc1Dims,int *nSrc2Dims, int *nDims, int *nNodes); 100 101 void FMTXVDIMVECDAXPY ( double *srcpnt, double *dstpnt, double *coef, 102 int *nSrcDims, int *nDstDims, int *nDims, 103 int *nNodes); 104 105 void FLESSPARSEAPG ( int *colm, int *rowp, 106 double *lhsP, double *lesP, double *lesQ, 107 int *nNodes, int *nnz); 108 109 void FLESSPARSEAPNGT ( int *colm, int *rowp, 110 double *lhsP, double *lesP, double *lesQ, 111 int *nNodes, int *nnz); 112 113 void FLESSPARSEAPNGTC ( int *colm, int *rowp, 114 double *lhsP, double *lesP, double *lesQ, 115 int *nNodes, int *nnz); 116 117 void FLESSPARSEAPKG ( int *colm, int *rowp, double *lhsK, 118 double *lhsP, double *lesP, double *lesQ, 119 int *nNodes, int *nnz); 120 void RAMG_INTERFACE ( int *colm, int *rowp, double *lhsK,double *lhsP, 121 double *mcgR,double *mcgZ, 122 int *ilwork, double *BC, int *iBC,int *iper); 123 #endif 124 /*---------------------------------------------------------------------- 125 * 126 * lesPreDiag 127 * 128 * operation : a = 1/sqrt(abs(a)) 129 * 130 *---------------------------------------------------------------------- 131 */ 132 void lesPrepDiag( UsrHd usrHd ) 133 { 134 char* funcName = "lesPrepDiag" ; /* function name */ 135 136 if ( (usrHd->eqnType) == 1 ) { /* continuity & momentum */ 137 138 drvlesPrepDiag( usrHd->flowDiag, 139 usrHd->ilwork, usrHd->iBC, 140 usrHd->BC, usrHd->iper, 141 usrHd->rowp, usrHd->colm, 142 usrHd->lhsK, usrHd->lhsP) ; 143 } 144 145 if ( (usrHd->eqnType) == 2 ) { /* temperature or scalar variable */ 146 147 drvsclrDiag ( usrHd->sclrDiag, usrHd->ilwork, 148 usrHd->iBC, usrHd->BC, 149 usrHd->iper, usrHd->rowp, 150 usrHd->colm, usrHd->lhsS ) ; 151 152 } 153 154 } 155 156 /*---------------------------------------------------------------------- 157 * 158 * lesDiagScaleCp 159 * 160 * operation : c = a * b 161 * 162 *---------------------------------------------------------------------- 163 */ 164 void lesDiagScaleCp ( UsrHd usrHd, 165 Integer srcId, 166 Integer dstId, 167 Integer nSrcDims, 168 Integer srcOff, 169 Integer nDstDims, 170 Integer dstOff, 171 Integer diagOff, 172 Integer nDims ) 173 { 174 char* funcName = "lesDiagScaleCp" ; /* function name */ 175 Integer nDofs ; /* No. of Dofs */ 176 Real* dstpnt ; /* destination */ 177 Real* srcpnt ; /* source */ 178 179 if ( (usrHd->eqnType) == 1 ) { 180 181 nDofs = 4 ; 182 183 srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; 184 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; 185 186 fMtxVdimVecMult( srcpnt, 187 usrHd->flowDiag+diagOff*usrHd->nNodes, 188 dstpnt, 189 &nSrcDims, 190 &nDofs, 191 &nDstDims, 192 &nDims, 193 &(usrHd->nNodes) ) ; 194 } 195 196 if ( (usrHd->eqnType) == 2 ) { 197 198 nDofs = 1 ; 199 200 srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; 201 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; 202 203 fMtxVdimVecMult( srcpnt, 204 usrHd->sclrDiag+diagOff*usrHd->nNodes, 205 dstpnt, 206 &nSrcDims, 207 &nDofs, 208 &nDstDims, 209 &nDims, 210 &(usrHd->nNodes) ) ; 211 } 212 213 } 214 215 /*---------------------------------------------------------------------- 216 * 217 * lesZero 218 * 219 * operation : a = 0 220 * 221 *---------------------------------------------------------------------- 222 */ 223 void lesZero ( UsrHd usrHd, 224 Integer dstId, 225 Integer nDims ) 226 { 227 char* funcName = "lesZero" ; /* function namea */ 228 Real* dstpnt ; /* destination */ 229 Integer dstOff ; /* destination offset */ 230 231 232 dstOff = 0 ; 233 234 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ); 235 236 flesZero( dstpnt, &nDims, &(usrHd->nNodes) ); 237 } 238 239 /*---------------------------------------------------------------------- 240 * 241 * lesCp 242 * 243 * operation : b = a 244 * 245 *----------------------------------------------------------------------- 246 */ 247 void lesCp ( UsrHd usrHd, 248 Integer srcId, 249 Integer dstId, 250 Integer nDims ) 251 { 252 char* funcName = "lesCp" ; /* function name */ 253 Real* srcpnt ; /* source */ 254 Real* dstpnt ; /* destination */ 255 Integer dim ; /* a simple dimension */ 256 Integer srcOff ; /* source offset */ 257 Integer dstOff ; /* destination offset */ 258 259 srcOff = 0 ; 260 dstOff = 0 ; 261 262 srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; 263 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; 264 265 dim = usrHd->nNodes ; 266 267 flesCp( srcpnt, dstpnt, &nDims, &dim ) ; 268 } 269 270 /*----------------------------------------------------------------------- 271 * 272 * lesScale 273 * 274 * operation : a = a * coef 275 * 276 *----------------------------------------------------------------------- 277 */ 278 void lesScale ( UsrHd usrHd, 279 Integer dstId, 280 Real coef, 281 Integer nDims ) 282 { 283 char* funcName = "lesScale" ; /* function name */ 284 Real* dstpnt ; /* destination */ 285 Integer dstOff ; /* destination offset */ 286 Integer dim ; /* a simple dimension */ 287 288 dstOff = 0 ; 289 290 dim = usrHd->nNodes ; 291 292 dstpnt = usrPointer( usrHd, dstId, dstOff, nDims ) ; 293 294 flesScale( dstpnt, 295 &coef, 296 &nDims, 297 &dim ) ; 298 } 299 /*----------------------------------------------------------------------- 300 * 301 * lesScaleCp 302 * 303 * operation : b = a * coef 304 * 305 *----------------------------------------------------------------------- 306 */ 307 void lesScaleCp ( UsrHd usrHd, 308 Integer srcId, 309 Integer dstId, 310 Real coef, 311 Integer nDims ) 312 { 313 char* funcName = "lesScaleCp" ; /* function name */ 314 Real* srcpnt ; /* source */ 315 Real* dstpnt ; /* destination */ 316 Integer dim ; /* a simple dimension */ 317 Integer srcOff ; /* source offset */ 318 Integer dstOff ; /* destination offset */ 319 320 srcOff = 0 ; 321 dstOff = 0 ; 322 323 srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; 324 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; 325 326 dim = usrHd->nNodes ; 327 328 flesScaleCp( srcpnt, 329 dstpnt, 330 &coef, 331 &nDims, 332 &dim ) ; 333 } 334 335 /*----------------------------------------------------------------------- 336 * 337 * lesAdd 338 * 339 * operation : b = b + a 340 * 341 *----------------------------------------------------------------------- 342 */ 343 void lesAdd ( UsrHd usrHd, 344 Integer srcId, 345 Integer dstId, 346 Integer nDims ) 347 { 348 char* funcName = "lesAdd" ; /* function name */ 349 Real* srcpnt ; /* source */ 350 Real* dstpnt ; /* destination */ 351 Integer srcOff ; /* source offset */ 352 Integer dstOff ; /* destination offset */ 353 Integer dim ; /* a simple dimension */ 354 355 srcOff = 0 ; 356 dstOff = 0 ; 357 358 srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; 359 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; 360 361 dim = usrHd->nNodes ; 362 363 flesAdd( srcpnt, 364 dstpnt, 365 &nDims, 366 &dim ) ; 367 } 368 369 /*----------------------------------------------------------------------- 370 * 371 * lesSub 372 * 373 * operation : b = b - a 374 * 375 *----------------------------------------------------------------------- 376 */ 377 void lesSub ( UsrHd usrHd, 378 Integer srcId, 379 Integer dstId, 380 Integer nDims ) 381 { 382 char* funcName = "lesSub" ; /* function name */ 383 Real* srcpnt ; /* source */ 384 Real* dstpnt ; /* destination */ 385 Integer srcOff ; /* source offset */ 386 Integer dstOff ; /* destination offset */ 387 Integer dim ; /* a simple dimension */ 388 389 srcOff = 0 ; 390 dstOff = 0 ; 391 392 srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; 393 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; 394 395 dim = usrHd->nNodes ; 396 397 flesSub( srcpnt, 398 dstpnt, 399 &nDims, 400 &dim ) ; 401 } 402 403 /*----------------------------------------------------------------------- 404 * 405 * lesDot1 406 * 407 * operation : tmp = tmp + a * a 408 * 409 *----------------------------------------------------------------------- 410 */ 411 Real lesDot1 ( UsrHd usrHd, 412 Integer srcId, 413 Integer nDims ) 414 { 415 char* funcName = "lesDot1" ; /* function name */ 416 Real* srcpnt ; /* source */ 417 Integer srcOff ; /* source offset */ 418 Integer dim ; /* a simple dimension */ 419 Real tmp ; /* a temporary value */ 420 Real tmpp ; /* a temporary value on each proc. */ 421 422 srcOff = 0 ; 423 424 srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; 425 426 dim = usrHd->nNodes ; 427 428 tmpp = flesDot1( srcpnt, 429 &nDims, 430 &dim ) ; 431 432 drvAllreducesclr ( &tmpp, 433 &tmp ) ; 434 #ifdef AMG 435 ramg_normcheck(&tmp); 436 #endif 437 return ( tmp ) ; 438 } 439 440 /*----------------------------------------------------------------------- 441 * 442 * lesDot2 443 * 444 * operation : tmp = tmp + a * b 445 * 446 *----------------------------------------------------------------------- 447 */ 448 Real lesDot2 ( UsrHd usrHd, 449 Integer src1Id, 450 Integer src2Id, 451 Integer nDims ) 452 { 453 char* funcName = "lesDot2" ; /* function name */ 454 Real* src1pnt ; /* source 1 */ 455 Real* src2pnt ; /* source 2 */ 456 Integer src1Off ; /* source 1 offset */ 457 Integer src2Off ; /* source 2 offset */ 458 Integer dim ; /* a simple dimension */ 459 Real tmp ; /* a temporary value */ 460 Real tmpp ; /* a temporary value on each proc. */ 461 462 src1Off = 0 ; 463 src2Off = 0 ; 464 465 src1pnt = usrPointer ( usrHd, src1Id, src1Off, nDims ); 466 src2pnt = usrPointer ( usrHd, src2Id, src2Off, nDims ); 467 468 dim = usrHd->nNodes ; 469 470 tmpp = flesDot2( src1pnt, 471 src2pnt, 472 &nDims, 473 &dim ) ; 474 475 drvAllreducesclr ( &tmpp, 476 &tmp ) ; 477 478 return ( tmp ) ; 479 } 480 481 /*----------------------------------------------------------------------- 482 * 483 * lesDaxpy 484 * 485 * operation : y = y + coef * x 486 * 487 *----------------------------------------------------------------------- 488 */ 489 void lesDaxpy ( UsrHd usrHd, 490 Integer srcId, 491 Integer dstId, 492 Real coef, 493 Integer nDims ) 494 { 495 char* funcName = "lesDapxy" ; /* function name */ 496 Real* srcpnt ; /* source */ 497 Real* dstpnt ; /* destination */ 498 Integer srcOff ; /* source offset */ 499 Integer dstOff ; /* destination offset */ 500 Integer dim ; /* a simple dimension */ 501 502 srcOff = 0 ; 503 dstOff = 0 ; 504 505 srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; 506 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; 507 508 dim = usrHd->nNodes ; 509 510 flesDaxpy( srcpnt, 511 dstpnt, 512 &coef, 513 &nDims, 514 &dim ) ; 515 } 516 517 /*----------------------------------------------------------------------- 518 * 519 * lesDxpay 520 * 521 * operation : y = coef * y + x 522 * 523 *----------------------------------------------------------------------- 524 */ 525 void lesDxpay ( UsrHd usrHd, 526 Integer srcId, 527 Integer dstId, 528 Real coef, 529 Integer nDims ) 530 { 531 char* funcName = "lesDxpay" ; /* function name */ 532 Real* srcpnt ; /* source */ 533 Real* dstpnt ; /* destination */ 534 Integer srcOff ; /* source offset */ 535 Integer dstOff ; /* destination offset */ 536 Integer dim ; /* a simple dimension */ 537 538 srcOff = 0 ; 539 dstOff = 0 ; 540 541 srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; 542 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; 543 544 dim = usrHd->nNodes ; 545 546 flesDxpay( srcpnt, 547 dstpnt, 548 &coef, 549 &nDims, 550 &dim ) ; 551 } 552 553 /*----------------------------------------------------------------------- 554 * 555 * lesInv 556 * 557 * operation : x = 1. / x 558 * 559 *----------------------------------------------------------------------- 560 */ 561 void lesInv ( UsrHd usrHd, 562 Integer dstId, 563 Integer nDims ) 564 { 565 char* funcName = "lesInv" ; /* function name */ 566 Integer dim ; /* a simple dimension */ 567 Real* dstpnt ; /* destination */ 568 Integer dstOff ; /* destination offset */ 569 570 dstOff = 0 ; 571 572 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; 573 574 dim = usrHd->nNodes ; 575 576 flesInv( dstpnt, &nDims, &dim ) ; 577 } 578 579 /*----------------------------------------------------------------------- 580 * 581 * lesBlkDot2 582 * 583 * operation : 584 * 585 *----------------------------------------------------------------------- 586 */ 587 void lesBlkDot2 ( UsrHd usrHd, 588 Integer src1Id, 589 Integer src2Id, 590 Real* values, 591 Integer mDims, 592 Integer nDims ) 593 { 594 char* funcName = "lesBlkDot2" ; /* function name */ 595 Real* src1pnt ; /* source 1 */ 596 Real* src2pnt ; /* source 2 */ 597 Integer src1Off ; /* source 1 offset */ 598 Integer src2Off ; /* source 2 offset */ 599 Integer dim ; /* a simple dimension */ 600 Real* valuesp ; /* temporary values on proc */ 601 602 if ( mDims * nDims == 0 ) { 603 return ; 604 } 605 606 valuesp = (double *) malloc (mDims * sizeof(double)) ; 607 608 src1Off = 0 ; 609 src2Off = 0 ; 610 611 src1pnt = usrPointer ( usrHd, src1Id, src1Off, nDims ) ; 612 src2pnt = usrPointer ( usrHd, src2Id, src2Off, nDims ) ; 613 614 dim = nDims * usrHd->nNodes ; 615 616 fMtxBlkDot2( src1pnt, 617 src2pnt, 618 valuesp, 619 &mDims, 620 &dim ) ; 621 622 drvAllreduce ( valuesp, 623 values, 624 &mDims ) ; 625 626 free(valuesp); 627 } 628 629 /*----------------------------------------------------------------------- 630 * 631 * lesBlkDaxpy 632 * 633 * operation : 634 * 635 *----------------------------------------------------------------------- 636 */ 637 void lesBlkDaxpy ( UsrHd usrHd, 638 Integer srcId, 639 Integer dstId, 640 Real* coef, 641 Integer mDims, 642 Integer nDims ) 643 { 644 char* funcName = "lesBlkDaxpy" ; /* function name */ 645 Real* srcpnt ; /* source */ 646 Real* dstpnt ; /* destination */ 647 Integer srcOff ; /* source offset */ 648 Integer dstOff ; /* destination offset */ 649 Integer dim ; /* a simple dimension */ 650 651 if ( mDims * nDims == 0 ) { 652 return ; 653 } 654 655 srcOff = 0 ; 656 dstOff = 0 ; 657 658 srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; 659 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; 660 661 dim = nDims * usrHd->nNodes ; 662 663 fMtxBlkDaxpy( srcpnt, 664 dstpnt, 665 coef, 666 &mDims, 667 &dim ); 668 } 669 670 /*----------------------------------------------------------------------- 671 * 672 * lesBlkDyeax 673 * 674 * operation : 675 * 676 *----------------------------------------------------------------------- 677 */ 678 void lesBlkDyeax ( UsrHd usrHd, 679 Integer srcId, 680 Integer dstId, 681 Real* coef, 682 Integer mDims, 683 Integer nDims ) 684 { 685 char* funcName = "lesBlkDyeax" ; /* function name */ 686 Real* srcpnt ; /* source */ 687 Real* dstpnt ; /* destination */ 688 Integer srcOff ; /* source offset */ 689 Integer dstOff ; /* destination offset */ 690 Integer dim ; /* a simple dimension */ 691 692 if ( mDims * nDims == 0 ) { 693 lesZero ( usrHd, dstId, nDims ) ; 694 return ; 695 } 696 697 srcOff = 0 ; 698 dstOff = 0 ; 699 700 srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; 701 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; 702 703 dim = nDims * usrHd->nNodes ; 704 705 fMtxBlkDyeax( srcpnt, 706 dstpnt, 707 coef, 708 &mDims, 709 &dim ) ; 710 } 711 712 /*----------------------------------------------------------------------- 713 * 714 * lesBlkDmaxpy 715 * 716 * operation : 717 * 718 *----------------------------------------------------------------------- 719 */ 720 void lesBlkDmaxpy ( UsrHd usrHd, 721 Integer srcId, 722 Integer dstId, 723 Real* coef, 724 Integer mDims, 725 Integer nDims ) 726 { 727 char* funcName = "lesBlkDmaxpy" ; /* function name */ 728 Real* srcpnt ; /* source */ 729 Real* dstpnt ; /* destination */ 730 Integer srcOff ; /* source offset */ 731 Integer dstOff ; /* destination offset */ 732 Integer dim ; /* a simple dimension */ 733 734 if ( mDims * nDims == 0 ) { 735 return ; 736 } 737 738 srcOff = 0 ; 739 dstOff = 0 ; 740 741 srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; 742 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; 743 744 dim = nDims * usrHd->nNodes ; 745 746 fMtxBlkDmaxpy( srcpnt, 747 dstpnt, 748 coef, 749 &mDims, 750 &dim ); 751 } 752 753 /*----------------------------------------------------------------------- 754 * 755 * lesVdimCp 756 * 757 * operation : 758 * 759 *----------------------------------------------------------------------- 760 */ 761 void lesVdimCp ( UsrHd usrHd, 762 Integer srcId, 763 Integer dstId, 764 Integer nSrcDims, 765 Integer srcOff, 766 Integer nDstDims, 767 Integer dstOff, 768 Integer nDims ) 769 { 770 char* funcName = "lesVdimCp"; /* function name */ 771 Real* srcpnt ; /* source */ 772 Real* dstpnt ; /* destination */ 773 774 srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; 775 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; 776 777 fMtxVdimVecCp( srcpnt, 778 dstpnt, 779 &nSrcDims, 780 &nDstDims, 781 &nDims, 782 &(usrHd->nNodes) ); 783 } 784 785 /*----------------------------------------------------------------------- 786 * 787 * lesVdimDot2 788 * 789 * operation : 790 * 791 *----------------------------------------------------------------------- 792 */ 793 void lesVdimDot2 ( UsrHd usrHd, 794 Integer src1Id, 795 Integer src2Id, 796 Real* coef, 797 Integer nSrc1Dims, 798 Integer src1Off, 799 Integer nSrc2Dims, 800 Integer src2Off, 801 Integer nDims ) 802 { 803 char* funcName = "lesVdimDot2" ; /* function name */ 804 Real* src1pnt ; /* source 1 */ 805 Real* src2pnt ; /* source 2 */ 806 Real* coefp ; /* temporary coef on proc */ 807 808 809 if ( nDims == 0 ) { 810 return ; 811 } 812 813 coefp = (double *) malloc (nDims * sizeof(double)) ; 814 815 src1pnt = usrPointer ( usrHd, src1Id, src1Off, nSrc1Dims ) ; 816 src2pnt = usrPointer ( usrHd, src2Id, src2Off, nSrc2Dims ) ; 817 818 fMtxVdimVecDot2( src1pnt, 819 src2pnt, 820 coefp, 821 &nSrc1Dims, 822 &nSrc2Dims, 823 &nDims, 824 &(usrHd->nNodes) ); 825 826 drvAllreduce ( coefp, 827 coef, 828 &nDims ) ; 829 } 830 831 /*----------------------------------------------------------------------- 832 * 833 * lesVdimDaxpy 834 * 835 * operation : 836 * 837 *----------------------------------------------------------------------- 838 */ 839 void lesVdimDaxpy ( UsrHd usrHd, 840 Integer srcId, 841 Integer dstId, 842 Real* coef, 843 Integer nSrcDims, 844 Integer srcOff, 845 Integer nDstDims, 846 Integer dstOff, 847 Integer nDims ) 848 { 849 char* funcName = "lesVdimDaxpy" ; /* function name */ 850 Real* srcpnt ; /* source */ 851 Real* dstpnt ; /* destination */ 852 853 if ( nDims == 0 ) { 854 return ; 855 } 856 857 srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; 858 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; 859 860 fMtxVdimVecDaxpy( srcpnt, 861 dstpnt, 862 coef, 863 &nSrcDims, 864 &nDstDims, 865 &nDims, 866 &(usrHd->nNodes) ) ; 867 } 868 869 /*----------------------------------------------------------------------- 870 * 871 * lesApG 872 * 873 * operation : G(:,3*nenl,nenl) * Dp(:,nenl,1) = lesQ(;,nenl,3) 874 * 875 *----------------------------------------------------------------------- 876 */ 877 void lesApG ( UsrHd usrHd, 878 Integer srcId, 879 Integer dstId, 880 Integer nSrcDims, 881 Integer srcOff, 882 Integer nDstDims, 883 Integer dstOff ) 884 { 885 char* funcName = "lesApG" ; /* function name */ 886 Integer nDofs ; /* No. of dofs */ 887 Integer nPs ; /* No. of P dimension */ 888 Integer nQs ; /* No. of Q dimension */ 889 Integer pOff ; /* offset */ 890 Integer qOff ; /* offset */ 891 Real* srcpnt ; /* source */ 892 Real* dstpnt ; /* destin */ 893 894 nDofs = 4 ; 895 nPs = 1 ; 896 nQs = 3 ; 897 pOff = 3 * usrHd->nNodes ; 898 qOff = 0 * usrHd->nNodes ; 899 900 srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; 901 902 fMtxVdimVecMult( srcpnt, 903 usrHd->flowDiag+pOff, 904 usrHd->lesP, 905 &nSrcDims, 906 &nDofs, 907 &nPs, 908 &nPs, 909 &(usrHd->nNodes) ) ; 910 911 commOut ( usrHd->lesP, usrHd->ilwork, &nPs, 912 usrHd->iper, usrHd->iBC, usrHd->BC ); 913 914 fLesSparseApG( usrHd->colm, usrHd->rowp, usrHd->lhsP, 915 usrHd->lesP, usrHd->lesQ, &(usrHd->nNodes), 916 &(usrHd->nnz_tot)); 917 918 commIn ( usrHd->lesQ, usrHd->ilwork, &nQs, 919 usrHd->iper, usrHd->iBC, usrHd->BC ); 920 921 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; 922 923 fMtxVdimVecMult( usrHd->lesQ, 924 usrHd->flowDiag+qOff, 925 dstpnt, 926 &nQs, 927 &nDofs, 928 &nDstDims, 929 &nQs, 930 &(usrHd->nNodes) ) ; 931 } 932 933 /*----------------------------------------------------------------------- 934 * 935 * lesApKG 936 * 937 * operation : K(:,3*nenl,3*nenl) * Du(:,nenl,3) = lesQ(:,nenl,3) 938 * G(:,3*nenl, nenl) * Dp(:,nenl,1) = lesQ(:,nenl,3) 939 * 940 *----------------------------------------------------------------------- 941 */ 942 void lesApKG ( UsrHd usrHd, 943 Integer src1Id, 944 Integer src2Id, 945 Integer dstId, 946 Integer nSrc1Dims, 947 Integer src1Off, 948 Integer nSrc2Dims, 949 Integer src2Off, 950 Integer nDstDims, 951 Integer dstOff ) 952 { 953 char* funcName = "lesApKG" ; /* function name */ 954 Integer nDofs ; /* No. of Dofs */ 955 Integer p1Off ; /* Diag offset for P */ 956 Integer p2Off ; /* Diag offset for Q */ 957 Integer nP1s ; /* No. of P dimension */ 958 Integer nP2s ; /* No. of P dimension */ 959 Integer nPs ; /* No. of P dimension */ 960 Integer nQs ; /* No. of Q dimension */ 961 Integer qOff ; /* offset */ 962 Real* src1pnt ; /* Source 1 */ 963 Real* src2pnt ; /* Source 2 */ 964 Real* dstpnt ; /* destination */ 965 966 nDofs = 4 ; 967 nP1s = 3 ; 968 nP2s = 1 ; 969 nPs = nP1s + nP2s ; 970 nQs = 3 ; 971 p1Off = 0 * usrHd->nNodes ; 972 p2Off = 3 * usrHd->nNodes ; 973 qOff = 0 * usrHd->nNodes ; 974 975 src1pnt = usrPointer ( usrHd, src1Id, src1Off, nSrc1Dims ) ; 976 src2pnt = usrPointer ( usrHd, src2Id, src2Off, nSrc2Dims ) ; 977 978 fMtxVdimVecMult( src1pnt, 979 usrHd->flowDiag+p1Off, 980 usrHd->lesP+p1Off, 981 &nSrc1Dims, 982 &nDofs, 983 &nPs, 984 &nP1s, 985 &(usrHd->nNodes) ) ; 986 987 fMtxVdimVecMult( src2pnt, 988 usrHd->flowDiag+p2Off, 989 usrHd->lesP+p2Off, 990 &nSrc2Dims, 991 &nDofs, 992 &nPs, 993 &nP2s, 994 &(usrHd->nNodes) ); 995 996 commOut ( usrHd->lesP, usrHd->ilwork, &nPs, 997 usrHd->iper, usrHd->iBC, usrHd->BC ); 998 999 fLesSparseApKG( usrHd->colm, usrHd->rowp, usrHd->lhsK, 1000 usrHd->lhsP, usrHd->lesP, usrHd->lesQ, 1001 &(usrHd->nNodes),&(usrHd->nnz_tot)); 1002 1003 commIn ( usrHd->lesQ, usrHd->ilwork, &nQs, 1004 usrHd->iper, usrHd->iBC, usrHd->BC ); 1005 1006 1007 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; 1008 1009 fMtxVdimVecMult( usrHd->lesQ, 1010 usrHd->flowDiag+qOff, 1011 dstpnt, 1012 &nQs, 1013 &nDofs, 1014 &nDstDims, 1015 &nQs, 1016 &(usrHd->nNodes) ) ; 1017 } 1018 1019 /*----------------------------------------------------------------------- 1020 * 1021 * lesApNGt 1022 * 1023 * operation : -G^t(:,nenl,3*nenl) * Du(:,nenl,3) = lesQ(:,nenl,1) 1024 * 1025 *----------------------------------------------------------------------- 1026 */ 1027 void lesApNGt ( UsrHd usrHd, 1028 Integer srcId, 1029 Integer dstId, 1030 Integer nSrcDims, 1031 Integer srcOff, 1032 Integer nDstDims, 1033 Integer dstOff ) 1034 { 1035 char* funcName = "lesApNGt" ; /* function name */ 1036 Integer pOff ; /* Diag offset for P */ 1037 Integer qOff ; /* Diag offset for Q */ 1038 Integer nDofs ; /* No. of Dofs */ 1039 Integer nPs ; /* No. of P dimension */ 1040 Integer nQs ; /* No. of Q dimension */ 1041 Real* srcpnt ; /* Source */ 1042 Real* dstpnt ; /* Destination */ 1043 1044 nDofs = 4 ; 1045 nPs = 3 ; 1046 nQs = 1 ; 1047 pOff = 0 * usrHd->nNodes ; 1048 qOff = 3 * usrHd->nNodes ; 1049 1050 srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; 1051 1052 fMtxVdimVecMult( srcpnt, 1053 usrHd->flowDiag+pOff, 1054 usrHd->lesP, 1055 &nSrcDims, 1056 &nDofs, 1057 &nPs, 1058 &nPs, 1059 &(usrHd->nNodes) ) ; 1060 1061 commOut ( usrHd->lesP, usrHd->ilwork, &nPs, 1062 usrHd->iper, usrHd->iBC, usrHd->BC ); 1063 1064 fLesSparseApNGt( usrHd->colm, usrHd->rowp, 1065 usrHd->lhsP, usrHd->lesP, usrHd->lesQ, 1066 &(usrHd->nNodes),&(usrHd->nnz_tot)); 1067 1068 commIn ( usrHd->lesQ, usrHd->ilwork, &nQs, 1069 usrHd->iper, usrHd->iBC, usrHd->BC ); 1070 1071 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; 1072 1073 fMtxVdimVecMult( usrHd->lesQ, 1074 usrHd->flowDiag+qOff, 1075 dstpnt, 1076 &nQs, 1077 &nDofs, 1078 &nDstDims, 1079 &nQs, 1080 &(usrHd->nNodes) ) ; 1081 } 1082 1083 /*----------------------------------------------------------------------- 1084 * 1085 * lesApNGtC 1086 * 1087 * operation : -G^t(:,nenl,3*nenl) * Du(:,nenl,3) = lesQ(:,nenl,1) 1088 * C(:,nenl, nenl) * Dp(:,nenl,1) = lesQ(:,nenl,1) 1089 * 1090 *----------------------------------------------------------------------- 1091 */ 1092 void lesApNGtC ( UsrHd usrHd, 1093 Integer src1Id, 1094 Integer src2Id, 1095 Integer dstId, 1096 Integer nSrc1Dims, 1097 Integer src1Off, 1098 Integer nSrc2Dims, 1099 Integer src2Off, 1100 Integer nDstDims, 1101 Integer dstOff ) 1102 { 1103 char* funcName = "lesApNGtC" ; /* function name */ 1104 Integer p1Off ; /* Diag offset for P */ 1105 Integer p2Off ; /* Diag offset for P */ 1106 Integer qOff ; /* Diag offset for Q */ 1107 Integer nDofs ; /* No. of Dofs */ 1108 Integer nP1s ; /* No. of P dimension */ 1109 Integer nP2s ; /* No. of P dimension */ 1110 Integer nPs ; /* No. of P dimension */ 1111 Integer nQs ; /* No. of Q dimension */ 1112 Real* dstpnt ; /* Destination */ 1113 Real* src1pnt ; /* Source 1 */ 1114 Real* src2pnt ; /* Source 2 */ 1115 1116 nDofs = 4 ; 1117 nP1s = 3 ; 1118 nP2s = 1 ; 1119 nPs = nP1s + nP2s ; 1120 nQs = 1 ; 1121 p1Off = 0 * usrHd->nNodes ; 1122 p2Off = 3 * usrHd->nNodes ; 1123 qOff = 3 * usrHd->nNodes ; 1124 1125 src1pnt = usrPointer ( usrHd, src1Id, src1Off, nSrc1Dims ) ; 1126 src2pnt = usrPointer ( usrHd, src2Id, src2Off, nSrc2Dims ) ; 1127 1128 fMtxVdimVecMult( src1pnt, 1129 usrHd->flowDiag+p1Off, 1130 usrHd->lesP+p1Off, 1131 &nSrc1Dims, 1132 &nDofs, 1133 &nPs, 1134 &nP1s, 1135 &(usrHd->nNodes) ) ; 1136 1137 fMtxVdimVecMult( src2pnt, 1138 usrHd->flowDiag+p2Off, 1139 usrHd->lesP+p2Off, 1140 &nSrc2Dims, 1141 &nDofs, 1142 &nPs, 1143 &nP2s, 1144 &(usrHd->nNodes) ) ; 1145 commOut ( usrHd->lesP, usrHd->ilwork, &nPs, 1146 usrHd->iper, usrHd->iBC, usrHd->BC ); 1147 1148 fLesSparseApNGtC( usrHd->colm, usrHd->rowp, 1149 usrHd->lhsP, usrHd->lesP, usrHd->lesQ, 1150 &(usrHd->nNodes),&(usrHd->nnz_tot)); 1151 1152 commIn ( usrHd->lesQ, usrHd->ilwork, &nQs, 1153 usrHd->iper, usrHd->iBC, usrHd->BC ); 1154 1155 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; 1156 1157 fMtxVdimVecMult( usrHd->lesQ, 1158 usrHd->flowDiag+qOff, 1159 dstpnt, 1160 &nQs, 1161 &nDofs, 1162 &nDstDims, 1163 &nQs, 1164 &(usrHd->nNodes) ) ; 1165 } 1166 1167 /*----------------------------------------------------------------------- 1168 * 1169 * lesApFull 1170 * 1171 * operation : K * Du + G * Dp = lesQ(:,nenl,1:3) 1172 * -G^t * Du + C * Dp = lesQ(:,nenl,4:4) 1173 * 1174 *----------------------------------------------------------------------- 1175 */ 1176 void lesApFull ( UsrHd usrHd, 1177 Integer srcId, 1178 Integer dstId, 1179 Integer nSrcDims, 1180 Integer srcOff, 1181 Integer nDstDims, 1182 Integer dstOff ) 1183 { 1184 char* funcName = "lesApFull" ; /* function name */ 1185 Integer pOff ; /* Diag offset for P */ 1186 Integer qOff ; /* Diag offset for Q */ 1187 Integer nDofs ; /* No. of Dofs */ 1188 Integer nPs ; /* No. of P dimension */ 1189 Integer nQs ; /* No. of Q dimension */ 1190 Real* srcpnt ; /* Source */ 1191 Real* dstpnt ; /* Destination */ 1192 1193 nDofs = 4 ; 1194 nPs = 4 ; 1195 nQs = 4 ; 1196 pOff = 0 * usrHd->nNodes ; 1197 qOff = 0 * usrHd->nNodes ; 1198 1199 srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; 1200 1201 fMtxVdimVecMult( srcpnt, 1202 usrHd->flowDiag+pOff, 1203 usrHd->lesP, 1204 &nSrcDims, 1205 &nDofs, 1206 &nPs, 1207 &nPs, 1208 &(usrHd->nNodes) ) ; 1209 commOut ( usrHd->lesP, usrHd->ilwork, &nPs, 1210 usrHd->iper, usrHd->iBC, usrHd->BC ); 1211 1212 fLesSparseApFull( usrHd->colm, usrHd->rowp, usrHd->lhsK, 1213 usrHd->lhsP, usrHd->lesP, usrHd->lesQ, 1214 &(usrHd->nNodes),&(usrHd->nnz_tot)); 1215 1216 commIn ( usrHd->lesQ, usrHd->ilwork, &nQs, 1217 usrHd->iper, usrHd->iBC, usrHd->BC ); 1218 1219 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; 1220 1221 fMtxVdimVecMult( usrHd->lesQ, 1222 usrHd->flowDiag+qOff, 1223 dstpnt, 1224 &nQs, 1225 &nDofs, 1226 &nDstDims, 1227 &nQs, 1228 &(usrHd->nNodes) ) ; 1229 } 1230 1231 /*----------------------------------------------------------------------- 1232 * 1233 * lesApSclr 1234 * 1235 * operation : M(:,nenl,nenl) * Ds(:,nenl,1) = lesQ(:,nenl,1) 1236 * 1237 *----------------------------------------------------------------------- 1238 */ 1239 void lesApSclr ( UsrHd usrHd, 1240 Integer srcId, 1241 Integer dstId, 1242 Integer nSrcDims, 1243 Integer srcOff, 1244 Integer nDstDims, 1245 Integer dstOff ) 1246 { 1247 char* funcName = "lesApSclr" ; /* function name */ 1248 Integer pOff ; /* Diag offset for P */ 1249 Integer qOff ; /* Diag offset for Q */ 1250 Integer nDofs ; /* No. of Dofs */ 1251 Integer nPs ; /* No. of P dimension */ 1252 Integer nQs ; /* No. of Q dimension */ 1253 Real* srcpnt ; /* Source */ 1254 Real* dstpnt ; /* Destination */ 1255 Integer lhsStiffFlag ; 1256 double sclrRegFct ; 1257 1258 nDofs = 1 ; 1259 nPs = 1 ; 1260 nQs = 1 ; 1261 pOff = 0 ; 1262 qOff = 0 ; 1263 1264 lhsStiffFlag = 0 ; 1265 sclrRegFct = 0 ; 1266 1267 srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; 1268 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; 1269 1270 1271 fMtxVdimVecMult ( srcpnt, 1272 usrHd->sclrDiag+pOff, 1273 usrHd->lesP, 1274 &nSrcDims, 1275 &nDofs, 1276 &nPs, 1277 &nPs, 1278 &(usrHd->nNodes) ) ; 1279 commOut ( usrHd->lesP, usrHd->ilwork, &nPs, 1280 usrHd->iper, usrHd->iBC, usrHd->BC ); 1281 1282 fLesSparseApSclr( usrHd->colm, usrHd->rowp, usrHd->lhsS, 1283 usrHd->lesP, usrHd->lesQ, 1284 &(usrHd->nNodes),&(usrHd->nnz_tot)); 1285 1286 commIn ( usrHd->lesQ, usrHd->ilwork, &nQs, 1287 usrHd->iper, usrHd->iBC, usrHd->BC ); 1288 1289 1290 if ( lhsStiffFlag && sclrRegFct != 0 ) { 1291 1292 fMtxVdimVecMult ( usrHd->lesQ, 1293 usrHd->sclrDiag+qOff, 1294 usrHd->lesP, 1295 &nQs, 1296 &nDofs, 1297 &nDstDims, 1298 &nQs, 1299 &(usrHd->nNodes) ) ; 1300 1301 flesDaxpy ( srcpnt, 1302 usrHd->lesP, 1303 &sclrRegFct, 1304 &nDstDims, 1305 &(usrHd->nNodes) ) ; 1306 1307 flesCp ( usrHd->lesP, 1308 dstpnt, 1309 &nDstDims, 1310 &(usrHd->nNodes) ) ; 1311 1312 } else { 1313 1314 fMtxVdimVecMult ( usrHd->lesQ, 1315 usrHd->sclrDiag+qOff, 1316 dstpnt, 1317 &nQs, 1318 &nDofs, 1319 &nDstDims, 1320 &nQs, 1321 &(usrHd->nNodes) ) ; 1322 } 1323 1324 } 1325 1326 /********************************************************************* 1327 * lesPrecPPE 1328 * outer routine to solve PPE 1329 *******************************************************************/ 1330 1331 void lesPrecPPE(UsrHd usrHd, 1332 Integer srcId, 1333 Integer dstId, 1334 Integer nSrcDims, 1335 Integer srcOff, 1336 Integer nDstDims, 1337 Integer dstOff) 1338 { 1339 Real* srcpnt ; /* Source R */ 1340 Real* dstpnt ; /* Destination Z */ 1341 srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; 1342 dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; 1343 #ifdef AMG 1344 ramg_interface( usrHd->colm, 1345 usrHd->rowp,usrHd->lhsK,usrHd->lhsP,usrHd->flowDiag, 1346 srcpnt,dstpnt, 1347 usrHd->ilwork,usrHd->BC,usrHd->iBC,usrHd->iper 1348 ); 1349 #endif 1350 return; 1351 } 1352 1353