1 2 static char help[] = "Tests MatTranspose() and MatEqual() for MPIAIJ matrices.\n\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **argv) 7 { 8 Mat A,B; 9 PetscInt m = 7,n,i,rstart,rend,cols[3]; 10 PetscScalar v[3]; 11 PetscBool equal; 12 const char *eq[2]; 13 14 PetscFunctionBeginUser; 15 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 16 PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 17 n = m; 18 19 /* ------- Assemble matrix, --------- */ 20 21 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,m,n,0,0,0,0,&A)); 22 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE)); 23 PetscCall(MatGetOwnershipRange(A,&rstart,&rend)); 24 if (!rstart) { 25 cols[0] = 0; 26 cols[1] = 1; 27 v[0] = 2.0; v[1] = -1.0; 28 PetscCall(MatSetValues(A,1,&rstart,2,cols,v,INSERT_VALUES)); 29 rstart++; 30 } 31 if (rend == m) { 32 rend--; 33 cols[0] = rend-1; 34 cols[1] = rend; 35 v[0] = -1.0; v[1] = 2.0; 36 PetscCall(MatSetValues(A,1,&rend,2,cols,v,INSERT_VALUES)); 37 } 38 v[0] = -1.0; v[1] = 2.0; v[2] = -1.0; 39 for (i=rstart; i<rend; i++) { 40 cols[0] = i-1; 41 cols[1] = i; 42 cols[2] = i+1; 43 PetscCall(MatSetValues(A,1,&i,3,cols,v,INSERT_VALUES)); 44 } 45 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 46 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 47 48 PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX,&B)); 49 50 PetscCall(MatEqual(A,B,&equal)); 51 52 eq[0] = "not equal"; 53 eq[1] = "equal"; 54 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Matrices are %s\n",eq[equal])); 55 56 PetscCall(MatTranspose(A,MAT_REUSE_MATRIX,&B)); 57 PetscCall(MatEqual(A,B,&equal)); 58 if (!equal) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatTranspose with MAT_REUSE_MATRIX failed")); 59 60 /* Free data structures */ 61 PetscCall(MatDestroy(&A)); 62 PetscCall(MatDestroy(&B)); 63 64 PetscCall(PetscFinalize()); 65 return 0; 66 } 67 68 /*TEST 69 70 test: 71 72 test: 73 suffix: 2 74 nsize: 2 75 output_file: output/ex58_1.out 76 77 TEST*/ 78