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