Actual source code: ex69.c
petsc-3.14.6 2021-03-30
1: static char help[] = "Tests MatCreateDenseCUDA(), MatDenseCUDAPlaceArray(), MatDenseCUDAReplaceArray(), MatDenseCUDAResetArray()\n";
3: #include <petscmat.h>
5: static PetscErrorCode MatMult_S(Mat S,Vec x,Vec y)
6: {
8: Mat A;
11: MatShellGetContext(S,(void**)&A);
12: MatMult(A,x,y);
13: return(0);
14: }
16: static PetscBool test_cusparse_transgen = PETSC_FALSE;
18: static PetscErrorCode MatMultTranspose_S(Mat S,Vec x,Vec y)
19: {
21: Mat A;
24: MatShellGetContext(S,(void**)&A);
25: MatMultTranspose(A,x,y);
27: /* alternate transgen true and false to test code logic */
28: MatSeqAIJCUSPARSESetGenerateTranspose(A,test_cusparse_transgen);
29: test_cusparse_transgen = (PetscBool)!test_cusparse_transgen;
30: return(0);
31: }
33: int main(int argc,char **argv)
34: {
35: Mat A,B,C,S;
36: Vec t,v;
37: PetscScalar *vv,*aa;
38: PetscInt n=30,k=6,l=0,i,Istart,Iend,nloc,bs,test=1;
39: PetscBool flg,reset,use_shell = PETSC_FALSE;
41: VecType vtype;
43: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
44: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
45: PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL);
46: PetscOptionsGetInt(NULL,NULL,"-l",&l,NULL);
47: PetscOptionsGetInt(NULL,NULL,"-test",&test,NULL);
48: PetscOptionsGetBool(NULL,NULL,"-use_shell",&use_shell,NULL);
49: if (k < 0) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_USER,"k %D must be positive",k);
50: if (l < 0) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_USER,"l %D must be positive",l);
51: if (l > k) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_USER,"l %D must be smaller or equal than k %D\n",l,k);
53: /* sparse matrix */
54: MatCreate(PETSC_COMM_WORLD,&A);
55: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
56: MatSetType(A,MATAIJCUSPARSE);
57: MatSetOptionsPrefix(A,"A_");
58: MatSetFromOptions(A);
59: MatSetUp(A);
61: /* test special case for SeqAIJCUSPARSE to generate explicit transpose (not default) */
62: MatSeqAIJCUSPARSESetGenerateTranspose(A,PETSC_TRUE);
64: MatGetOwnershipRange(A,&Istart,&Iend);
65: for (i=Istart;i<Iend;i++) {
66: if (i>0) { MatSetValue(A,i,i-1,-1.0,INSERT_VALUES); }
67: if (i<n-1) { MatSetValue(A,i,i+1,-1.0,INSERT_VALUES); }
68: MatSetValue(A,i,i,2.0,INSERT_VALUES);
69: }
70: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
71: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
73: /* template vector */
74: MatCreateVecs(A,NULL,&t);
75: VecGetType(t,&vtype);
77: /* long vector, contains the stacked columns of an nxk dense matrix */
78: VecGetLocalSize(t,&nloc);
79: VecGetBlockSize(t,&bs);
80: VecCreate(PetscObjectComm((PetscObject)t),&v);
81: VecSetType(v,vtype);
82: VecSetSizes(v,k*nloc,k*n);
83: VecSetBlockSize(v,bs);
84: VecSetRandom(v,NULL);
86: /* dense matrix that contains the columns of v */
87: VecCUDAGetArray(v,&vv);
89: /* test few cases for MatDenseCUDA handling pointers */
90: switch (test) {
91: case 1:
92: MatCreateDenseCUDA(PetscObjectComm((PetscObject)v),nloc,PETSC_DECIDE,n,k-l,vv,&B); /* pass a pointer to avoid allocation of storage */
93: MatDenseCUDAReplaceArray(B,NULL); /* replace with a null pointer, the value after BVRestoreMat */
94: MatDenseCUDAPlaceArray(B,vv+l*nloc); /* set the actual pointer */
95: reset = PETSC_TRUE;
96: break;
97: case 2:
98: MatCreateDenseCUDA(PetscObjectComm((PetscObject)v),nloc,PETSC_DECIDE,n,k-l,NULL,&B);
99: MatDenseCUDAPlaceArray(B,vv+l*nloc); /* set the actual pointer */
100: reset = PETSC_TRUE;
101: break;
102: default:
103: MatCreateDenseCUDA(PetscObjectComm((PetscObject)v),nloc,PETSC_DECIDE,n,k-l,vv+l*nloc,&B);
104: reset = PETSC_FALSE;
105: break;
106: }
107: VecCUDARestoreArray(v,&vv);
109: /* Test MatMatMult */
110: if (use_shell) {
111: /* we could have called the general convertor below, but we explicit set the operations
112: ourselves to test MatProductSymbolic_X_Dense, MatProductNumeric_X_Dense code */
113: /* MatConvert(A,MATSHELL,MAT_INITIAL_MATRIX,&S); */
114: MatCreateShell(PetscObjectComm((PetscObject)v),nloc,nloc,n,n,A,&S);
115: MatShellSetOperation(S,MATOP_MULT,(void(*)(void))MatMult_S);
116: MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_S);
117: MatShellSetVecType(S,vtype);
118: } else {
119: PetscObjectReference((PetscObject)A);
120: S = A;
121: }
123: MatCreateDenseCUDA(PetscObjectComm((PetscObject)v),nloc,PETSC_DECIDE,n,k-l,NULL,&C);
124: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
125: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
126: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
127: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
129: /* test MatMatMult */
130: MatProductCreateWithMat(S,B,NULL,C);
131: MatProductSetType(C,MATPRODUCT_AB);
132: MatProductSetFromOptions(C);
133: MatProductSymbolic(C);
134: MatProductNumeric(C);
135: MatMatMultEqual(S,B,C,10,&flg);
136: if (!flg) { PetscPrintf(PETSC_COMM_WORLD,"Error MatMatMult\n"); }
138: /* test MatTransposeMatMult */
139: MatProductCreateWithMat(S,B,NULL,C);
140: MatProductSetType(C,MATPRODUCT_AtB);
141: MatProductSetFromOptions(C);
142: MatProductSymbolic(C);
143: MatProductNumeric(C);
144: MatTransposeMatMultEqual(S,B,C,10,&flg);
145: if (!flg) { PetscPrintf(PETSC_COMM_WORLD,"Error MatTransposeMatMult\n"); }
147: MatDestroy(&C);
148: MatDestroy(&S);
150: /* finished using B */
151: MatDenseCUDAGetArray(B,&aa);
152: if (vv != aa-l*nloc) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong array");
153: MatDenseCUDARestoreArray(B,&aa);
154: if (reset) {
155: MatDenseCUDAResetArray(B);
156: }
157: VecCUDARestoreArray(v,&vv);
159: if (test == 1) {
160: MatDenseCUDAGetArray(B,&aa);
161: if (aa) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expected a null pointer");
162: MatDenseCUDARestoreArray(B,&aa);
163: }
165: /* free work space */
166: MatDestroy(&B);
167: MatDestroy(&A);
168: VecDestroy(&t);
169: VecDestroy(&v);
170: PetscFinalize();
171: return ierr;
172: }
174: /*TEST
176: build:
177: requires: cuda
179: test:
180: requires: cuda
181: suffix: 1
182: nsize: {{1 2}}
183: args: -A_mat_type {{aij aijcusparse}} -test {{0 1 2}} -k 6 -l {{0 5}} -use_shell {{0 1}}
185: TEST*/