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