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 *, void *, void *, void *, void *, const void *, const void *, const void *, const void *, PetscInt *); 28 PETSC_EXTERN void fortranmaxpy3_(void *, void *, void *, void *, const void *, const void *, const void *, PetscInt *); 29 PETSC_EXTERN void fortranmaxpy2_(void *, void *, void *, const void *, const void *, 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 #elif defined(PETSC_USE_FOR_KERNELS) 198 199 #define PetscKernelAXPY(U, a1, p1, n) \ 200 { \ 201 PetscInt __i; \ 202 PetscScalar __s1, __s2; \ 203 for (__i = 0; __i < n - 1; __i += 2) { \ 204 __s1 = a1 * p1[__i]; \ 205 __s2 = a1 * p1[__i + 1]; \ 206 __s1 += U[__i]; \ 207 __s2 += U[__i + 1]; \ 208 U[__i] = __s1; \ 209 U[__i + 1] = __s2; \ 210 } \ 211 if (n & 0x1) U[__i] += a1 * p1[__i]; \ 212 } 213 #define PetscKernelAXPY2(U, a1, a2, p1, p2, n) \ 214 { \ 215 PetscInt __i; \ 216 for (__i = 0; __i < n; __i++) U[__i] += a1 * p1[__i] + a2 * p2[__i]; \ 217 } 218 #define PetscKernelAXPY3(U, a1, a2, a3, p1, p2, p3, n) \ 219 { \ 220 PetscInt __i; \ 221 for (__i = 0; __i < n; __i++) U[__i] += a1 * p1[__i] + a2 * p2[__i] + a3 * p3[__i]; \ 222 } 223 #define PetscKernelAXPY4(U, a1, a2, a3, a4, p1, p2, p3, p4, n) \ 224 { \ 225 PetscInt __i; \ 226 for (__i = 0; __i < n; __i++) U[__i] += a1 * p1[__i] + a2 * p2[__i] + a3 * p3[__i] + a4 * p4[__i]; \ 227 } 228 229 #else 230 231 #define PetscKernelAXPY(U, a1, p1, n) \ 232 { \ 233 PetscInt __i; \ 234 PetscScalar _a1 = a1; \ 235 for (__i = 0; __i < n; __i++) U[__i] += _a1 * p1[__i]; \ 236 } 237 #define PetscKernelAXPY2(U, a1, a2, p1, p2, n) \ 238 { \ 239 PetscInt __i; \ 240 for (__i = 0; __i < n; __i++) U[__i] += a1 * p1[__i] + a2 * p2[__i]; \ 241 } 242 #define PetscKernelAXPY3(U, a1, a2, a3, p1, p2, p3, n) \ 243 { \ 244 PetscInt __i; \ 245 for (__i = 0; __i < n; __i++) U[__i] += a1 * p1[__i] + a2 * p2[__i] + a3 * p3[__i]; \ 246 } 247 #define PetscKernelAXPY4(U, a1, a2, a3, a4, p1, p2, p3, p4, n) \ 248 { \ 249 PetscInt __i; \ 250 for (__i = 0; __i < n; __i++) U[__i] += a1 * p1[__i] + a2 * p2[__i] + a3 * p3[__i] + a4 * p4[__i]; \ 251 } 252 253 #endif 254 255 #endif // PETSC_KERNELS_PETSCAXPY_H 256