1 2 static char help[] = "Test interface of Elemental. \n\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **args) 7 { 8 Mat C,Caij; 9 PetscInt i,j,m = 5,n,nrows,ncols; 10 const PetscInt *rows,*cols; 11 IS isrows,iscols; 12 PetscBool flg,Test_MatMatMult=PETSC_FALSE,mats_view=PETSC_FALSE; 13 PetscScalar *v; 14 PetscMPIInt rank,size; 15 16 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 17 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 18 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 19 PetscCall(PetscOptionsHasName(NULL,NULL,"-mats_view",&mats_view)); 20 21 /* Get local block or element size*/ 22 PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 23 n = m; 24 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 25 26 PetscCall(MatCreate(PETSC_COMM_WORLD,&C)); 27 PetscCall(MatSetSizes(C,m,n,PETSC_DECIDE,PETSC_DECIDE)); 28 PetscCall(MatSetType(C,MATELEMENTAL)); 29 PetscCall(MatSetFromOptions(C)); 30 PetscCall(MatSetUp(C)); 31 32 PetscCall(PetscOptionsHasName(NULL,NULL,"-row_oriented",&flg)); 33 if (flg) PetscCall(MatSetOption(C,MAT_ROW_ORIENTED,PETSC_TRUE)); 34 PetscCall(MatGetOwnershipIS(C,&isrows,&iscols)); 35 PetscCall(PetscOptionsHasName(NULL,NULL,"-Cexp_view_ownership",&flg)); 36 if (flg) { /* View ownership of explicit C */ 37 IS tmp; 38 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Ownership of explicit C:\n")); 39 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Row index set:\n")); 40 PetscCall(ISOnComm(isrows,PETSC_COMM_WORLD,PETSC_USE_POINTER,&tmp)); 41 PetscCall(ISView(tmp,PETSC_VIEWER_STDOUT_WORLD)); 42 PetscCall(ISDestroy(&tmp)); 43 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Column index set:\n")); 44 PetscCall(ISOnComm(iscols,PETSC_COMM_WORLD,PETSC_USE_POINTER,&tmp)); 45 PetscCall(ISView(tmp,PETSC_VIEWER_STDOUT_WORLD)); 46 PetscCall(ISDestroy(&tmp)); 47 } 48 49 /* Set local matrix entries */ 50 PetscCall(ISGetLocalSize(isrows,&nrows)); 51 PetscCall(ISGetIndices(isrows,&rows)); 52 PetscCall(ISGetLocalSize(iscols,&ncols)); 53 PetscCall(ISGetIndices(iscols,&cols)); 54 PetscCall(PetscMalloc1(nrows*ncols,&v)); 55 for (i=0; i<nrows; i++) { 56 for (j=0; j<ncols; j++) { 57 /*v[i*ncols+j] = (PetscReal)(rank);*/ 58 v[i*ncols+j] = (PetscReal)(rank*10000+100*rows[i]+cols[j]); 59 } 60 } 61 PetscCall(MatSetValues(C,nrows,rows,ncols,cols,v,INSERT_VALUES)); 62 PetscCall(ISRestoreIndices(isrows,&rows)); 63 PetscCall(ISRestoreIndices(iscols,&cols)); 64 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 65 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 66 67 /* Test MatView() */ 68 if (mats_view) PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 69 70 /* Test MatMissingDiagonal() */ 71 PetscCall(MatMissingDiagonal(C,&flg,NULL)); 72 PetscCheck(!flg,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"MatMissingDiagonal() did not return false!"); 73 74 /* Set unowned matrix entries - add subdiagonals and diagonals from proc[0] */ 75 if (rank == 0) { 76 PetscInt M,N,cols[2]; 77 PetscCall(MatGetSize(C,&M,&N)); 78 for (i=0; i<M; i++) { 79 cols[0] = i; v[0] = i + 0.5; 80 cols[1] = i-1; v[1] = 0.5; 81 if (i) { 82 PetscCall(MatSetValues(C,1,&i,2,cols,v,ADD_VALUES)); 83 } else { 84 PetscCall(MatSetValues(C,1,&i,1,&i,v,ADD_VALUES)); 85 } 86 } 87 } 88 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 89 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 90 91 /* Test MatMult() */ 92 PetscCall(MatComputeOperator(C,MATAIJ,&Caij)); 93 PetscCall(MatMultEqual(C,Caij,5,&flg)); 94 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultEqual() fails"); 95 PetscCall(MatMultTransposeEqual(C,Caij,5,&flg)); 96 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultTransposeEqual() fails"); 97 98 /* Test MatMultAdd() and MatMultTransposeAddEqual() */ 99 PetscCall(MatMultAddEqual(C,Caij,5,&flg)); 100 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultAddEqual() fails"); 101 PetscCall(MatMultTransposeAddEqual(C,Caij,5,&flg)); 102 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultTransposeAddEqual() fails"); 103 104 /* Test MatMatMult() */ 105 PetscCall(PetscOptionsHasName(NULL,NULL,"-test_matmatmult",&Test_MatMatMult)); 106 if (Test_MatMatMult) { 107 Mat CCelem,CCaij; 108 PetscCall(MatMatMult(C,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CCelem)); 109 PetscCall(MatMatMult(Caij,Caij,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CCaij)); 110 PetscCall(MatMultEqual(CCelem,CCaij,5,&flg)); 111 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"CCelem != CCaij. MatMatMult() fails"); 112 PetscCall(MatDestroy(&CCaij)); 113 PetscCall(MatDestroy(&CCelem)); 114 } 115 116 PetscCall(PetscFree(v)); 117 PetscCall(ISDestroy(&isrows)); 118 PetscCall(ISDestroy(&iscols)); 119 PetscCall(MatDestroy(&C)); 120 PetscCall(MatDestroy(&Caij)); 121 PetscCall(PetscFinalize()); 122 return 0; 123 } 124 125 /*TEST 126 127 test: 128 nsize: 2 129 requires: elemental 130 args: -mat_type elemental -m 2 -n 3 131 132 test: 133 suffix: 2 134 nsize: 6 135 requires: elemental 136 args: -mat_type elemental -m 2 -n 2 137 138 test: 139 suffix: 3 140 nsize: 6 141 requires: elemental 142 args: -mat_type elemental -m 2 -n 2 -test_matmatmult 143 output_file: output/ex38_2.out 144 145 TEST*/ 146