1 2 static char help[] = "Tests MPI parallel matrix creation. Test MatCreateRedundantMatrix() \n\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **args) 7 { 8 Mat C,Credundant; 9 MatInfo info; 10 PetscMPIInt rank,size,subsize; 11 PetscInt i,j,m = 3,n = 2,low,high,iglobal; 12 PetscInt Ii,J,ldim,nsubcomms; 13 PetscBool flg_info,flg_mat; 14 PetscScalar v,one = 1.0; 15 Vec x,y; 16 MPI_Comm subcomm; 17 18 PetscFunctionBeginUser; 19 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 20 PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 21 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 22 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 23 n = 2*size; 24 25 PetscCall(MatCreate(PETSC_COMM_WORLD,&C)); 26 PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n)); 27 PetscCall(MatSetFromOptions(C)); 28 PetscCall(MatSetUp(C)); 29 30 /* Create the matrix for the five point stencil, YET AGAIN */ 31 for (i=0; i<m; i++) { 32 for (j=2*rank; j<2*rank+2; j++) { 33 v = -1.0; Ii = j + n*i; 34 if (i>0) {J = Ii - n; PetscCall(MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES));} 35 if (i<m-1) {J = Ii + n; PetscCall(MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES));} 36 if (j>0) {J = Ii - 1; PetscCall(MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES));} 37 if (j<n-1) {J = Ii + 1; PetscCall(MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES));} 38 v = 4.0; PetscCall(MatSetValues(C,1,&Ii,1,&Ii,&v,INSERT_VALUES)); 39 } 40 } 41 42 /* Add extra elements (to illustrate variants of MatGetInfo) */ 43 Ii = n; J = n-2; v = 100.0; 44 PetscCall(MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES)); 45 Ii = n-2; J = n; v = 100.0; 46 PetscCall(MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES)); 47 48 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 49 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 50 51 /* Form vectors */ 52 PetscCall(MatCreateVecs(C,&x,&y)); 53 PetscCall(VecGetLocalSize(x,&ldim)); 54 PetscCall(VecGetOwnershipRange(x,&low,&high)); 55 for (i=0; i<ldim; i++) { 56 iglobal = i + low; 57 v = one*((PetscReal)i) + 100.0*rank; 58 PetscCall(VecSetValues(x,1,&iglobal,&v,INSERT_VALUES)); 59 } 60 PetscCall(VecAssemblyBegin(x)); 61 PetscCall(VecAssemblyEnd(x)); 62 63 PetscCall(MatMult(C,x,y)); 64 65 PetscCall(PetscOptionsHasName(NULL,NULL,"-view_info",&flg_info)); 66 if (flg_info) { 67 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO)); 68 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 69 70 PetscCall(MatGetInfo(C,MAT_GLOBAL_SUM,&info)); 71 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"matrix information (global sums):\nnonzeros = %" PetscInt_FMT ", allocated nonzeros = %" PetscInt_FMT "\n",(PetscInt)info.nz_used,(PetscInt)info.nz_allocated)); 72 PetscCall(MatGetInfo (C,MAT_GLOBAL_MAX,&info)); 73 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"matrix information (global max):\nnonzeros = %" PetscInt_FMT ", allocated nonzeros = %" PetscInt_FMT "\n",(PetscInt)info.nz_used,(PetscInt)info.nz_allocated)); 74 } 75 76 PetscCall(PetscOptionsHasName(NULL,NULL,"-view_mat",&flg_mat)); 77 if (flg_mat) PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 78 79 /* Test MatCreateRedundantMatrix() */ 80 nsubcomms = size; 81 PetscCall(PetscOptionsGetInt(NULL,NULL,"-nsubcomms",&nsubcomms,NULL)); 82 PetscCall(MatCreateRedundantMatrix(C,nsubcomms,MPI_COMM_NULL,MAT_INITIAL_MATRIX,&Credundant)); 83 PetscCall(MatCreateRedundantMatrix(C,nsubcomms,MPI_COMM_NULL,MAT_REUSE_MATRIX,&Credundant)); 84 85 PetscCall(PetscObjectGetComm((PetscObject)Credundant,&subcomm)); 86 PetscCallMPI(MPI_Comm_size(subcomm,&subsize)); 87 88 if (subsize==2 && flg_mat) { 89 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(subcomm),"\n[%d] Credundant:\n",rank)); 90 PetscCall(MatView(Credundant,PETSC_VIEWER_STDOUT_(subcomm))); 91 } 92 PetscCall(MatDestroy(&Credundant)); 93 94 /* Test MatCreateRedundantMatrix() with user-provided subcomm */ 95 { 96 PetscSubcomm psubcomm; 97 98 PetscCall(PetscSubcommCreate(PETSC_COMM_WORLD,&psubcomm)); 99 PetscCall(PetscSubcommSetNumber(psubcomm,nsubcomms)); 100 PetscCall(PetscSubcommSetType(psubcomm,PETSC_SUBCOMM_CONTIGUOUS)); 101 /* enable runtime switch of psubcomm type, e.g., '-psubcomm_type interlaced */ 102 PetscCall(PetscSubcommSetFromOptions(psubcomm)); 103 104 PetscCall(MatCreateRedundantMatrix(C,nsubcomms,PetscSubcommChild(psubcomm),MAT_INITIAL_MATRIX,&Credundant)); 105 PetscCall(MatCreateRedundantMatrix(C,nsubcomms,PetscSubcommChild(psubcomm),MAT_REUSE_MATRIX,&Credundant)); 106 107 PetscCall(PetscSubcommDestroy(&psubcomm)); 108 PetscCall(MatDestroy(&Credundant)); 109 } 110 111 PetscCall(VecDestroy(&x)); 112 PetscCall(VecDestroy(&y)); 113 PetscCall(MatDestroy(&C)); 114 PetscCall(PetscFinalize()); 115 return 0; 116 } 117 118 /*TEST 119 120 test: 121 nsize: 3 122 args: -view_info 123 124 test: 125 suffix: 2 126 nsize: 3 127 args: -nsubcomms 2 -view_mat -psubcomm_type interlaced 128 129 test: 130 suffix: 3 131 nsize: 3 132 args: -nsubcomms 2 -view_mat -psubcomm_type contiguous 133 134 test: 135 suffix: 3_baij 136 nsize: 3 137 args: -mat_type baij -nsubcomms 2 -view_mat 138 139 test: 140 suffix: 3_sbaij 141 nsize: 3 142 args: -mat_type sbaij -nsubcomms 2 -view_mat 143 144 test: 145 suffix: 3_dense 146 nsize: 3 147 args: -mat_type dense -nsubcomms 2 -view_mat 148 149 test: 150 suffix: 4_baij 151 nsize: 3 152 args: -mat_type baij -nsubcomms 2 -view_mat -psubcomm_type interlaced 153 154 test: 155 suffix: 4_sbaij 156 nsize: 3 157 args: -mat_type sbaij -nsubcomms 2 -view_mat -psubcomm_type interlaced 158 159 test: 160 suffix: 4_dense 161 nsize: 3 162 args: -mat_type dense -nsubcomms 2 -view_mat -psubcomm_type interlaced 163 164 TEST*/ 165