Actual source code: ex69.c

  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: {
  7:   Mat A;

 10:   MatShellGetContext(S, &A);
 11:   MatMult(A, x, y);
 12:   return 0;
 13: }

 15: static PetscBool test_cusparse_transgen = PETSC_FALSE;

 17: static PetscErrorCode MatMultTranspose_S(Mat S, Vec x, Vec y)
 18: {
 19:   Mat A;

 22:   MatShellGetContext(S, &A);
 23:   MatMultTranspose(A, x, y);

 25:   /* alternate transgen true and false to test code logic */
 26:   MatSetOption(A, MAT_FORM_EXPLICIT_TRANSPOSE, test_cusparse_transgen);
 27:   test_cusparse_transgen = (PetscBool)!test_cusparse_transgen;
 28:   return 0;
 29: }

 31: int main(int argc, char **argv)
 32: {
 33:   Mat          A, B, C, S;
 34:   Vec          t, v;
 35:   PetscScalar *vv, *aa;
 36:   PetscInt     n = 30, k = 6, l = 0, i, Istart, Iend, nloc, bs, test = 1;
 37:   PetscBool    flg, reset, use_shell = PETSC_FALSE;
 38:   VecType      vtype;

 41:   PetscInitialize(&argc, &argv, (char *)0, help);
 42:   PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
 43:   PetscOptionsGetInt(NULL, NULL, "-k", &k, NULL);
 44:   PetscOptionsGetInt(NULL, NULL, "-l", &l, NULL);
 45:   PetscOptionsGetInt(NULL, NULL, "-test", &test, NULL);
 46:   PetscOptionsGetBool(NULL, NULL, "-use_shell", &use_shell, NULL);

 51:   /* sparse matrix */
 52:   MatCreate(PETSC_COMM_WORLD, &A);
 53:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n);
 54:   MatSetType(A, MATAIJCUSPARSE);
 55:   MatSetOptionsPrefix(A, "A_");
 56:   MatSetFromOptions(A);
 57:   MatSetUp(A);

 59:   /* test special case for SeqAIJCUSPARSE to generate explicit transpose (not default) */
 60:   MatSetOption(A, MAT_FORM_EXPLICIT_TRANSPOSE, PETSC_TRUE);

 62:   MatGetOwnershipRange(A, &Istart, &Iend);
 63:   for (i = Istart; i < Iend; i++) {
 64:     if (i > 0) MatSetValue(A, i, i - 1, -1.0, INSERT_VALUES);
 65:     if (i < n - 1) MatSetValue(A, i, i + 1, -1.0, INSERT_VALUES);
 66:     MatSetValue(A, i, i, 2.0, INSERT_VALUES);
 67:   }
 68:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 69:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

 71:   /* template vector */
 72:   MatCreateVecs(A, NULL, &t);
 73:   VecGetType(t, &vtype);

 75:   /* long vector, contains the stacked columns of an nxk dense matrix */
 76:   VecGetLocalSize(t, &nloc);
 77:   VecGetBlockSize(t, &bs);
 78:   VecCreate(PetscObjectComm((PetscObject)t), &v);
 79:   VecSetType(v, vtype);
 80:   VecSetSizes(v, k * nloc, k * n);
 81:   VecSetBlockSize(v, bs);
 82:   VecSetRandom(v, NULL);

 84:   /* dense matrix that contains the columns of v */
 85:   VecCUDAGetArray(v, &vv);

 87:   /* test few cases for MatDenseCUDA handling pointers */
 88:   switch (test) {
 89:   case 1:
 90:     MatCreateDenseCUDA(PetscObjectComm((PetscObject)v), nloc, PETSC_DECIDE, n, k - l, vv, &B); /* pass a pointer to avoid allocation of storage */
 91:     MatDenseCUDAReplaceArray(B, NULL);                                                         /* replace with a null pointer, the value after BVRestoreMat */
 92:     MatDenseCUDAPlaceArray(B, vv + l * nloc);                                                  /* set the actual pointer */
 93:     reset = PETSC_TRUE;
 94:     break;
 95:   case 2:
 96:     MatCreateDenseCUDA(PetscObjectComm((PetscObject)v), nloc, PETSC_DECIDE, n, k - l, NULL, &B);
 97:     MatDenseCUDAPlaceArray(B, vv + l * nloc); /* set the actual pointer */
 98:     reset = PETSC_TRUE;
 99:     break;
100:   default:
101:     MatCreateDenseCUDA(PetscObjectComm((PetscObject)v), nloc, PETSC_DECIDE, n, k - l, vv + l * nloc, &B);
102:     reset = PETSC_FALSE;
103:     break;
104:   }
105:   VecCUDARestoreArray(v, &vv);

107:   /* Test MatMatMult */
108:   if (use_shell) {
109:     /* we could have called the general convertor below, but we explicit set the operations
110:        ourselves to test MatProductSymbolic_X_Dense, MatProductNumeric_X_Dense code */
111:     /* MatConvert(A,MATSHELL,MAT_INITIAL_MATRIX,&S); */
112:     MatCreateShell(PetscObjectComm((PetscObject)v), nloc, nloc, n, n, A, &S);
113:     MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MatMult_S);
114:     MatShellSetOperation(S, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_S);
115:     MatShellSetVecType(S, vtype);
116:   } else {
117:     PetscObjectReference((PetscObject)A);
118:     S = A;
119:   }

121:   MatCreateDenseCUDA(PetscObjectComm((PetscObject)v), nloc, PETSC_DECIDE, n, k - l, NULL, &C);

123:   /* test MatMatMult */
124:   MatProductCreateWithMat(S, B, NULL, C);
125:   MatProductSetType(C, MATPRODUCT_AB);
126:   MatProductSetFromOptions(C);
127:   MatProductSymbolic(C);
128:   MatProductNumeric(C);
129:   MatMatMultEqual(S, B, C, 10, &flg);
130:   if (!flg) PetscPrintf(PETSC_COMM_WORLD, "Error MatMatMult\n");

132:   /* test MatTransposeMatMult */
133:   MatProductCreateWithMat(S, B, NULL, C);
134:   MatProductSetType(C, MATPRODUCT_AtB);
135:   MatProductSetFromOptions(C);
136:   MatProductSymbolic(C);
137:   MatProductNumeric(C);
138:   MatTransposeMatMultEqual(S, B, C, 10, &flg);
139:   if (!flg) PetscPrintf(PETSC_COMM_WORLD, "Error MatTransposeMatMult\n");

141:   MatDestroy(&C);
142:   MatDestroy(&S);

144:   /* finished using B */
145:   MatDenseCUDAGetArray(B, &aa);
147:   MatDenseCUDARestoreArray(B, &aa);
148:   if (reset) MatDenseCUDAResetArray(B);
149:   VecCUDARestoreArray(v, &vv);

151:   if (test == 1) {
152:     MatDenseCUDAGetArray(B, &aa);
154:     MatDenseCUDARestoreArray(B, &aa);
155:   }

157:   /* free work space */
158:   MatDestroy(&B);
159:   MatDestroy(&A);
160:   VecDestroy(&t);
161:   VecDestroy(&v);
162:   PetscFinalize();
163:   return 0;
164: }

166: /*TEST

168:   build:
169:     requires: cuda

171:   test:
172:     requires: cuda
173:     suffix: 1
174:     nsize: {{1 2}}
175:     args: -A_mat_type {{aij aijcusparse}} -test {{0 1 2}} -k 6 -l {{0 5}} -use_shell {{0 1}}

177: TEST*/