1 /* 2 PetscKernelAXPY - X = X + alpha * Y 3 4 Input Parameters: 5 + X, Y - arrays 6 . alpha - scalar 7 - n - length of arrays 8 9 Also PetscKernelAXPY2(), PetscKernelAXPY3(), PetscKernelAXPY4() 10 11 */ 12 13 #pragma once 14 15 #if defined(PETSC_USE_FORTRAN_KERNEL_MAXPY) 16 #if defined(PETSC_HAVE_FORTRAN_CAPS) 17 #define fortranmaxpy4_ FORTRANMAXPY4 18 #define fortranmaxpy3_ FORTRANMAXPY3 19 #define fortranmaxpy2_ FORTRANMAXPY2 20 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 21 #define fortranmaxpy4_ fortranmaxpy4 22 #define fortranmaxpy3_ fortranmaxpy3 23 #define fortranmaxpy2_ fortranmaxpy2 24 #endif 25 PETSC_EXTERN void fortranmaxpy4_(void *, const void *, const void *, const void *, const void *, const void *, const void *, const void *, const void *, const PetscInt *); 26 PETSC_EXTERN void fortranmaxpy3_(void *, const void *, const void *, const void *, const void *, const void *, const void *, const PetscInt *); 27 PETSC_EXTERN void fortranmaxpy2_(void *, const void *, const void *, const void *, const void *, const PetscInt *); 28 #endif 29 #include <petscblaslapack.h> 30 31 #if defined(PETSC_USE_FORTRAN_KERNEL_MAXPY) 32 #define PetscKernelAXPY(U, a1, p1, n) \ 33 { \ 34 PetscBLASInt one = 1; \ 35 PetscBLASInt nn = (PetscBLASInt)n; \ 36 PetscCallBLAS("BLASaxpy", BLASaxpy_(&nn, &a1, p1, &one, U, &one)); \ 37 } 38 #define PetscKernelAXPY2(U, a1, a2, p1, p2, n) \ 39 { \ 40 fortranmaxpy2_(U, &a1, &a2, p1, p2, &n); \ 41 } 42 #define PetscKernelAXPY3(U, a1, a2, a3, p1, p2, p3, n) \ 43 { \ 44 fortranmaxpy3_(U, &a1, &a2, &a3, p1, p2, p3, &n); \ 45 } 46 #define PetscKernelAXPY4(U, a1, a2, a3, a4, p1, p2, p3, p4, n) \ 47 { \ 48 fortranmaxpy4_(U, &a1, &a2, &a3, &a4, p1, p2, p3, p4, &n); \ 49 } 50 51 #elif defined(PETSC_USE_UNROLL_KERNELS) 52 53 #define PetscKernelAXPY(U, Alpha, P, n) \ 54 { \ 55 switch (n & 0x3) { \ 56 case 3: \ 57 *U++ += Alpha * *P++; \ 58 case 2: \ 59 *U++ += Alpha * *P++; \ 60 case 1: \ 61 *U++ += Alpha * *P++; \ 62 n -= 4; \ 63 case 0: \ 64 break; \ 65 } \ 66 while (n > 0) { \ 67 U[0] += Alpha * P[0]; \ 68 U[1] += Alpha * P[1]; \ 69 U[2] += Alpha * P[2]; \ 70 U[3] += Alpha * P[3]; \ 71 U += 4; \ 72 P += 4; \ 73 n -= 4; \ 74 } \ 75 } 76 #define PetscKernelAXPY2(U, a1, a2, p1, p2, n) \ 77 { \ 78 switch (n & 0x3) { \ 79 case 3: \ 80 *U++ += a1 * *p1++ + a2 * *p2++; \ 81 case 2: \ 82 *U++ += a1 * *p1++ + a2 * *p2++; \ 83 case 1: \ 84 *U++ += a1 * *p1++ + a2 * *p2++; \ 85 n -= 4; \ 86 case 0: \ 87 break; \ 88 } \ 89 while (n > 0) { \ 90 U[0] += a1 * p1[0] + a2 * p2[0]; \ 91 U[1] += a1 * p1[1] + a2 * p2[1]; \ 92 U[2] += a1 * p1[2] + a2 * p2[2]; \ 93 U[3] += a1 * p1[3] + a2 * p2[3]; \ 94 U += 4; \ 95 p1 += 4; \ 96 p2 += 4; \ 97 n -= 4; \ 98 } \ 99 } 100 #define PetscKernelAXPY3(U, a1, a2, a3, p1, p2, p3, n) \ 101 { \ 102 switch (n & 0x3) { \ 103 case 3: \ 104 *U++ += a1 * *p1++ + a2 * *p2++ + a3 * *p3++; \ 105 case 2: \ 106 *U++ += a1 * *p1++ + a2 * *p2++ + a3 * *p3++; \ 107 case 1: \ 108 *U++ += a1 * *p1++ + a2 * *p2++ + a3 * *p3++; \ 109 n -= 4; \ 110 case 0: \ 111 break; \ 112 } \ 113 while (n > 0) { \ 114 U[0] += a1 * p1[0] + a2 * p2[0] + a3 * p3[0]; \ 115 U[1] += a1 * p1[1] + a2 * p2[1] + a3 * p3[1]; \ 116 U[2] += a1 * p1[2] + a2 * p2[2] + a3 * p3[2]; \ 117 U[3] += a1 * p1[3] + a2 * p2[3] + a3 * p3[3]; \ 118 U += 4; \ 119 p1 += 4; \ 120 p2 += 4; \ 121 p3 += 4; \ 122 n -= 4; \ 123 } \ 124 } 125 #define PetscKernelAXPY4(U, a1, a2, a3, a4, p1, p2, p3, p4, n) \ 126 { \ 127 switch (n & 0x3) { \ 128 case 3: \ 129 *U++ += a1 * *p1++ + a2 * *p2++ + a3 * *p3++ + a4 * *p4++; \ 130 case 2: \ 131 *U++ += a1 * *p1++ + a2 * *p2++ + a3 * *p3++ + a4 * *p4++; \ 132 case 1: \ 133 *U++ += a1 * *p1++ + a2 * *p2++ + a3 * *p3++ + a4 * *p4++; \ 134 n -= 4; \ 135 case 0: \ 136 break; \ 137 } \ 138 while (n > 0) { \ 139 U[0] += a1 * p1[0] + a2 * p2[0] + a3 * p3[0] + a4 * p4[0]; \ 140 U[1] += a1 * p1[1] + a2 * p2[1] + a3 * p3[1] + a4 * p4[1]; \ 141 U[2] += a1 * p1[2] + a2 * p2[2] + a3 * p3[2] + a4 * p4[2]; \ 142 U[3] += a1 * p1[3] + a2 * p2[3] + a3 * p3[3] + a4 * p4[3]; \ 143 U += 4; \ 144 p1 += 4; \ 145 p2 += 4; \ 146 p3 += 4; \ 147 p4 += 4; \ 148 n -= 4; \ 149 } \ 150 } 151 152 #elif defined(PETSC_USE_WHILE_KERNELS) 153 154 #define PetscKernelAXPY(U, a1, p1, n) \ 155 { \ 156 while (n--) *U++ += a1 * *p1++; \ 157 } 158 #define PetscKernelAXPY2(U, a1, a2, p1, p2, n) \ 159 { \ 160 while (n--) *U++ += a1 * *p1++ + a2 * *p2++; \ 161 } 162 #define PetscKernelAXPY3(U, a1, a2, a3, p1, p2, p3, n) \ 163 { \ 164 while (n--) *U++ += a1 * *p1++ + a2 * *p2++ + a3 * *p3++; \ 165 } 166 #define PetscKernelAXPY4(U, a1, a2, a3, a4, p1, p2, p3, p4, n) \ 167 { \ 168 while (n--) *U++ += a1 * *p1++ + a2 * *p2++ + a3 * *p3++ + a4 * *p4++; \ 169 } 170 171 #elif defined(PETSC_USE_BLAS_KERNELS) 172 173 #define PetscKernelAXPY(U, a1, p1, n) \ 174 { \ 175 PetscBLASInt one = 1; \ 176 PetscBLASInt nn = (PetscBLASInt)n; \ 177 PetscCallBLAS("BLASaxpy", BLASaxpy_(&nn, &a1, p1, &one, U, &one)); \ 178 } 179 #define PetscKernelAXPY2(U, a1, a2, p1, p2, n) \ 180 { \ 181 PetscKernelAXPY(U, a1, p1, n); \ 182 PetscKernelAXPY(U, a2, p2, n); \ 183 } 184 #define PetscKernelAXPY3(U, a1, a2, a3, p1, p2, p3, n) \ 185 { \ 186 PetscKernelAXPY2(U, a1, a2, p1, p2, n); \ 187 PetscKernelAXPY(U, a3, p3, n); \ 188 } 189 #define PetscKernelAXPY4(U, a1, a2, a3, a4, p1, p2, p3, p4, n) \ 190 { \ 191 PetscKernelAXPY2(U, a1, a2, p1, p2, n); \ 192 PetscKernelAXPY2(U, a3, a4, p3, p4, n); \ 193 } 194 195 #else 196 197 #define PetscKernelAXPY(U, a1, p1, n) \ 198 do { \ 199 const PetscInt _n = n; \ 200 const PetscScalar _a1 = a1; \ 201 const PetscScalar *PETSC_RESTRICT _p1 = p1; \ 202 PetscScalar *PETSC_RESTRICT _U = U; \ 203 PetscPragmaUseOMPKernels(parallel for) \ 204 for (PetscInt __i = 0; __i < _n; __i++) _U[__i] += _a1 * _p1[__i]; \ 205 } while (0) 206 #define PetscKernelAXPY2(U, a1, a2, p1, p2, n) \ 207 do { \ 208 const PetscInt _n = n; \ 209 const PetscScalar _a1 = a1; \ 210 const PetscScalar _a2 = a2; \ 211 const PetscScalar *PETSC_RESTRICT _p1 = p1; \ 212 const PetscScalar *PETSC_RESTRICT _p2 = p2; \ 213 PetscScalar *PETSC_RESTRICT _U = U; \ 214 PetscPragmaUseOMPKernels(parallel for) \ 215 for (PetscInt __i = 0; __i < _n; __i++) _U[__i] += _a1 * _p1[__i] + _a2 * _p2[__i]; \ 216 } while (0) 217 #define PetscKernelAXPY3(U, a1, a2, a3, p1, p2, p3, n) \ 218 do { \ 219 const PetscInt _n = n; \ 220 const PetscScalar _a1 = a1; \ 221 const PetscScalar _a2 = a2; \ 222 const PetscScalar _a3 = a3; \ 223 const PetscScalar *PETSC_RESTRICT _p1 = p1; \ 224 const PetscScalar *PETSC_RESTRICT _p2 = p2; \ 225 const PetscScalar *PETSC_RESTRICT _p3 = p3; \ 226 PetscScalar *PETSC_RESTRICT _U = U; \ 227 PetscPragmaUseOMPKernels(parallel for) \ 228 for (PetscInt __i = 0; __i < _n; __i++) _U[__i] += _a1 * _p1[__i] + _a2 * _p2[__i] + _a3 * _p3[__i]; \ 229 } while (0) 230 #define PetscKernelAXPY4(U, a1, a2, a3, a4, p1, p2, p3, p4, n) \ 231 do { \ 232 const PetscInt _n = n; \ 233 const PetscScalar _a1 = a1; \ 234 const PetscScalar _a2 = a2; \ 235 const PetscScalar _a3 = a3; \ 236 const PetscScalar _a4 = a4; \ 237 const PetscScalar *PETSC_RESTRICT _p1 = p1; \ 238 const PetscScalar *PETSC_RESTRICT _p2 = p2; \ 239 const PetscScalar *PETSC_RESTRICT _p3 = p3; \ 240 const PetscScalar *PETSC_RESTRICT _p4 = p4; \ 241 PetscScalar *PETSC_RESTRICT _U = U; \ 242 PetscPragmaUseOMPKernels(parallel for) \ 243 for (PetscInt __i = 0; __i < _n; __i++) _U[__i] += _a1 * _p1[__i] + _a2 * _p2[__i] + _a3 * _p3[__i] + _a4 * _p4[__i]; \ 244 } while (0) 245 246 #endif 247