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