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 */
lesPrepDiag(UsrHd usrHd)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 */
lesDiagScaleCp(UsrHd usrHd,Integer srcId,Integer dstId,Integer nSrcDims,Integer srcOff,Integer nDstDims,Integer dstOff,Integer diagOff,Integer nDims)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 */
lesZero(UsrHd usrHd,Integer dstId,Integer nDims)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 */
lesCp(UsrHd usrHd,Integer srcId,Integer dstId,Integer nDims)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 */
lesScale(UsrHd usrHd,Integer dstId,Real coef,Integer nDims)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 */
lesScaleCp(UsrHd usrHd,Integer srcId,Integer dstId,Real coef,Integer nDims)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 */
lesAdd(UsrHd usrHd,Integer srcId,Integer dstId,Integer nDims)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 */
lesSub(UsrHd usrHd,Integer srcId,Integer dstId,Integer nDims)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 */
lesDot1(UsrHd usrHd,Integer srcId,Integer nDims)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 */
lesDot2(UsrHd usrHd,Integer src1Id,Integer src2Id,Integer nDims)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 */
lesDaxpy(UsrHd usrHd,Integer srcId,Integer dstId,Real coef,Integer nDims)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 */
lesDxpay(UsrHd usrHd,Integer srcId,Integer dstId,Real coef,Integer nDims)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 */
lesInv(UsrHd usrHd,Integer dstId,Integer nDims)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 */
lesBlkDot2(UsrHd usrHd,Integer src1Id,Integer src2Id,Real * values,Integer mDims,Integer nDims)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 */
lesBlkDaxpy(UsrHd usrHd,Integer srcId,Integer dstId,Real * coef,Integer mDims,Integer nDims)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 */
lesBlkDyeax(UsrHd usrHd,Integer srcId,Integer dstId,Real * coef,Integer mDims,Integer nDims)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 */
lesBlkDmaxpy(UsrHd usrHd,Integer srcId,Integer dstId,Real * coef,Integer mDims,Integer nDims)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 */
lesVdimCp(UsrHd usrHd,Integer srcId,Integer dstId,Integer nSrcDims,Integer srcOff,Integer nDstDims,Integer dstOff,Integer nDims)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 */
lesVdimDot2(UsrHd usrHd,Integer src1Id,Integer src2Id,Real * coef,Integer nSrc1Dims,Integer src1Off,Integer nSrc2Dims,Integer src2Off,Integer nDims)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 */
lesVdimDaxpy(UsrHd usrHd,Integer srcId,Integer dstId,Real * coef,Integer nSrcDims,Integer srcOff,Integer nDstDims,Integer dstOff,Integer nDims)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 */
lesApG(UsrHd usrHd,Integer srcId,Integer dstId,Integer nSrcDims,Integer srcOff,Integer nDstDims,Integer dstOff)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 */
lesApKG(UsrHd usrHd,Integer src1Id,Integer src2Id,Integer dstId,Integer nSrc1Dims,Integer src1Off,Integer nSrc2Dims,Integer src2Off,Integer nDstDims,Integer dstOff)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 */
lesApNGt(UsrHd usrHd,Integer srcId,Integer dstId,Integer nSrcDims,Integer srcOff,Integer nDstDims,Integer dstOff)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 */
lesApNGtC(UsrHd usrHd,Integer src1Id,Integer src2Id,Integer dstId,Integer nSrc1Dims,Integer src1Off,Integer nSrc2Dims,Integer src2Off,Integer nDstDims,Integer dstOff)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 */
lesApFull(UsrHd usrHd,Integer srcId,Integer dstId,Integer nSrcDims,Integer srcOff,Integer nDstDims,Integer dstOff)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 */
lesApSclr(UsrHd usrHd,Integer srcId,Integer dstId,Integer nSrcDims,Integer srcOff,Integer nDstDims,Integer dstOff)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
lesPrecPPE(UsrHd usrHd,Integer srcId,Integer dstId,Integer nSrcDims,Integer srcOff,Integer nDstDims,Integer dstOff)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