1c4762a1bSJed Brown 2c4762a1bSJed Brown static char help[] = "Tests for bugs in A->offloadmask consistency for GPU matrices\n\n"; 3c4762a1bSJed Brown 4c4762a1bSJed Brown #include <petscmat.h> 5c4762a1bSJed Brown 6c4762a1bSJed Brown int main(int argc,char **args) 7c4762a1bSJed Brown { 8c4762a1bSJed Brown Mat A; 9c4762a1bSJed Brown PetscInt i,j,rstart,rend,m = 3; 10c4762a1bSJed Brown PetscScalar one = 1.0,zero = 0.0,negativeone = -1.0; 11c4762a1bSJed Brown PetscReal norm; 12c4762a1bSJed Brown Vec x,y; 13c4762a1bSJed Brown 14*9566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 15*9566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 16c4762a1bSJed Brown 17c4762a1bSJed Brown for (i=0; i<2; i++) { 18c4762a1bSJed Brown /* Create the matrix and set it to contain explicit zero entries on the diagonal. */ 19*9566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 20*9566063dSJacob Faibussowitsch PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*m,m*m)); 21*9566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(A)); 22*9566063dSJacob Faibussowitsch PetscCall(MatSetUp(A)); 23*9566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRange(A,&rstart,&rend)); 24*9566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(A,&x,&y)); 25*9566063dSJacob Faibussowitsch PetscCall(VecSet(x,one)); 26*9566063dSJacob Faibussowitsch PetscCall(VecSet(y,zero)); 27*9566063dSJacob Faibussowitsch PetscCall(MatDiagonalSet(A,y,INSERT_VALUES)); 28c4762a1bSJed Brown 29c4762a1bSJed Brown /* Now set A to be the identity using various approaches. 30c4762a1bSJed Brown * Note that there may be other approaches that should be added here. */ 31c4762a1bSJed Brown switch (i) { 32c4762a1bSJed Brown case 0: 33*9566063dSJacob Faibussowitsch PetscCall(MatDiagonalSet(A,x,INSERT_VALUES)); 34c4762a1bSJed Brown break; 35c4762a1bSJed Brown case 1: 36c4762a1bSJed Brown for (j=rstart; j<rend; j++) { 37*9566063dSJacob Faibussowitsch PetscCall(MatSetValue(A,j,j,one,INSERT_VALUES)); 38c4762a1bSJed Brown } 39*9566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 40*9566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 41c4762a1bSJed Brown break; 42c4762a1bSJed Brown case 2: 43c4762a1bSJed Brown for (j=rstart; j<rend; j++) { 44*9566063dSJacob Faibussowitsch PetscCall(MatSetValuesRow(A,j,&one)); 45c4762a1bSJed Brown } 46*9566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 47*9566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 48c4762a1bSJed Brown default: 49c4762a1bSJed Brown break; 50c4762a1bSJed Brown } 51c4762a1bSJed Brown 52c4762a1bSJed Brown /* Compute y <- A*x and verify that the difference between y and x is negligible, as it should be since A is the identity. */ 53*9566063dSJacob Faibussowitsch PetscCall(MatMult(A,x,y)); 54*9566063dSJacob Faibussowitsch PetscCall(VecAXPY(y,negativeone,x)); 55*9566063dSJacob Faibussowitsch PetscCall(VecNorm(y,NORM_2,&norm)); 56c4762a1bSJed Brown if (norm > PETSC_SQRT_MACHINE_EPSILON) { 57*9566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test %" PetscInt_FMT ": Norm of error is %g, but should be near 0.\n",i,(double)norm)); 58c4762a1bSJed Brown } 59c4762a1bSJed Brown 60*9566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A)); 61*9566063dSJacob Faibussowitsch PetscCall(VecDestroy(&x)); 62*9566063dSJacob Faibussowitsch PetscCall(VecDestroy(&y)); 63c4762a1bSJed Brown } 64c4762a1bSJed Brown 65*9566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 66b122ec5aSJacob Faibussowitsch return 0; 67c4762a1bSJed Brown } 68c4762a1bSJed Brown 69c4762a1bSJed Brown /*TEST 70c4762a1bSJed Brown 71c4762a1bSJed Brown test: 72c4762a1bSJed Brown suffix: aijviennacl_1 73c4762a1bSJed Brown nsize: 1 74c4762a1bSJed Brown args: -mat_type aijviennacl 75c4762a1bSJed Brown requires: viennacl 76c4762a1bSJed Brown 77c4762a1bSJed Brown test: 78c4762a1bSJed Brown suffix: aijviennacl_2 79c4762a1bSJed Brown nsize: 2 80c4762a1bSJed Brown args: -mat_type aijviennacl 81c4762a1bSJed Brown requires: viennacl 82c4762a1bSJed Brown 83c4762a1bSJed Brown test: 84c4762a1bSJed Brown suffix: aijcusparse_1 85c4762a1bSJed Brown nsize: 1 86c4762a1bSJed Brown args: -mat_type aijcusparse 87c4762a1bSJed Brown requires: cuda 88c4762a1bSJed Brown 89c4762a1bSJed Brown test: 90c4762a1bSJed Brown suffix: aijcusparse_2 91c4762a1bSJed Brown nsize: 2 92c4762a1bSJed Brown args: -mat_type aijcusparse 93c4762a1bSJed Brown requires: cuda 94c4762a1bSJed Brown TEST*/ 95