/*====================================================================== * * lestools.c : Linear Algebra Solver Tools * * small single character : vector or matrix * *====================================================================== */ #include "les.h" #include "usr.h" #ifdef intel void DRVSCLRDIAG( double *sclrDiag, int *ilwork, int *iBC, double *BC, int *iper, int *rowp, int *colm, double *lhsS); void FMTXBLKDAXPY(double *srcpnt, double *dstpnt, double *coef, int *mDims, int *dim); void FMTXBLKDYEAX(double *srcpnt, double *dstpnt, double *coef, int *mDims, int *dim); void FMTXBLKDMAXPY(double *srcpnt, double *dstpnt, double *coef, int *mDims, int *dim); void FMTXVDIMVECCP(double *srcpnt, double *dstpnt, int *nSrcDims, int *nDstDims, int *nDims, int *nNodes ); void DRVLESPREPDIAG( double *flowDiag, int *ilwork, int *iBC, double *BC, int *iper, int *rowp, int *colm, double *lhsK, double *lhsP) ; void FMTXVDIMVECMULT( double* , double*, double *dstpnt, int *nSrcDims, int *nDofs, int *nDstDims, int *nDims, int *nNodes); void FLESZERO( double *dstpnt, int *nDims, int *nNodes); void FLESCP( double *srcpnt, double *dstpnt, int *nDims, int *dim ) ; void FLESSCALE( double *dstpnt, double *coef, int *nDims, int *dim ) ; void FLESSCALECP( double *srcpnt, double *dstpnt, double *coef, int *nDims, int *dim ) ; void FLESADD( double *srcpnt, double *dstpnt, int *nDims, int *dim ) ; void FLESSUB( double *srcpnt, double *dstpnt, int *nDims, int *dim ) ; void DRVALLREDUCESCLR(double *tmpp, double *tmp); void DRVALLREDUCE(double *valuesp, double *values, int* mDims); double FLESDOT1( double *srcpnt, int* nDims, int* dim ); double FLESDOT2(double *src1pnt, double *src2pnt, int* nDims, int* dim ); void FLESDAXPY( double *srcpnt, double *lesP, double *sclrRegFct, int *nDstDims, int *nNodes); void FLESDXPAY( double *srcpnt, double *dstpnt, double *coef, int *nDims, int *dim); void FLESINV( double *dstpnt, int *nDims, int *dim ) ; void FMTXBLKDOT2 (double *src1pnt, double *src2pnt, double *valuesp, int* mDims, int* dim); void COMMIN( double *lesQ, int *ilwork, int *nPs, int *iper, int *iBC, double *BC); void COMMOUT( double *lesP, int *ilwork, int *nQs, int *iper, int *iBC, double *BC); void FLESSPARSEAPFULL(int *colm, int *rowp, double *lhsK, double *lhsP, double *lesP, double *lesQ, int *nNodes, int *nnz); void FLESSPARSEAPSCLR(int *colm, int *rowp, double *lhsS, double *lesP, double *lesQ, int *nNodes, int *nnz); void FMTXVDIMVECDOT2 (double *src1pnt, double *src2pnt, double *coefp, int *nSrc1Dims,int *nSrc2Dims, int *nDims, int *nNodes); void FMTXVDIMVECDAXPY ( double *srcpnt, double *dstpnt, double *coef, int *nSrcDims, int *nDstDims, int *nDims, int *nNodes); void FLESSPARSEAPG ( int *colm, int *rowp, double *lhsP, double *lesP, double *lesQ, int *nNodes, int *nnz); void FLESSPARSEAPNGT ( int *colm, int *rowp, double *lhsP, double *lesP, double *lesQ, int *nNodes, int *nnz); void FLESSPARSEAPNGTC ( int *colm, int *rowp, double *lhsP, double *lesP, double *lesQ, int *nNodes, int *nnz); void FLESSPARSEAPKG ( int *colm, int *rowp, double *lhsK, double *lhsP, double *lesP, double *lesQ, int *nNodes, int *nnz); void RAMG_INTERFACE ( int *colm, int *rowp, double *lhsK,double *lhsP, double *mcgR,double *mcgZ, int *ilwork, double *BC, int *iBC,int *iper); #endif /*---------------------------------------------------------------------- * * lesPreDiag * * operation : a = 1/sqrt(abs(a)) * *---------------------------------------------------------------------- */ void lesPrepDiag( UsrHd usrHd ) { char* funcName = "lesPrepDiag" ; /* function name */ if ( (usrHd->eqnType) == 1 ) { /* continuity & momentum */ drvlesPrepDiag( usrHd->flowDiag, usrHd->ilwork, usrHd->iBC, usrHd->BC, usrHd->iper, usrHd->rowp, usrHd->colm, usrHd->lhsK, usrHd->lhsP) ; } if ( (usrHd->eqnType) == 2 ) { /* temperature or scalar variable */ drvsclrDiag ( usrHd->sclrDiag, usrHd->ilwork, usrHd->iBC, usrHd->BC, usrHd->iper, usrHd->rowp, usrHd->colm, usrHd->lhsS ) ; } } /*---------------------------------------------------------------------- * * lesDiagScaleCp * * operation : c = a * b * *---------------------------------------------------------------------- */ void lesDiagScaleCp ( UsrHd usrHd, Integer srcId, Integer dstId, Integer nSrcDims, Integer srcOff, Integer nDstDims, Integer dstOff, Integer diagOff, Integer nDims ) { char* funcName = "lesDiagScaleCp" ; /* function name */ Integer nDofs ; /* No. of Dofs */ Real* dstpnt ; /* destination */ Real* srcpnt ; /* source */ if ( (usrHd->eqnType) == 1 ) { nDofs = 4 ; srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; fMtxVdimVecMult( srcpnt, usrHd->flowDiag+diagOff*usrHd->nNodes, dstpnt, &nSrcDims, &nDofs, &nDstDims, &nDims, &(usrHd->nNodes) ) ; } if ( (usrHd->eqnType) == 2 ) { nDofs = 1 ; srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; fMtxVdimVecMult( srcpnt, usrHd->sclrDiag+diagOff*usrHd->nNodes, dstpnt, &nSrcDims, &nDofs, &nDstDims, &nDims, &(usrHd->nNodes) ) ; } } /*---------------------------------------------------------------------- * * lesZero * * operation : a = 0 * *---------------------------------------------------------------------- */ void lesZero ( UsrHd usrHd, Integer dstId, Integer nDims ) { char* funcName = "lesZero" ; /* function namea */ Real* dstpnt ; /* destination */ Integer dstOff ; /* destination offset */ dstOff = 0 ; dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ); flesZero( dstpnt, &nDims, &(usrHd->nNodes) ); } /*---------------------------------------------------------------------- * * lesCp * * operation : b = a * *----------------------------------------------------------------------- */ void lesCp ( UsrHd usrHd, Integer srcId, Integer dstId, Integer nDims ) { char* funcName = "lesCp" ; /* function name */ Real* srcpnt ; /* source */ Real* dstpnt ; /* destination */ Integer dim ; /* a simple dimension */ Integer srcOff ; /* source offset */ Integer dstOff ; /* destination offset */ srcOff = 0 ; dstOff = 0 ; srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; dim = usrHd->nNodes ; flesCp( srcpnt, dstpnt, &nDims, &dim ) ; } /*----------------------------------------------------------------------- * * lesScale * * operation : a = a * coef * *----------------------------------------------------------------------- */ void lesScale ( UsrHd usrHd, Integer dstId, Real coef, Integer nDims ) { char* funcName = "lesScale" ; /* function name */ Real* dstpnt ; /* destination */ Integer dstOff ; /* destination offset */ Integer dim ; /* a simple dimension */ dstOff = 0 ; dim = usrHd->nNodes ; dstpnt = usrPointer( usrHd, dstId, dstOff, nDims ) ; flesScale( dstpnt, &coef, &nDims, &dim ) ; } /*----------------------------------------------------------------------- * * lesScaleCp * * operation : b = a * coef * *----------------------------------------------------------------------- */ void lesScaleCp ( UsrHd usrHd, Integer srcId, Integer dstId, Real coef, Integer nDims ) { char* funcName = "lesScaleCp" ; /* function name */ Real* srcpnt ; /* source */ Real* dstpnt ; /* destination */ Integer dim ; /* a simple dimension */ Integer srcOff ; /* source offset */ Integer dstOff ; /* destination offset */ srcOff = 0 ; dstOff = 0 ; srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; dim = usrHd->nNodes ; flesScaleCp( srcpnt, dstpnt, &coef, &nDims, &dim ) ; } /*----------------------------------------------------------------------- * * lesAdd * * operation : b = b + a * *----------------------------------------------------------------------- */ void lesAdd ( UsrHd usrHd, Integer srcId, Integer dstId, Integer nDims ) { char* funcName = "lesAdd" ; /* function name */ Real* srcpnt ; /* source */ Real* dstpnt ; /* destination */ Integer srcOff ; /* source offset */ Integer dstOff ; /* destination offset */ Integer dim ; /* a simple dimension */ srcOff = 0 ; dstOff = 0 ; srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; dim = usrHd->nNodes ; flesAdd( srcpnt, dstpnt, &nDims, &dim ) ; } /*----------------------------------------------------------------------- * * lesSub * * operation : b = b - a * *----------------------------------------------------------------------- */ void lesSub ( UsrHd usrHd, Integer srcId, Integer dstId, Integer nDims ) { char* funcName = "lesSub" ; /* function name */ Real* srcpnt ; /* source */ Real* dstpnt ; /* destination */ Integer srcOff ; /* source offset */ Integer dstOff ; /* destination offset */ Integer dim ; /* a simple dimension */ srcOff = 0 ; dstOff = 0 ; srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; dim = usrHd->nNodes ; flesSub( srcpnt, dstpnt, &nDims, &dim ) ; } /*----------------------------------------------------------------------- * * lesDot1 * * operation : tmp = tmp + a * a * *----------------------------------------------------------------------- */ Real lesDot1 ( UsrHd usrHd, Integer srcId, Integer nDims ) { char* funcName = "lesDot1" ; /* function name */ Real* srcpnt ; /* source */ Integer srcOff ; /* source offset */ Integer dim ; /* a simple dimension */ Real tmp ; /* a temporary value */ Real tmpp ; /* a temporary value on each proc. */ srcOff = 0 ; srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; dim = usrHd->nNodes ; tmpp = flesDot1( srcpnt, &nDims, &dim ) ; drvAllreducesclr ( &tmpp, &tmp ) ; #ifdef AMG ramg_normcheck(&tmp); #endif return ( tmp ) ; } /*----------------------------------------------------------------------- * * lesDot2 * * operation : tmp = tmp + a * b * *----------------------------------------------------------------------- */ Real lesDot2 ( UsrHd usrHd, Integer src1Id, Integer src2Id, Integer nDims ) { char* funcName = "lesDot2" ; /* function name */ Real* src1pnt ; /* source 1 */ Real* src2pnt ; /* source 2 */ Integer src1Off ; /* source 1 offset */ Integer src2Off ; /* source 2 offset */ Integer dim ; /* a simple dimension */ Real tmp ; /* a temporary value */ Real tmpp ; /* a temporary value on each proc. */ src1Off = 0 ; src2Off = 0 ; src1pnt = usrPointer ( usrHd, src1Id, src1Off, nDims ); src2pnt = usrPointer ( usrHd, src2Id, src2Off, nDims ); dim = usrHd->nNodes ; tmpp = flesDot2( src1pnt, src2pnt, &nDims, &dim ) ; drvAllreducesclr ( &tmpp, &tmp ) ; return ( tmp ) ; } /*----------------------------------------------------------------------- * * lesDaxpy * * operation : y = y + coef * x * *----------------------------------------------------------------------- */ void lesDaxpy ( UsrHd usrHd, Integer srcId, Integer dstId, Real coef, Integer nDims ) { char* funcName = "lesDapxy" ; /* function name */ Real* srcpnt ; /* source */ Real* dstpnt ; /* destination */ Integer srcOff ; /* source offset */ Integer dstOff ; /* destination offset */ Integer dim ; /* a simple dimension */ srcOff = 0 ; dstOff = 0 ; srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; dim = usrHd->nNodes ; flesDaxpy( srcpnt, dstpnt, &coef, &nDims, &dim ) ; } /*----------------------------------------------------------------------- * * lesDxpay * * operation : y = coef * y + x * *----------------------------------------------------------------------- */ void lesDxpay ( UsrHd usrHd, Integer srcId, Integer dstId, Real coef, Integer nDims ) { char* funcName = "lesDxpay" ; /* function name */ Real* srcpnt ; /* source */ Real* dstpnt ; /* destination */ Integer srcOff ; /* source offset */ Integer dstOff ; /* destination offset */ Integer dim ; /* a simple dimension */ srcOff = 0 ; dstOff = 0 ; srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; dim = usrHd->nNodes ; flesDxpay( srcpnt, dstpnt, &coef, &nDims, &dim ) ; } /*----------------------------------------------------------------------- * * lesInv * * operation : x = 1. / x * *----------------------------------------------------------------------- */ void lesInv ( UsrHd usrHd, Integer dstId, Integer nDims ) { char* funcName = "lesInv" ; /* function name */ Integer dim ; /* a simple dimension */ Real* dstpnt ; /* destination */ Integer dstOff ; /* destination offset */ dstOff = 0 ; dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; dim = usrHd->nNodes ; flesInv( dstpnt, &nDims, &dim ) ; } /*----------------------------------------------------------------------- * * lesBlkDot2 * * operation : * *----------------------------------------------------------------------- */ void lesBlkDot2 ( UsrHd usrHd, Integer src1Id, Integer src2Id, Real* values, Integer mDims, Integer nDims ) { char* funcName = "lesBlkDot2" ; /* function name */ Real* src1pnt ; /* source 1 */ Real* src2pnt ; /* source 2 */ Integer src1Off ; /* source 1 offset */ Integer src2Off ; /* source 2 offset */ Integer dim ; /* a simple dimension */ Real* valuesp ; /* temporary values on proc */ if ( mDims * nDims == 0 ) { return ; } valuesp = (double *) malloc (mDims * sizeof(double)) ; src1Off = 0 ; src2Off = 0 ; src1pnt = usrPointer ( usrHd, src1Id, src1Off, nDims ) ; src2pnt = usrPointer ( usrHd, src2Id, src2Off, nDims ) ; dim = nDims * usrHd->nNodes ; fMtxBlkDot2( src1pnt, src2pnt, valuesp, &mDims, &dim ) ; drvAllreduce ( valuesp, values, &mDims ) ; free(valuesp); } /*----------------------------------------------------------------------- * * lesBlkDaxpy * * operation : * *----------------------------------------------------------------------- */ void lesBlkDaxpy ( UsrHd usrHd, Integer srcId, Integer dstId, Real* coef, Integer mDims, Integer nDims ) { char* funcName = "lesBlkDaxpy" ; /* function name */ Real* srcpnt ; /* source */ Real* dstpnt ; /* destination */ Integer srcOff ; /* source offset */ Integer dstOff ; /* destination offset */ Integer dim ; /* a simple dimension */ if ( mDims * nDims == 0 ) { return ; } srcOff = 0 ; dstOff = 0 ; srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; dim = nDims * usrHd->nNodes ; fMtxBlkDaxpy( srcpnt, dstpnt, coef, &mDims, &dim ); } /*----------------------------------------------------------------------- * * lesBlkDyeax * * operation : * *----------------------------------------------------------------------- */ void lesBlkDyeax ( UsrHd usrHd, Integer srcId, Integer dstId, Real* coef, Integer mDims, Integer nDims ) { char* funcName = "lesBlkDyeax" ; /* function name */ Real* srcpnt ; /* source */ Real* dstpnt ; /* destination */ Integer srcOff ; /* source offset */ Integer dstOff ; /* destination offset */ Integer dim ; /* a simple dimension */ if ( mDims * nDims == 0 ) { lesZero ( usrHd, dstId, nDims ) ; return ; } srcOff = 0 ; dstOff = 0 ; srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; dim = nDims * usrHd->nNodes ; fMtxBlkDyeax( srcpnt, dstpnt, coef, &mDims, &dim ) ; } /*----------------------------------------------------------------------- * * lesBlkDmaxpy * * operation : * *----------------------------------------------------------------------- */ void lesBlkDmaxpy ( UsrHd usrHd, Integer srcId, Integer dstId, Real* coef, Integer mDims, Integer nDims ) { char* funcName = "lesBlkDmaxpy" ; /* function name */ Real* srcpnt ; /* source */ Real* dstpnt ; /* destination */ Integer srcOff ; /* source offset */ Integer dstOff ; /* destination offset */ Integer dim ; /* a simple dimension */ if ( mDims * nDims == 0 ) { return ; } srcOff = 0 ; dstOff = 0 ; srcpnt = usrPointer ( usrHd, srcId, srcOff, nDims ) ; dstpnt = usrPointer ( usrHd, dstId, dstOff, nDims ) ; dim = nDims * usrHd->nNodes ; fMtxBlkDmaxpy( srcpnt, dstpnt, coef, &mDims, &dim ); } /*----------------------------------------------------------------------- * * lesVdimCp * * operation : * *----------------------------------------------------------------------- */ void lesVdimCp ( UsrHd usrHd, Integer srcId, Integer dstId, Integer nSrcDims, Integer srcOff, Integer nDstDims, Integer dstOff, Integer nDims ) { char* funcName = "lesVdimCp"; /* function name */ Real* srcpnt ; /* source */ Real* dstpnt ; /* destination */ srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; fMtxVdimVecCp( srcpnt, dstpnt, &nSrcDims, &nDstDims, &nDims, &(usrHd->nNodes) ); } /*----------------------------------------------------------------------- * * lesVdimDot2 * * operation : * *----------------------------------------------------------------------- */ void lesVdimDot2 ( UsrHd usrHd, Integer src1Id, Integer src2Id, Real* coef, Integer nSrc1Dims, Integer src1Off, Integer nSrc2Dims, Integer src2Off, Integer nDims ) { char* funcName = "lesVdimDot2" ; /* function name */ Real* src1pnt ; /* source 1 */ Real* src2pnt ; /* source 2 */ Real* coefp ; /* temporary coef on proc */ if ( nDims == 0 ) { return ; } coefp = (double *) malloc (nDims * sizeof(double)) ; src1pnt = usrPointer ( usrHd, src1Id, src1Off, nSrc1Dims ) ; src2pnt = usrPointer ( usrHd, src2Id, src2Off, nSrc2Dims ) ; fMtxVdimVecDot2( src1pnt, src2pnt, coefp, &nSrc1Dims, &nSrc2Dims, &nDims, &(usrHd->nNodes) ); drvAllreduce ( coefp, coef, &nDims ) ; } /*----------------------------------------------------------------------- * * lesVdimDaxpy * * operation : * *----------------------------------------------------------------------- */ void lesVdimDaxpy ( UsrHd usrHd, Integer srcId, Integer dstId, Real* coef, Integer nSrcDims, Integer srcOff, Integer nDstDims, Integer dstOff, Integer nDims ) { char* funcName = "lesVdimDaxpy" ; /* function name */ Real* srcpnt ; /* source */ Real* dstpnt ; /* destination */ if ( nDims == 0 ) { return ; } srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; fMtxVdimVecDaxpy( srcpnt, dstpnt, coef, &nSrcDims, &nDstDims, &nDims, &(usrHd->nNodes) ) ; } /*----------------------------------------------------------------------- * * lesApG * * operation : G(:,3*nenl,nenl) * Dp(:,nenl,1) = lesQ(;,nenl,3) * *----------------------------------------------------------------------- */ void lesApG ( UsrHd usrHd, Integer srcId, Integer dstId, Integer nSrcDims, Integer srcOff, Integer nDstDims, Integer dstOff ) { char* funcName = "lesApG" ; /* function name */ Integer nDofs ; /* No. of dofs */ Integer nPs ; /* No. of P dimension */ Integer nQs ; /* No. of Q dimension */ Integer pOff ; /* offset */ Integer qOff ; /* offset */ Real* srcpnt ; /* source */ Real* dstpnt ; /* destin */ nDofs = 4 ; nPs = 1 ; nQs = 3 ; pOff = 3 * usrHd->nNodes ; qOff = 0 * usrHd->nNodes ; srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; fMtxVdimVecMult( srcpnt, usrHd->flowDiag+pOff, usrHd->lesP, &nSrcDims, &nDofs, &nPs, &nPs, &(usrHd->nNodes) ) ; commOut ( usrHd->lesP, usrHd->ilwork, &nPs, usrHd->iper, usrHd->iBC, usrHd->BC ); fLesSparseApG( usrHd->colm, usrHd->rowp, usrHd->lhsP, usrHd->lesP, usrHd->lesQ, &(usrHd->nNodes), &(usrHd->nnz_tot)); commIn ( usrHd->lesQ, usrHd->ilwork, &nQs, usrHd->iper, usrHd->iBC, usrHd->BC ); dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; fMtxVdimVecMult( usrHd->lesQ, usrHd->flowDiag+qOff, dstpnt, &nQs, &nDofs, &nDstDims, &nQs, &(usrHd->nNodes) ) ; } /*----------------------------------------------------------------------- * * lesApKG * * operation : K(:,3*nenl,3*nenl) * Du(:,nenl,3) = lesQ(:,nenl,3) * G(:,3*nenl, nenl) * Dp(:,nenl,1) = lesQ(:,nenl,3) * *----------------------------------------------------------------------- */ void lesApKG ( UsrHd usrHd, Integer src1Id, Integer src2Id, Integer dstId, Integer nSrc1Dims, Integer src1Off, Integer nSrc2Dims, Integer src2Off, Integer nDstDims, Integer dstOff ) { char* funcName = "lesApKG" ; /* function name */ Integer nDofs ; /* No. of Dofs */ Integer p1Off ; /* Diag offset for P */ Integer p2Off ; /* Diag offset for Q */ Integer nP1s ; /* No. of P dimension */ Integer nP2s ; /* No. of P dimension */ Integer nPs ; /* No. of P dimension */ Integer nQs ; /* No. of Q dimension */ Integer qOff ; /* offset */ Real* src1pnt ; /* Source 1 */ Real* src2pnt ; /* Source 2 */ Real* dstpnt ; /* destination */ nDofs = 4 ; nP1s = 3 ; nP2s = 1 ; nPs = nP1s + nP2s ; nQs = 3 ; p1Off = 0 * usrHd->nNodes ; p2Off = 3 * usrHd->nNodes ; qOff = 0 * usrHd->nNodes ; src1pnt = usrPointer ( usrHd, src1Id, src1Off, nSrc1Dims ) ; src2pnt = usrPointer ( usrHd, src2Id, src2Off, nSrc2Dims ) ; fMtxVdimVecMult( src1pnt, usrHd->flowDiag+p1Off, usrHd->lesP+p1Off, &nSrc1Dims, &nDofs, &nPs, &nP1s, &(usrHd->nNodes) ) ; fMtxVdimVecMult( src2pnt, usrHd->flowDiag+p2Off, usrHd->lesP+p2Off, &nSrc2Dims, &nDofs, &nPs, &nP2s, &(usrHd->nNodes) ); commOut ( usrHd->lesP, usrHd->ilwork, &nPs, usrHd->iper, usrHd->iBC, usrHd->BC ); fLesSparseApKG( usrHd->colm, usrHd->rowp, usrHd->lhsK, usrHd->lhsP, usrHd->lesP, usrHd->lesQ, &(usrHd->nNodes),&(usrHd->nnz_tot)); commIn ( usrHd->lesQ, usrHd->ilwork, &nQs, usrHd->iper, usrHd->iBC, usrHd->BC ); dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; fMtxVdimVecMult( usrHd->lesQ, usrHd->flowDiag+qOff, dstpnt, &nQs, &nDofs, &nDstDims, &nQs, &(usrHd->nNodes) ) ; } /*----------------------------------------------------------------------- * * lesApNGt * * operation : -G^t(:,nenl,3*nenl) * Du(:,nenl,3) = lesQ(:,nenl,1) * *----------------------------------------------------------------------- */ void lesApNGt ( UsrHd usrHd, Integer srcId, Integer dstId, Integer nSrcDims, Integer srcOff, Integer nDstDims, Integer dstOff ) { char* funcName = "lesApNGt" ; /* function name */ Integer pOff ; /* Diag offset for P */ Integer qOff ; /* Diag offset for Q */ Integer nDofs ; /* No. of Dofs */ Integer nPs ; /* No. of P dimension */ Integer nQs ; /* No. of Q dimension */ Real* srcpnt ; /* Source */ Real* dstpnt ; /* Destination */ nDofs = 4 ; nPs = 3 ; nQs = 1 ; pOff = 0 * usrHd->nNodes ; qOff = 3 * usrHd->nNodes ; srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; fMtxVdimVecMult( srcpnt, usrHd->flowDiag+pOff, usrHd->lesP, &nSrcDims, &nDofs, &nPs, &nPs, &(usrHd->nNodes) ) ; commOut ( usrHd->lesP, usrHd->ilwork, &nPs, usrHd->iper, usrHd->iBC, usrHd->BC ); fLesSparseApNGt( usrHd->colm, usrHd->rowp, usrHd->lhsP, usrHd->lesP, usrHd->lesQ, &(usrHd->nNodes),&(usrHd->nnz_tot)); commIn ( usrHd->lesQ, usrHd->ilwork, &nQs, usrHd->iper, usrHd->iBC, usrHd->BC ); dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; fMtxVdimVecMult( usrHd->lesQ, usrHd->flowDiag+qOff, dstpnt, &nQs, &nDofs, &nDstDims, &nQs, &(usrHd->nNodes) ) ; } /*----------------------------------------------------------------------- * * lesApNGtC * * operation : -G^t(:,nenl,3*nenl) * Du(:,nenl,3) = lesQ(:,nenl,1) * C(:,nenl, nenl) * Dp(:,nenl,1) = lesQ(:,nenl,1) * *----------------------------------------------------------------------- */ void lesApNGtC ( UsrHd usrHd, Integer src1Id, Integer src2Id, Integer dstId, Integer nSrc1Dims, Integer src1Off, Integer nSrc2Dims, Integer src2Off, Integer nDstDims, Integer dstOff ) { char* funcName = "lesApNGtC" ; /* function name */ Integer p1Off ; /* Diag offset for P */ Integer p2Off ; /* Diag offset for P */ Integer qOff ; /* Diag offset for Q */ Integer nDofs ; /* No. of Dofs */ Integer nP1s ; /* No. of P dimension */ Integer nP2s ; /* No. of P dimension */ Integer nPs ; /* No. of P dimension */ Integer nQs ; /* No. of Q dimension */ Real* dstpnt ; /* Destination */ Real* src1pnt ; /* Source 1 */ Real* src2pnt ; /* Source 2 */ nDofs = 4 ; nP1s = 3 ; nP2s = 1 ; nPs = nP1s + nP2s ; nQs = 1 ; p1Off = 0 * usrHd->nNodes ; p2Off = 3 * usrHd->nNodes ; qOff = 3 * usrHd->nNodes ; src1pnt = usrPointer ( usrHd, src1Id, src1Off, nSrc1Dims ) ; src2pnt = usrPointer ( usrHd, src2Id, src2Off, nSrc2Dims ) ; fMtxVdimVecMult( src1pnt, usrHd->flowDiag+p1Off, usrHd->lesP+p1Off, &nSrc1Dims, &nDofs, &nPs, &nP1s, &(usrHd->nNodes) ) ; fMtxVdimVecMult( src2pnt, usrHd->flowDiag+p2Off, usrHd->lesP+p2Off, &nSrc2Dims, &nDofs, &nPs, &nP2s, &(usrHd->nNodes) ) ; commOut ( usrHd->lesP, usrHd->ilwork, &nPs, usrHd->iper, usrHd->iBC, usrHd->BC ); fLesSparseApNGtC( usrHd->colm, usrHd->rowp, usrHd->lhsP, usrHd->lesP, usrHd->lesQ, &(usrHd->nNodes),&(usrHd->nnz_tot)); commIn ( usrHd->lesQ, usrHd->ilwork, &nQs, usrHd->iper, usrHd->iBC, usrHd->BC ); dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; fMtxVdimVecMult( usrHd->lesQ, usrHd->flowDiag+qOff, dstpnt, &nQs, &nDofs, &nDstDims, &nQs, &(usrHd->nNodes) ) ; } /*----------------------------------------------------------------------- * * lesApFull * * operation : K * Du + G * Dp = lesQ(:,nenl,1:3) * -G^t * Du + C * Dp = lesQ(:,nenl,4:4) * *----------------------------------------------------------------------- */ void lesApFull ( UsrHd usrHd, Integer srcId, Integer dstId, Integer nSrcDims, Integer srcOff, Integer nDstDims, Integer dstOff ) { char* funcName = "lesApFull" ; /* function name */ Integer pOff ; /* Diag offset for P */ Integer qOff ; /* Diag offset for Q */ Integer nDofs ; /* No. of Dofs */ Integer nPs ; /* No. of P dimension */ Integer nQs ; /* No. of Q dimension */ Real* srcpnt ; /* Source */ Real* dstpnt ; /* Destination */ nDofs = 4 ; nPs = 4 ; nQs = 4 ; pOff = 0 * usrHd->nNodes ; qOff = 0 * usrHd->nNodes ; srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; fMtxVdimVecMult( srcpnt, usrHd->flowDiag+pOff, usrHd->lesP, &nSrcDims, &nDofs, &nPs, &nPs, &(usrHd->nNodes) ) ; commOut ( usrHd->lesP, usrHd->ilwork, &nPs, usrHd->iper, usrHd->iBC, usrHd->BC ); fLesSparseApFull( usrHd->colm, usrHd->rowp, usrHd->lhsK, usrHd->lhsP, usrHd->lesP, usrHd->lesQ, &(usrHd->nNodes),&(usrHd->nnz_tot)); commIn ( usrHd->lesQ, usrHd->ilwork, &nQs, usrHd->iper, usrHd->iBC, usrHd->BC ); dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; fMtxVdimVecMult( usrHd->lesQ, usrHd->flowDiag+qOff, dstpnt, &nQs, &nDofs, &nDstDims, &nQs, &(usrHd->nNodes) ) ; } /*----------------------------------------------------------------------- * * lesApSclr * * operation : M(:,nenl,nenl) * Ds(:,nenl,1) = lesQ(:,nenl,1) * *----------------------------------------------------------------------- */ void lesApSclr ( UsrHd usrHd, Integer srcId, Integer dstId, Integer nSrcDims, Integer srcOff, Integer nDstDims, Integer dstOff ) { char* funcName = "lesApSclr" ; /* function name */ Integer pOff ; /* Diag offset for P */ Integer qOff ; /* Diag offset for Q */ Integer nDofs ; /* No. of Dofs */ Integer nPs ; /* No. of P dimension */ Integer nQs ; /* No. of Q dimension */ Real* srcpnt ; /* Source */ Real* dstpnt ; /* Destination */ Integer lhsStiffFlag ; double sclrRegFct ; nDofs = 1 ; nPs = 1 ; nQs = 1 ; pOff = 0 ; qOff = 0 ; lhsStiffFlag = 0 ; sclrRegFct = 0 ; srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; fMtxVdimVecMult ( srcpnt, usrHd->sclrDiag+pOff, usrHd->lesP, &nSrcDims, &nDofs, &nPs, &nPs, &(usrHd->nNodes) ) ; commOut ( usrHd->lesP, usrHd->ilwork, &nPs, usrHd->iper, usrHd->iBC, usrHd->BC ); fLesSparseApSclr( usrHd->colm, usrHd->rowp, usrHd->lhsS, usrHd->lesP, usrHd->lesQ, &(usrHd->nNodes),&(usrHd->nnz_tot)); commIn ( usrHd->lesQ, usrHd->ilwork, &nQs, usrHd->iper, usrHd->iBC, usrHd->BC ); if ( lhsStiffFlag && sclrRegFct != 0 ) { fMtxVdimVecMult ( usrHd->lesQ, usrHd->sclrDiag+qOff, usrHd->lesP, &nQs, &nDofs, &nDstDims, &nQs, &(usrHd->nNodes) ) ; flesDaxpy ( srcpnt, usrHd->lesP, &sclrRegFct, &nDstDims, &(usrHd->nNodes) ) ; flesCp ( usrHd->lesP, dstpnt, &nDstDims, &(usrHd->nNodes) ) ; } else { fMtxVdimVecMult ( usrHd->lesQ, usrHd->sclrDiag+qOff, dstpnt, &nQs, &nDofs, &nDstDims, &nQs, &(usrHd->nNodes) ) ; } } /********************************************************************* * lesPrecPPE * outer routine to solve PPE *******************************************************************/ void lesPrecPPE(UsrHd usrHd, Integer srcId, Integer dstId, Integer nSrcDims, Integer srcOff, Integer nDstDims, Integer dstOff) { Real* srcpnt ; /* Source R */ Real* dstpnt ; /* Destination Z */ srcpnt = usrPointer ( usrHd, srcId, srcOff, nSrcDims ) ; dstpnt = usrPointer ( usrHd, dstId, dstOff, nDstDims ) ; #ifdef AMG ramg_interface( usrHd->colm, usrHd->rowp,usrHd->lhsK,usrHd->lhsP,usrHd->flowDiag, srcpnt,dstpnt, usrHd->ilwork,usrHd->BC,usrHd->iBC,usrHd->iper ); #endif return; }