1 2 static char help[] = "Test repeated LU factorizations. Used for checking memory leak\n\ 3 -m <size> : problem size\n\ 4 -mat_nonsym : use nonsymmetric matrix (default is symmetric)\n\n"; 5 6 #include <petscmat.h> 7 int main(int argc,char **args) 8 { 9 Mat C,F; /* matrix */ 10 Vec x,u,b; /* approx solution, RHS, exact solution */ 11 PetscReal norm; /* norm of solution error */ 12 PetscScalar v,none = -1.0; 13 PetscInt I,J,ldim,low,high,iglobal,Istart,Iend; 14 PetscErrorCode ierr; 15 PetscInt i,j,m = 3,n = 2,its; 16 PetscMPIInt size,rank; 17 PetscBool mat_nonsymmetric; 18 PetscInt its_max; 19 MatFactorInfo factinfo; 20 IS perm,iperm; 21 22 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 23 ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr); 24 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 25 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 26 n = 2*size; 27 28 /* 29 Set flag if we are doing a nonsymmetric problem; the default is symmetric. 30 */ 31 ierr = PetscOptionsHasName(NULL,NULL,"-mat_nonsym",&mat_nonsymmetric);CHKERRQ(ierr); 32 33 /* 34 Create parallel matrix, specifying only its global dimensions. 35 When using MatCreate(), the matrix format can be specified at 36 runtime. Also, the parallel partitioning of the matrix is 37 determined by PETSc at runtime. 38 */ 39 ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr); 40 ierr = MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr); 41 ierr = MatSetFromOptions(C);CHKERRQ(ierr); 42 ierr = MatGetOwnershipRange(C,&Istart,&Iend);CHKERRQ(ierr); 43 44 /* 45 Set matrix entries matrix in parallel. 46 - Each processor needs to insert only elements that it owns 47 locally (but any non-local elements will be sent to the 48 appropriate processor during matrix assembly). 49 - Always specify global row and columns of matrix entries. 50 */ 51 for (I=Istart; I<Iend; I++) { 52 v = -1.0; i = I/n; j = I - i*n; 53 if (i>0) {J = I - n; ierr = MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} 54 if (i<m-1) {J = I + n; ierr = MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} 55 if (j>0) {J = I - 1; ierr = MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} 56 if (j<n-1) {J = I + 1; ierr = MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} 57 v = 4.0; ierr = MatSetValues(C,1,&I,1,&I,&v,ADD_VALUES); 58 } 59 60 /* 61 Make the matrix nonsymmetric if desired 62 */ 63 if (mat_nonsymmetric) { 64 for (I=Istart; I<Iend; I++) { 65 v = -1.5; i = I/n; 66 if (i>1) {J = I-n-1; ierr = MatSetValues(C,1,&I,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} 67 } 68 } else { 69 ierr = MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 70 ierr = MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);CHKERRQ(ierr); 71 } 72 73 /* 74 Assemble matrix, using the 2-step process: 75 MatAssemblyBegin(), MatAssemblyEnd() 76 Computations can be done while messages are in transition 77 by placing code between these two statements. 78 */ 79 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 80 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 81 82 its_max=1000; 83 /* 84 Create parallel vectors. 85 - When using VecSetSizes(), we specify only the vector's global 86 dimension; the parallel partitioning is determined at runtime. 87 - Note: We form 1 vector from scratch and then duplicate as needed. 88 */ 89 ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr); 90 ierr = VecSetSizes(u,PETSC_DECIDE,m*n);CHKERRQ(ierr); 91 ierr = VecSetFromOptions(u);CHKERRQ(ierr); 92 ierr = VecDuplicate(u,&b);CHKERRQ(ierr); 93 ierr = VecDuplicate(b,&x);CHKERRQ(ierr); 94 95 /* 96 Currently, all parallel PETSc vectors are partitioned by 97 contiguous chunks across the processors. Determine which 98 range of entries are locally owned. 99 */ 100 ierr = VecGetOwnershipRange(x,&low,&high);CHKERRQ(ierr); 101 102 /* 103 Set elements within the exact solution vector in parallel. 104 - Each processor needs to insert only elements that it owns 105 locally (but any non-local entries will be sent to the 106 appropriate processor during vector assembly). 107 - Always specify global locations of vector entries. 108 */ 109 ierr = VecGetLocalSize(x,&ldim);CHKERRQ(ierr); 110 for (i=0; i<ldim; i++) { 111 iglobal = i + low; 112 v = (PetscScalar)(i + 100*rank); 113 ierr = VecSetValues(u,1,&iglobal,&v,INSERT_VALUES);CHKERRQ(ierr); 114 } 115 116 /* 117 Assemble vector, using the 2-step process: 118 VecAssemblyBegin(), VecAssemblyEnd() 119 Computations can be done while messages are in transition, 120 by placing code between these two statements. 121 */ 122 ierr = VecAssemblyBegin(u);CHKERRQ(ierr); 123 ierr = VecAssemblyEnd(u);CHKERRQ(ierr); 124 125 /* Compute right-hand-side vector */ 126 ierr = MatMult(C,u,b);CHKERRQ(ierr); 127 128 ierr = MatGetOrdering(C,MATORDERINGNATURAL,&perm,&iperm);CHKERRQ(ierr); 129 its_max = 2000; 130 for (i=0; i<its_max; i++) { 131 ierr = MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_LU,&F);CHKERRQ(ierr); 132 ierr = MatLUFactorSymbolic(F,C,perm,iperm,&factinfo);CHKERRQ(ierr); 133 for (j=0; j<1; j++) { 134 ierr = MatLUFactorNumeric(F,C,&factinfo);CHKERRQ(ierr); 135 } 136 ierr = MatSolve(F,b,x);CHKERRQ(ierr); 137 ierr = MatDestroy(&F);CHKERRQ(ierr); 138 } 139 ierr = ISDestroy(&perm);CHKERRQ(ierr); 140 ierr = ISDestroy(&iperm);CHKERRQ(ierr); 141 142 /* Check the error */ 143 ierr = VecAXPY(x,none,u);CHKERRQ(ierr); 144 ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr); 145 ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %t\n",(double)norm);CHKERRQ(ierr); 146 147 /* Free work space. */ 148 ierr = VecDestroy(&u);CHKERRQ(ierr); 149 ierr = VecDestroy(&x);CHKERRQ(ierr); 150 ierr = VecDestroy(&b);CHKERRQ(ierr); 151 ierr = MatDestroy(&C);CHKERRQ(ierr); 152 ierr = PetscFinalize(); 153 return ierr; 154 } 155 156 157