Actual source code: blockinvert.h

      1: /*
      2:     Kernels used in sparse ILU (and LU) and in the resulting triangular
      3:  solves. These are for block algorithms where the block sizes are on
      4:  the order of 2-6+.

      6:     There are TWO versions of the macros below.
      7:     1) standard for MatScalar == PetscScalar use the standard BLAS for
      8:        block size larger than 7 and
      9:     2) handcoded Fortran single precision for the matrices, since BLAS
     10:        does not have some arguments in single and some in double.

     12: */
     15: #include <petscblaslapack.h>

     17: /*
     18:       These are C kernels,they are contained in
     19:    src/mat/impls/baij/seq
     20: */

     22: PETSC_INTERN PetscErrorCode PetscLINPACKgefa(MatScalar*,PetscInt,PetscInt*,PetscBool,PetscBool*);
     23: PETSC_INTERN PetscErrorCode PetscLINPACKgedi(MatScalar*,PetscInt,PetscInt*,MatScalar*);
     24: PETSC_EXTERN PetscErrorCode PetscKernel_A_gets_inverse_A_2(MatScalar*,PetscReal,PetscBool,PetscBool*);
     25: PETSC_EXTERN PetscErrorCode PetscKernel_A_gets_inverse_A_3(MatScalar*,PetscReal,PetscBool,PetscBool*);

     27: #define PetscKernel_A_gets_inverse_A_4_nopivot(mat) 0; \
     28:   { \
     29:     MatScalar d, di; \
     30: \
     31:     di       = mat[0]; \
     32:     mat[0]   = d = 1.0 / di; \
     33:     mat[4]  *= -d; \
     34:     mat[8]  *= -d; \
     35:     mat[12] *= -d; \
     36:     mat[1]  *= d; \
     37:     mat[2]  *= d; \
     38:     mat[3]  *= d; \
     39:     mat[5]  += mat[4] * mat[1] * di; \
     40:     mat[6]  += mat[4] * mat[2] * di; \
     41:     mat[7]  += mat[4] * mat[3] * di; \
     42:     mat[9]  += mat[8] * mat[1] * di; \
     43:     mat[10] += mat[8] * mat[2] * di; \
     44:     mat[11] += mat[8] * mat[3] * di; \
     45:     mat[13] += mat[12] * mat[1] * di; \
     46:     mat[14] += mat[12] * mat[2] * di; \
     47:     mat[15] += mat[12] * mat[3] * di; \
     48:     di       = mat[5]; \
     49:     mat[5]   = d = 1.0 / di; \
     50:     mat[1]  *= -d; \
     51:     mat[9]  *= -d; \
     52:     mat[13] *= -d; \
     53:     mat[4]  *= d; \
     54:     mat[6]  *= d; \
     55:     mat[7]  *= d; \
     56:     mat[0]  += mat[1] * mat[4] * di; \
     57:     mat[2]  += mat[1] * mat[6] * di; \
     58:     mat[3]  += mat[1] * mat[7] * di; \
     59:     mat[8]  += mat[9] * mat[4] * di; \
     60:     mat[10] += mat[9] * mat[6] * di; \
     61:     mat[11] += mat[9] * mat[7] * di; \
     62:     mat[12] += mat[13] * mat[4] * di; \
     63:     mat[14] += mat[13] * mat[6] * di; \
     64:     mat[15] += mat[13] * mat[7] * di; \
     65:     di       = mat[10]; \
     66:     mat[10]  = d = 1.0 / di; \
     67:     mat[2]  *= -d; \
     68:     mat[6]  *= -d; \
     69:     mat[14] *= -d; \
     70:     mat[8]  *= d; \
     71:     mat[9]  *= d; \
     72:     mat[11] *= d; \
     73:     mat[0]  += mat[2] * mat[8] * di; \
     74:     mat[1]  += mat[2] * mat[9] * di; \
     75:     mat[3]  += mat[2] * mat[11] * di; \
     76:     mat[4]  += mat[6] * mat[8] * di; \
     77:     mat[5]  += mat[6] * mat[9] * di; \
     78:     mat[7]  += mat[6] * mat[11] * di; \
     79:     mat[12] += mat[14] * mat[8] * di; \
     80:     mat[13] += mat[14] * mat[9] * di; \
     81:     mat[15] += mat[14] * mat[11] * di; \
     82:     di       = mat[15]; \
     83:     mat[15]  = d = 1.0 / di; \
     84:     mat[3]  *= -d; \
     85:     mat[7]  *= -d; \
     86:     mat[11] *= -d; \
     87:     mat[12] *= d; \
     88:     mat[13] *= d; \
     89:     mat[14] *= d; \
     90:     mat[0]  += mat[3] * mat[12] * di; \
     91:     mat[1]  += mat[3] * mat[13] * di; \
     92:     mat[2]  += mat[3] * mat[14] * di; \
     93:     mat[4]  += mat[7] * mat[12] * di; \
     94:     mat[5]  += mat[7] * mat[13] * di; \
     95:     mat[6]  += mat[7] * mat[14] * di; \
     96:     mat[8]  += mat[11] * mat[12] * di; \
     97:     mat[9]  += mat[11] * mat[13] * di; \
     98:     mat[10] += mat[11] * mat[14] * di; \
     99:   }

    101: PETSC_EXTERN PetscErrorCode PetscKernel_A_gets_inverse_A_4(MatScalar*,PetscReal,PetscBool,PetscBool*);
    102: # if defined(PETSC_HAVE_SSE)
    103: PETSC_EXTERN PetscErrorCode PetscKernel_A_gets_inverse_A_4_SSE(MatScalar*);
    104: # endif
    105: PETSC_EXTERN PetscErrorCode PetscKernel_A_gets_inverse_A_5(MatScalar*,PetscInt*,MatScalar*,PetscReal,PetscBool,PetscBool*);
    106: PETSC_EXTERN PetscErrorCode PetscKernel_A_gets_inverse_A_6(MatScalar*,PetscReal,PetscBool,PetscBool*);
    107: PETSC_EXTERN PetscErrorCode PetscKernel_A_gets_inverse_A_7(MatScalar*,PetscReal,PetscBool,PetscBool*);
    108: PETSC_EXTERN PetscErrorCode PetscKernel_A_gets_inverse_A_9(MatScalar*,PetscReal,PetscBool,PetscBool*);
    109: PETSC_EXTERN PetscErrorCode PetscKernel_A_gets_inverse_A_15(MatScalar*,PetscInt*,MatScalar*,PetscReal,PetscBool,PetscBool*);

    111: /*
    112:     A = inv(A)    A_gets_inverse_A

    114:    A      - square bs by bs array stored in column major order
    115:    pivots - integer work array of length bs
    116:    W      -  bs by 1 work array
    117: */
    118: #define PetscKernel_A_gets_inverse_A(bs,A,pivots,W,allowzeropivot,zeropivotdetected) (PetscLINPACKgefa((A),(bs),(pivots),(allowzeropivot),(zeropivotdetected)) || PetscLINPACKgedi((A),(bs),(pivots),(W)))

    120: /* -----------------------------------------------------------------------*/

    122: #if !defined(PETSC_USE_REAL_MAT_SINGLE)
    123: /*
    124:         Version that calls the BLAS directly
    125: */
    126: /*
    127:       A = A * B   A_gets_A_times_B

    129:    A, B - square bs by bs arrays stored in column major order
    130:    W    - square bs by bs work array

    132: */
    133: #define PetscKernel_A_gets_A_times_B(bs,A,B,W) \
    134:   {                                              \
    135:     PetscBLASInt   _bbs;                         \
    136:     PetscScalar    _one = 1.0,_zero = 0.0;       \
    137:     PetscErrorCode _ierr;                        \
    138:     _PetscBLASIntCast(bs,&_bbs); \
    139:     _PetscArraycpy((W),(A),(bs)*(bs));CHKERRQ(_ierr); \
    143:     PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&(_bbs),&(_bbs),&(_bbs),&_one,(W),&(_bbs),(B),&(_bbs),&_zero,(A),&(_bbs))); \
    144:   }

    146: /*

    148:     A = A - B * C  A_gets_A_minus_B_times_C

    150:    A, B, C - square bs by bs arrays stored in column major order
    151: */
    152: #define PetscKernel_A_gets_A_minus_B_times_C(bs,A,B,C) \
    153:   {                                                      \
    154:     PetscScalar    _mone = -1.0,_one = 1.0;              \
    155:     PetscBLASInt   _bbs;                                 \
    156:     PetscErrorCode _ierr;                                \
    157:     _PetscBLASIntCast(bs,&_bbs);CHKERRQ(_ierr);   \
    161:     PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&(_bbs),&(_bbs),&(_bbs),&_mone,(B),&(_bbs),(C),&(_bbs),&_one,(A),&(_bbs))); \
    162:   }

    164: /*

    166:     A = A + B^T * C  A_gets_A_plus_Btranspose_times_C

    168:    A, B, C - square bs by bs arrays stored in column major order
    169: */
    170: #define PetscKernel_A_gets_A_plus_Btranspose_times_C(bs,A,B,C) \
    171:   {                                                              \
    172:     PetscScalar    _one = 1.0;                                   \
    173:     PetscBLASInt   _bbs;                                         \
    174:     PetscErrorCode _ierr;                                        \
    175:     _PetscBLASIntCast(bs,&_bbs);CHKERRQ(_ierr);           \
    179:     PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&(_bbs),&(_bbs),&(_bbs),&_one,(B),&(_bbs),(C),&(_bbs),&_one,(A),&(_bbs))); \
    180:   }

    182: /*
    183:     v = v + A^T w  v_gets_v_plus_Atranspose_times_w

    185:    v - array of length bs
    186:    A - square bs by bs array
    187:    w - array of length bs
    188: */
    189: #define  PetscKernel_v_gets_v_plus_Atranspose_times_w(bs,v,A,w) \
    190:   {                                                               \
    191:     PetscScalar    _one  = 1.0;                                    \
    192:     PetscBLASInt   _ione = 1, _bbs;                               \
    193:     PetscErrorCode _ierr;                                         \
    194:     _PetscBLASIntCast(bs,&_bbs);CHKERRQ(_ierr);            \
    198:     PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&(_bbs),&(_bbs),&_one,A,&(_bbs),w,&_ione,&_one,v,&_ione)); \
    199:   }

    201: /*
    202:     v = v - A w  v_gets_v_minus_A_times_w

    204:    v - array of length bs
    205:    A - square bs by bs array
    206:    w - array of length bs
    207: */
    208: #define  PetscKernel_v_gets_v_minus_A_times_w(bs,v,A,w) \
    209:   {                                                       \
    210:     PetscScalar    _mone = -1.0,_one = 1.0;                 \
    211:     PetscBLASInt   _ione = 1,_bbs;                         \
    212:     PetscErrorCode _ierr;                                         \
    213:     _PetscBLASIntCast(bs,&_bbs);CHKERRQ(_ierr);            \
    217:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&(_bbs),&(_bbs),&_mone,A,&(_bbs),w,&_ione,&_one,v,&_ione)); \
    218:   }

    220: /*
    221:     v = v - A w  v_gets_v_minus_transA_times_w

    223:    v - array of length bs
    224:    A - square bs by bs array
    225:    w - array of length bs
    226: */
    227: #define  PetscKernel_v_gets_v_minus_transA_times_w(bs,v,A,w) \
    228:   {                                                            \
    229:     PetscScalar    _mone = -1.0,_one = 1.0;                      \
    230:     PetscBLASInt   _ione = 1,_bbs;                              \
    231:     PetscErrorCode _ierr;                                      \
    232:     _PetscBLASIntCast(bs,&_bbs);CHKERRQ(_ierr);         \
    236:     PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&(_bbs),&(_bbs),&_mone,A,&(_bbs),w,&_ione,&_one,v,&_ione)); \
    237:   }

    239: /*
    240:     v = v + A w  v_gets_v_plus_A_times_w

    242:    v - array of length bs
    243:    A - square bs by bs array
    244:    w - array of length bs
    245: */
    246: #define  PetscKernel_v_gets_v_plus_A_times_w(bs,v,A,w) \
    247:   {                                                      \
    248:     PetscScalar    _one  = 1.0;                             \
    249:     PetscBLASInt   _ione = 1,_bbs;                         \
    250:     PetscErrorCode _ierr;                                      \
    251:     _PetscBLASIntCast(bs,&_bbs);CHKERRQ(_ierr);         \
    255:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&(_bbs),&(_bbs),&_one,A,&(_bbs),w,&_ione,&_one,v,&_ione)); \
    256:   }

    258: /*
    259:     v = v + A w  v_gets_v_plus_Ar_times_w

    261:    v - array of length bs
    262:    A - square bs by bs array
    263:    w - array of length bs
    264: */
    265: #define  PetscKernel_w_gets_w_plus_Ar_times_v(bs,ncols,v,A,w) \
    266:   {                                                             \
    267:     PetscScalar    _one  = 1.0;                                    \
    268:     PetscBLASInt   _ione = 1,_bbs,_bncols;                        \
    269:     PetscErrorCode _ierr;                                      \
    270:     _PetscBLASIntCast(bs,&_bbs);CHKERRQ(_ierr);          \
    271:     _PetscBLASIntCast(ncols,&_bncols);CHKERRQ(_ierr);    \
    275:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&(_bbs),&(_bncols),&_one,A,&(_bbs),v,&_ione,&_one,w,&_ione)); \
    276:   }

    278: /*
    279:     w <- w - A v  w_gets_w_minus_Ar_times_v

    281:    v - array of length ncol
    282:    A(bs,ncols)
    283:    w - array of length bs
    284: */
    285: #define  PetscKernel_w_gets_w_minus_Ar_times_v(bs,ncols,w,A,v) \
    286:   {                                                              \
    287:     PetscScalar    _one  = 1.0,_mone = -1.0;                        \
    288:     PetscBLASInt   _ione = 1,_bbs,_bncols;                         \
    289:     PetscErrorCode _ierr;                                      \
    290:     _PetscBLASIntCast(bs,&_bbs);CHKERRQ(_ierr);          \
    291:     _PetscBLASIntCast(ncols,&_bncols);CHKERRQ(_ierr);    \
    295:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&(_bbs),&(_bncols),&_mone,A,&(_bbs),v,&_ione,&_one,w,&_ione)); \
    296:   }

    298: /*
    299:     w = A v   w_gets_A_times_v

    301:    v - array of length bs
    302:    A - square bs by bs array
    303:    w - array of length bs
    304: */
    305: #define PetscKernel_w_gets_A_times_v(bs,v,A,w) \
    306:   {                                              \
    307:     PetscScalar    _zero = 0.0,_one = 1.0;         \
    308:     PetscBLASInt   _ione = 1,_bbs;                 \
    309:     PetscErrorCode _ierr;                          \
    310:     _PetscBLASIntCast(bs,&_bbs);CHKERRQ(_ierr);\
    314:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&(_bbs),&(_bbs),&_one,A,&(_bbs),v,&_ione,&_zero,w,&_ione)); \
    315:   }

    317: /*
    318:     w = A' v   w_gets_transA_times_v

    320:    v - array of length bs
    321:    A - square bs by bs array
    322:    w - array of length bs
    323: */
    324: #define PetscKernel_w_gets_transA_times_v(bs,v,A,w)     \
    325:   {                                                       \
    326:     PetscScalar    _zero = 0.0,_one = 1.0;                  \
    327:     PetscBLASInt   _ione = 1,_bbs;                        \
    328:     PetscErrorCode _ierr;                          \
    329:     _PetscBLASIntCast(bs,&_bbs);CHKERRQ(_ierr); \
    333:     PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&(_bbs),&(_bbs),&_one,A,&(_bbs),v,&_ione,&_zero,w,&_ione)); \
    334:   }

    336: /*
    337:         z = A*x
    338: */
    339: #define PetscKernel_w_gets_Ar_times_v(bs,ncols,x,A,z) \
    340:   {                                                     \
    341:     PetscScalar    _one  = 1.0,_zero = 0.0;                 \
    342:     PetscBLASInt   _ione = 1,_bbs,_bncols;        \
    343:     PetscErrorCode _ierr; \
    344:     _PetscBLASIntCast(bs,&_bbs);CHKERRQ(_ierr); \
    345:     _PetscBLASIntCast(ncols,&_bncols);CHKERRQ(_ierr); \
    349:     PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&(_bbs),&_bncols,&_one,A,&(_bbs),x,&_ione,&_zero,z,&_ione)); \
    350:   }

    352: /*
    353:         z = trans(A)*x
    354: */
    355: #define PetscKernel_w_gets_w_plus_trans_Ar_times_v(bs,ncols,x,A,z) \
    356:   {                                                                  \
    357:     PetscScalar    _one  = 1.0;                                         \
    358:     PetscBLASInt   _ione = 1,_bbs,_bncols;                             \
    359:     PetscErrorCode _ierr; \
    360:     _PetscBLASIntCast(bs,&_bbs);CHKERRQ(_ierr); \
    361:     _PetscBLASIntCast(ncols,&_bncols);CHKERRQ(_ierr); \
    365:     PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&_bbs,&_bncols,&_one,A,&_bbs,x,&_ione,&_one,z,&_ione)); \
    366:   }

    368: #else
    369: /*
    370:        Version that calls Fortran routines; can handle different precision
    371:    of matrix (array) and vectors
    372: */
    373: /*
    374:      These are Fortran kernels: They replace certain BLAS routines but
    375:    have some arguments that may be single precision,rather than double
    376:    These routines are provided in src/fortran/kernels/sgemv.F
    377:    They are pretty pitiful but get the job done. The intention is
    378:    that for important block sizes (currently 1,2,3,4,5,6,7) custom inlined
    379:    code is used.
    380: */

    382: #if defined(PETSC_HAVE_FORTRAN_CAPS)
    383: #define msgemv_  MSGEMV
    384: #define msgemvp_ MSGEMVP
    385: #define msgemvm_ MSGEMVM
    386: #define msgemvt_ MSGEMVT
    387: #define msgemmi_ MSGEMMI
    388: #define msgemm_  MSGEMM
    389: #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
    390: #define msgemv_  msgemv
    391: #define msgemvp_ msgemvp
    392: #define msgemvm_ msgemvm
    393: #define msgemvt_ msgemvt
    394: #define msgemmi_ msgemmi
    395: #define msgemm_  msgemm
    396: #endif

    398: PETSC_EXTERN void msgemv_(PetscInt*,PetscInt*,MatScalar*,PetscScalar*,PetscScalar*);
    399: PETSC_EXTERN void msgemvp_(PetscInt*,PetscInt*,MatScalar*,PetscScalar*,PetscScalar*);
    400: PETSC_EXTERN void msgemvm_(PetscInt*,PetscInt*,MatScalar*,PetscScalar*,PetscScalar*);
    401: PETSC_EXTERN void msgemvt_(PetscInt*,PetscInt*,MatScalar*,PetscScalar*,PetscScalar*);
    402: PETSC_EXTERN void msgemmi_(PetscInt*,MatScalar*,MatScalar*,MatScalar*);
    403: PETSC_EXTERN void msgemm_(PetscInt*,MatScalar*,MatScalar*,MatScalar*);

    405: /*
    406:       A = A * B   A_gets_A_times_B

    408:    A, B - square bs by bs arrays stored in column major order
    409:    W    - square bs by bs work array

    411: */
    412: #define PetscKernel_A_gets_A_times_B(bs,A,B,W) \
    413:   {                                              \
    414:     PetscErrorCode _PetscArraycpy((W),(A),(bs)*(bs));CHKERRQ(_ierr); \
    415:     msgemmi_(&bs,A,B,W); \
    416:   }

    418: /*

    420:     A = A - B * C  A_gets_A_minus_B_times_C

    422:    A, B, C - square bs by bs arrays stored in column major order
    423: */
    424: #define PetscKernel_A_gets_A_minus_B_times_C(bs,A,B,C) \
    425:   {                                                      \
    426:     msgemm_(&bs,A,B,C);                                  \
    427:   }

    429: /*
    430:     v = v - A w  v_gets_v_minus_A_times_w

    432:    v - array of length bs
    433:    A - square bs by bs array
    434:    w - array of length bs
    435: */
    436: #define  PetscKernel_v_gets_v_minus_A_times_w(bs,v,A,w) \
    437:   {                                                       \
    438:     msgemvm_(&bs,&bs,A,w,v);                              \
    439:   }

    441: /*
    442:     v = v + A w  v_gets_v_plus_A_times_w

    444:    v - array of length bs
    445:    A - square bs by bs array
    446:    w - array of length bs
    447: */
    448: #define  PetscKernel_v_gets_v_plus_A_times_w(bs,v,A,w) \
    449:   {                                                      \
    450:     msgemvp_(&bs,&bs,A,w,v);                             \
    451:   }

    453: /*
    454:     v = v + A w  v_gets_v_plus_Ar_times_w

    456:    v - array of length bs
    457:    A - square bs by bs array
    458:    w - array of length bs
    459: */
    460: #define  PetscKernel_w_gets_w_plus_Ar_times_v(bs,ncol,v,A,w) \
    461:   {                                                            \
    462:     msgemvp_(&bs,&ncol,A,v,w);                                 \
    463:   }

    465: /*
    466:     w = A v   w_gets_A_times_v

    468:    v - array of length bs
    469:    A - square bs by bs array
    470:    w - array of length bs
    471: */
    472: #define PetscKernel_w_gets_A_times_v(bs,v,A,w) \
    473:   {  \
    474:     msgemv_(&bs,&bs,A,v,w); \
    475:   }

    477: /*
    478:         z = A*x
    479: */
    480: #define PetscKernel_w_gets_Ar_times_v(bs,ncols,x,A,z) \
    481:   { \
    482:     msgemv_(&bs,&ncols,A,x,z); \
    483:   }

    485: /*
    486:         z = trans(A)*x
    487: */
    488: #define PetscKernel_w_gets_w_plus_trans_Ar_times_v(bs,ncols,x,A,z) \
    489:   { \
    490:     msgemvt_(&bs,&ncols,A,x,z); \
    491:   }

    493: /* These do not work yet */
    494: #define PetscKernel_A_gets_A_plus_Btranspose_times_C(bs,A,B,C)
    495: #define PetscKernel_v_gets_v_plus_Atranspose_times_w(bs,v,A,w)

    498: #endif

    500: #endif