Actual source code: ex324.c
1: static const char help[] = "Test CUDA/HIP matrix storage and value updates.\n";
3: #include <petscmat.h>
4: #include <petscdevice.h>
5: #if PetscDefined(HAVE_CUDA)
6: #include <petscdevice_cuda.h>
7: #endif
8: #if PetscDefined(HAVE_HIP)
9: #include <petscdevice_hip.h>
10: #endif
12: static PetscErrorCode TestSELL(PetscDeviceType type, PetscDeviceContext dctx)
13: {
14: Mat A, B;
15: Vec x, y, z;
16: const PetscScalar *a;
17: const PetscScalar expected[] = {3, 3, 5, 9, 5, 8, 10};
18: #if PetscDefined(USE_LOG) && PetscDefined(HAVE_DEVICE)
19: PetscLogEvent event;
20: PetscEventPerfInfo before, after;
21: MatInfo info;
22: #endif
24: PetscFunctionBeginUser;
25: #if PetscDefined(USE_LOG) && PetscDefined(HAVE_DEVICE)
26: PetscCall(PetscLogDefaultBegin());
27: PetscCall(PetscLogEventRegister("SELLUpdate", MAT_CLASSID, &event));
28: #endif
29: PetscCall(PetscDeviceContextSynchronize(dctx));
30: PetscCall(MatCreate(PETSC_COMM_SELF, &A));
31: PetscCall(MatSetSizes(A, 32, 32, 32, 32));
32: PetscCall(MatSetType(A, MATSEQSELL));
33: PetscCall(MatSetUp(A));
34: for (PetscInt i = 0; i < 32; ++i) PetscCall(MatSetValue(A, i, i, 2, INSERT_VALUES));
35: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
36: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
37: PetscCall(MatConvert(A, type == PETSC_DEVICE_CUDA ? MATSEQSELLCUDA : MATSEQSELLHIP, MAT_INPLACE_MATRIX, &A));
38: PetscCall(MatSetFromOptions(A));
39: PetscCall(MatCreateVecs(A, &x, &y));
40: PetscCall(VecDuplicate(y, &z));
41: PetscCall(VecSet(x, 1));
42: PetscCall(VecSet(y, 1));
43: // Destroying a duplicate before its first upload must release its backend data.
44: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
45: PetscCall(MatDestroy(&B));
46: // MatMultAdd() must also work before any multiply has initialized the device arrays.
47: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
48: PetscCall(MatSetFromOptions(B));
49: for (size_t stage = 0; stage < PETSC_STATIC_ARRAY_LENGTH(expected); ++stage) {
50: Mat C = stage == 0 ? A : B;
52: #if PetscDefined(USE_LOG) && PetscDefined(HAVE_DEVICE)
53: PetscCall(PetscLogEventGetPerfInfo(PETSC_DETERMINE, event, &before));
54: PetscCall(PetscLogEventBegin(event, C, 0, 0, 0));
55: #endif
56: if (stage == 2) {
57: PetscCall(MatScale(A, 2));
58: PetscCall(MatCopy(A, B, SAME_NONZERO_PATTERN));
59: PetscCall(MatScale(A, 0.5));
60: } else if (stage == 3) {
61: PetscCall(MatSetOption(C, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE));
62: PetscCall(MatStoreValues(C));
63: PetscCall(MatScale(C, 2));
64: } else if (stage == 4) PetscCall(MatRetrieveValues(C));
65: else if (stage == 5) {
66: PetscCall(MatSetOption(C, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE));
67: PetscCall(MatSetOption(C, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
68: for (PetscInt i = 0; i < 32; ++i) PetscCall(MatSetValue(C, i, (i + 1) % 32, 3, INSERT_VALUES));
69: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
70: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
71: } else if (stage == 6) {
72: for (PetscInt i = 0; i < 32; ++i) PetscCall(MatSetValue(C, i, i, 6, INSERT_VALUES));
73: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
74: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
75: }
76: PetscCall(MatMultAdd(C, x, y, z));
77: PetscCall(PetscDeviceContextSynchronize(dctx));
78: #if PetscDefined(USE_LOG) && PetscDefined(HAVE_DEVICE)
79: PetscCall(PetscLogEventEnd(event, C, 0, 0, 0));
80: if (stage == 2 || stage == 3 || stage == 4 || stage == 6) {
81: PetscCall(PetscLogEventGetPerfInfo(PETSC_DETERMINE, event, &after));
82: PetscCall(MatGetInfo(C, MAT_LOCAL, &info));
83: PetscCheck(after.CpuToGpuSize - before.CpuToGpuSize == info.nz_used * sizeof(PetscScalar), PETSC_COMM_SELF, PETSC_ERR_PLIB, "SELL values-only update stage %zu uploaded %g bytes instead of %g", stage, after.CpuToGpuSize - before.CpuToGpuSize, info.nz_used * sizeof(PetscScalar));
84: }
85: #endif
86: PetscCall(VecGetArrayRead(z, &a));
87: for (PetscInt i = 0; i < 32; ++i) PetscCheck(a[i] == expected[stage], PETSC_COMM_SELF, PETSC_ERR_PLIB, "SELL MatMultAdd() stage %zu returned an incorrect entry %" PetscInt_FMT, stage, i);
88: PetscCall(VecRestoreArrayRead(z, &a));
89: }
90: PetscCall(MatDestroy(&B));
91: PetscCall(VecDestroy(&z));
92: PetscCall(VecDestroy(&y));
93: PetscCall(VecDestroy(&x));
94: PetscCall(MatDestroy(&A));
95: PetscFunctionReturn(PETSC_SUCCESS);
96: }
98: #if !PetscDefined(USE_COMPLEX)
99: static PetscErrorCode TestSELLChunks(PetscDeviceType type)
100: {
101: Mat A, B, C;
102: Vec x, y, z, x_ref, y_ref;
103: const PetscScalar *a, *a_ref;
104: PetscReal avgwidth;
105: PetscInt maxwidth;
107: PetscFunctionBeginUser;
108: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, 512, 512, 0, NULL, &A));
109: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
110: for (PetscInt i = 0; i < 512; ++i) PetscCall(MatSetValue(A, i, i, 2, INSERT_VALUES));
111: // One wide slice makes the default kernel selection use the chunk map.
112: for (PetscInt j = 1; j < 256; ++j) PetscCall(MatSetValue(A, 0, j, 1, INSERT_VALUES));
113: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
114: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
115: PetscCall(MatConvert(A, type == PETSC_DEVICE_CUDA ? MATSEQSELLCUDA : MATSEQSELLHIP, MAT_INITIAL_MATRIX, &B));
116: PetscCall(MatSeqSELLGetMaxSliceWidth(B, &maxwidth));
117: PetscCall(MatSeqSELLGetAvgSliceWidth(B, &avgwidth));
118: PetscCheck(maxwidth / avgwidth > 12 && maxwidth / avgwidth / 512 > 0.001, PETSC_COMM_SELF, PETSC_ERR_PLIB, "SELL regression must select the chunked kernel");
119: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &C));
120: PetscCall(MatDestroy(&B));
121: PetscCall(MatSetOption(C, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
122: PetscCall(MatCreateVecs(C, &x, &y));
123: PetscCall(VecDuplicate(y, &z));
124: PetscCall(MatCreateVecs(A, &x_ref, &y_ref));
125: PetscCall(VecSet(x, 1));
126: PetscCall(VecSet(x_ref, 1));
127: for (PetscInt stage = 0; stage < 3; ++stage) {
128: if (stage) {
129: PetscCall(MatSetValue(A, 16, 17, stage, INSERT_VALUES));
130: PetscCall(MatSetValue(C, 16, 17, stage, INSERT_VALUES));
131: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
132: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
133: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
134: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
135: }
136: PetscCall(MatMult(A, x_ref, y_ref));
137: PetscCall(MatMult(C, x, y));
138: PetscCall(MatMultAdd(C, x, y, z));
139: PetscCall(VecGetArrayRead(y_ref, &a_ref));
140: for (PetscInt add = 0; add < 2; ++add) {
141: PetscCall(VecGetArrayRead(add ? z : y, &a));
142: for (PetscInt i = 0; i < 512; ++i) PetscCheck(a[i] == (add + 1) * a_ref[i], PETSC_COMM_SELF, PETSC_ERR_PLIB, "SELL chunked multiply stage %" PetscInt_FMT " add %" PetscInt_FMT " returned an incorrect entry %" PetscInt_FMT, stage, add, i);
143: PetscCall(VecRestoreArrayRead(add ? z : y, &a));
144: }
145: PetscCall(VecRestoreArrayRead(y_ref, &a_ref));
146: }
147: PetscCall(VecDestroy(&x_ref));
148: PetscCall(VecDestroy(&y_ref));
149: PetscCall(VecDestroy(&x));
150: PetscCall(VecDestroy(&y));
151: PetscCall(VecDestroy(&z));
152: PetscCall(MatDestroy(&A));
153: PetscCall(MatDestroy(&C));
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
156: #endif
158: static PetscErrorCode TestSparseSolves(PetscDeviceType type)
159: {
160: Mat A, F;
161: Vec b, x;
162: IS perm;
163: MatFactorInfo info;
164: const PetscScalar *a;
166: PetscFunctionBeginUser;
167: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, 4, 4, 1, NULL, &A));
168: for (PetscInt i = 0; i < 4; ++i) PetscCall(MatSetValue(A, i, i, 2, INSERT_VALUES));
169: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
170: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
171: PetscCall(MatConvert(A, type == PETSC_DEVICE_CUDA ? MATSEQAIJCUSPARSE : MATSEQAIJHIPSPARSE, MAT_INPLACE_MATRIX, &A));
172: PetscCall(MatCreateVecs(A, &b, &x));
173: PetscCall(VecSet(b, 2));
174: PetscCall(ISCreateStride(PETSC_COMM_SELF, 4, 0, 1, &perm));
175: PetscCall(MatFactorInfoInitialize(&info));
176: for (PetscInt kind = 0; kind < 2; ++kind) {
177: PetscCall(MatGetFactor(A, type == PETSC_DEVICE_CUDA ? MATSOLVERCUSPARSE : MATSOLVERHIPSPARSE, kind ? MAT_FACTOR_CHOLESKY : MAT_FACTOR_LU, &F));
178: if (kind) {
179: PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
180: PetscCall(MatCholeskyFactorNumeric(F, A, &info));
181: } else {
182: PetscCall(MatLUFactorSymbolic(F, A, perm, perm, &info));
183: PetscCall(MatLUFactorNumeric(F, A, &info));
184: }
185: for (PetscInt transpose = 0; transpose < 2; ++transpose) {
186: if (transpose) PetscCall(MatSolveTranspose(F, b, x));
187: else PetscCall(MatSolve(F, b, x));
188: PetscCall(VecGetArrayRead(x, &a));
189: for (PetscInt i = 0; i < 4; ++i) PetscCheck(PetscAbsScalar(a[i] - 1) <= 10 * PETSC_MACHINE_EPSILON, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Incorrect sparse solve at entry %" PetscInt_FMT, i);
190: PetscCall(VecRestoreArrayRead(x, &a));
191: }
192: PetscCall(MatDestroy(&F));
193: }
194: PetscCall(ISDestroy(&perm));
195: PetscCall(VecDestroy(&x));
196: PetscCall(VecDestroy(&b));
197: PetscCall(MatDestroy(&A));
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: static PetscErrorCode TestTransposeUpdate(PetscDeviceType type)
202: {
203: Mat A;
204: Vec x, y;
205: PetscScalar *values;
206: const PetscScalar *a;
207: const PetscScalar ones[] = {1, 1, 1, 1};
209: PetscFunctionBeginUser;
210: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, 4, 4, 1, NULL, &A));
211: for (PetscInt i = 0; i < 4; ++i) PetscCall(MatSetValue(A, i, i, 2, INSERT_VALUES));
212: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
213: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
214: PetscCall(MatConvert(A, type == PETSC_DEVICE_CUDA ? MATSEQAIJCUSPARSE : MATSEQAIJHIPSPARSE, MAT_INPLACE_MATRIX, &A));
215: PetscCall(MatSetOption(A, MAT_FORM_EXPLICIT_TRANSPOSE, PETSC_TRUE));
216: PetscCall(MatCreateVecs(A, &x, &y));
217: PetscCall(VecSet(x, 2));
218: PetscCall(MatMultTranspose(A, x, y));
219: #if PetscDefined(HAVE_CUDA)
220: if (type == PETSC_DEVICE_CUDA) {
221: PetscCall(MatSeqAIJCUSPARSEGetArray(A, &values));
222: PetscCallCUDA(cudaMemcpy(values, ones, sizeof(ones), cudaMemcpyHostToDevice));
223: PetscCall(MatSeqAIJCUSPARSERestoreArray(A, &values));
224: }
225: #endif
226: #if PetscDefined(HAVE_HIP)
227: if (type == PETSC_DEVICE_HIP) {
228: PetscCall(MatSeqAIJHIPSPARSEGetArray(A, &values));
229: PetscCallHIP(hipMemcpy(values, ones, sizeof(ones), hipMemcpyHostToDevice));
230: PetscCall(MatSeqAIJHIPSPARSERestoreArray(A, &values));
231: }
232: #endif
233: PetscCall(MatMultTranspose(A, x, y));
234: PetscCall(VecGetArrayRead(y, &a));
235: for (PetscInt i = 0; i < 4; ++i) PetscCheck(a[i] == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Transpose multiply used stale values at entry %" PetscInt_FMT, i);
236: PetscCall(VecRestoreArrayRead(y, &a));
237: PetscCall(VecDestroy(&y));
238: PetscCall(VecDestroy(&x));
239: PetscCall(MatDestroy(&A));
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: int main(int argc, char **argv)
244: {
245: PetscDeviceContext dctx;
246: PetscDeviceType type;
248: PetscFunctionBeginUser;
249: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
250: PetscCall(PetscDeviceContextGetCurrentContext(&dctx));
251: PetscCall(PetscDeviceContextGetDeviceType(dctx, &type));
252: PetscCall(TestSELL(type, dctx));
253: #if !PetscDefined(USE_COMPLEX)
254: PetscCall(TestSELLChunks(type));
255: #endif
256: PetscCall(TestSparseSolves(type));
257: PetscCall(TestTransposeUpdate(type));
258: PetscCall(PetscFinalize());
259: return 0;
260: }
262: /*TEST
264: testset:
265: nsize: {{1 2}}
266: output_file: output/empty.out
267: test:
268: suffix: cuda
269: requires: cuda
270: args: -default_device_type cuda -mat_sell_spmv_cuda_kernel 1
271: test:
272: suffix: hip
273: requires: hip
274: args: -default_device_type hip -mat_sell_spmv_hip_kernel 1
276: TEST*/