20#ifndef OPM_CUISTL_DEVICE_BLOCK_OPERATIONS_HPP
21#define OPM_CUISTL_DEVICE_BLOCK_OPERATIONS_HPP
24#include <cuda_runtime.h>
35template <
class T,
int blocksize>
36__device__ __forceinline__
void
37invBlockOutOfPlace(
const T* __restrict__ srcBlock, T* __restrict__ dstBlock)
40 dstBlock[0] = 1.0 / (srcBlock[0]);
41 }
else if (blocksize == 2) {
42 T detInv = 1.0 / (srcBlock[0] * srcBlock[3] - srcBlock[1] * srcBlock[2]);
43 dstBlock[0] = srcBlock[3] * detInv;
44 dstBlock[1] = -srcBlock[1] * detInv;
45 dstBlock[2] = -srcBlock[2] * detInv;
46 dstBlock[3] = srcBlock[0] * detInv;
47 }
else if (blocksize == 3) {
49 T t4 = srcBlock[0] * srcBlock[4];
50 T t6 = srcBlock[0] * srcBlock[5];
51 T t8 = srcBlock[1] * srcBlock[3];
52 T t10 = srcBlock[2] * srcBlock[3];
53 T t12 = srcBlock[1] * srcBlock[6];
54 T t14 = srcBlock[2] * srcBlock[6];
57 / (t4 * srcBlock[8] - t6 * srcBlock[7] - t8 * srcBlock[8] + t10 * srcBlock[7] + t12 * srcBlock[5]
60 dstBlock[0] = (srcBlock[4] * srcBlock[8] - srcBlock[5] * srcBlock[7]) * t17;
61 dstBlock[1] = -(srcBlock[1] * srcBlock[8] - srcBlock[2] * srcBlock[7]) * t17;
62 dstBlock[2] = (srcBlock[1] * srcBlock[5] - srcBlock[2] * srcBlock[4]) * t17;
63 dstBlock[3] = -(srcBlock[3] * srcBlock[8] - srcBlock[5] * srcBlock[6]) * t17;
64 dstBlock[4] = (srcBlock[0] * srcBlock[8] - t14) * t17;
65 dstBlock[5] = -(t6 - t10) * t17;
66 dstBlock[6] = (srcBlock[3] * srcBlock[7] - srcBlock[4] * srcBlock[6]) * t17;
67 dstBlock[7] = -(srcBlock[0] * srcBlock[7] - t12) * t17;
68 dstBlock[8] = (t4 - t8) * t17;
73template <
class T,
int blocksize>
74__device__ __forceinline__
void
75invBlockInPlace(T* __restrict__ block)
78 block[0] = 1.0 / (block[0]);
79 }
else if (blocksize == 2) {
80 T detInv = 1.0 / (block[0] * block[3] - block[1] * block[2]);
83 block[0] = block[3] * detInv;
84 block[1] = -block[1] * detInv;
85 block[2] = -block[2] * detInv;
86 block[3] = temp * detInv;
87 }
else if (blocksize == 3) {
88 T t4 = block[0] * block[4];
89 T t6 = block[0] * block[5];
90 T t8 = block[1] * block[3];
91 T t10 = block[2] * block[3];
92 T t12 = block[1] * block[6];
93 T t14 = block[2] * block[6];
95 T det = (t4 * block[8] - t6 * block[7] - t8 * block[8] + t10 * block[7] + t12 * block[5] - t14 * block[4]);
98 T matrix01 = block[1];
99 T matrix00 = block[0];
100 T matrix10 = block[3];
101 T matrix11 = block[4];
103 block[0] = (block[4] * block[8] - block[5] * block[7]) * t17;
104 block[1] = -(block[1] * block[8] - block[2] * block[7]) * t17;
105 block[2] = (matrix01 * block[5] - block[2] * block[4]) * t17;
106 block[3] = -(block[3] * block[8] - block[5] * block[6]) * t17;
107 block[4] = (matrix00 * block[8] - t14) * t17;
108 block[5] = -(t6 - t10) * t17;
109 block[6] = (matrix10 * block[7] - matrix11 * block[6]) * t17;
110 block[7] = -(matrix00 * block[7] - t12) * t17;
111 block[8] = (t4 - t8) * t17;
115enum class MVType { SET, PLUS, MINUS };
119template <
class T,
int blocksize, MVType OP>
120__device__ __forceinline__
void
121matrixVectorProductWithAction(
const T* A,
const T* b, T* c)
123 for (
int i = 0; i < blocksize; ++i) {
124 if (OP == MVType::SET) {
128 for (
int j = 0; j < blocksize; ++j) {
129 if (OP == MVType::SET || OP == MVType::PLUS) {
130 c[i] += A[i * blocksize + j] * b[j];
131 }
else if (OP == MVType::MINUS) {
132 c[i] -= A[i * blocksize + j] * b[j];
138template <
class T,
int blocksize>
139__device__ __forceinline__
void
140mv(
const T* a,
const T* b, T* c)
142 matrixVectorProductWithAction<T, blocksize, MVType::SET>(a, b, c);
145template <
class T,
int blocksize>
146__device__ __forceinline__
void
147umv(
const T* a,
const T* b, T* c)
149 matrixVectorProductWithAction<T, blocksize, MVType::PLUS>(a, b, c);
152template <
class T,
int blocksize>
153__device__ __forceinline__
void
154mmv(
const T* a,
const T* b, T* c)
156 matrixVectorProductWithAction<T, blocksize, MVType::MINUS>(a, b, c);
160template <
class T,
int blocksize>
161__device__ __forceinline__
void
162mmx2Subtraction(T* A, T* B, T* C, T* dst)
165 T tmp[blocksize * blocksize] = {0};
167 for (
int i = 0; i < blocksize; ++i) {
168 for (
int k = 0; k < blocksize; ++k) {
169 for (
int j = 0; j < blocksize; ++j) {
170 tmp[i * blocksize + j] += A[i * blocksize + k] * B[k * blocksize + j];
176 for (
int i = 0; i < blocksize; ++i) {
177 for (
int k = 0; k < blocksize; ++k) {
178 for (
int j = 0; j < blocksize; ++j) {
179 dst[i * blocksize + j] -= tmp[i * blocksize + k] * C[k * blocksize + j];
186template <
class T,
int blocksize>
187__device__ __forceinline__
void
188mmNoOverlap(T* A, T* B, T* C)
190 for (
int i = 0; i < blocksize; ++i) {
191 for (
int k = 0; k < blocksize; ++k) {
192 for (
int j = 0; j < blocksize; ++j) {
193 C[i * blocksize + j] += A[i * blocksize + k] * B[k * blocksize + j];
200template <
class T,
int blocksize>
201__device__ __forceinline__
void
202mmOverlap(T* A, T* B, T* C)
204 T tmp[blocksize * blocksize] = {0};
205 for (
int i = 0; i < blocksize; ++i) {
206 for (
int k = 0; k < blocksize; ++k) {
207 for (
int j = 0; j < blocksize; ++j) {
208 tmp[i * blocksize + j] += A[i * blocksize + k] * B[k * blocksize + j];
213 for (
int i = 0; i < blocksize; ++i) {
214 for (
int j = 0; j < blocksize; ++j) {
215 C[i * blocksize + j] = tmp[i * blocksize + j];
221template <
class T,
int blocksize>
222__device__ __forceinline__
void
223matrixSubtraction(T* A, T* B)
226 for (
int i = 0; i < blocksize; ++i) {
227 for (
int j = 0; j < blocksize; ++j) {
228 A[i * blocksize + j] -= B[i * blocksize + j];