1 static const char help[] = "Test parallel assembly of SBAIJ matrices\n\n"; 2 3 #include <petscmat.h> 4 5 PetscErrorCode Assemble(MPI_Comm comm,PetscInt bs,MatType mtype) 6 { 7 const PetscInt rc[] = {0,1,2,3}; 8 const PetscScalar vals[] = {100, 2, 3, 4, 5, 600, 7, 8, 9 9,100,11,1200,13,14,15,1600, 10 17,18,19,20,21,22,23,24, 11 25,26,27,2800,29,30,31,32, 12 33,34,35,36,37,38,39,40, 13 41,42,43,44,45,46,47,48, 14 49,50,51,52,53,54,55,56, 15 57,58,49,60,61,62,63,64}; 16 Mat A; 17 #if defined(PETSC_HAVE_MUMPS) || defined(PETSC_HAVE_MKL_CPARDISO) 18 Mat F; 19 MatSolverType stype; 20 PetscRandom rdm; 21 Vec b,x,y; 22 PetscInt i,j; 23 PetscReal norm2,tol=100*PETSC_SMALL; 24 PetscBool issbaij; 25 #endif 26 PetscViewer viewer; 27 PetscErrorCode ierr; 28 29 PetscFunctionBegin; 30 ierr = MatCreate(comm,&A);CHKERRQ(ierr); 31 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,4*bs,4*bs);CHKERRQ(ierr); 32 ierr = MatSetType(A,mtype);CHKERRQ(ierr); 33 ierr = MatMPIBAIJSetPreallocation(A,bs,2,NULL,2,NULL);CHKERRQ(ierr); 34 ierr = MatMPISBAIJSetPreallocation(A,bs,2,NULL,2,NULL);CHKERRQ(ierr); 35 ierr = MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 36 /* All processes contribute a global matrix */ 37 ierr = MatSetValuesBlocked(A,4,rc,4,rc,vals,ADD_VALUES);CHKERRQ(ierr); 38 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 39 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 40 ierr = PetscPrintf(comm,"Matrix %s(%D)\n",mtype,bs);CHKERRQ(ierr); 41 ierr = PetscViewerASCIIGetStdout(comm,&viewer);CHKERRQ(ierr); 42 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL);CHKERRQ(ierr); 43 ierr = MatView(A,viewer);CHKERRQ(ierr); 44 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 45 ierr = MatView(A,viewer);CHKERRQ(ierr); 46 #if defined(PETSC_HAVE_MUMPS) || defined(PETSC_HAVE_MKL_CPARDISO) 47 ierr = PetscStrcmp(mtype,MATMPISBAIJ,&issbaij);CHKERRQ(ierr); 48 if(!issbaij) { 49 ierr = MatShift(A,10);CHKERRQ(ierr); 50 } 51 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rdm);CHKERRQ(ierr); 52 ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr); 53 ierr = MatCreateVecs(A,&x,&y);CHKERRQ(ierr); 54 ierr = VecDuplicate(x,&b);CHKERRQ(ierr); 55 for (j=0; j<2; j++) { 56 #if defined(PETSC_HAVE_MUMPS) 57 if (j==0) stype = MATSOLVERMUMPS; 58 #else 59 if (j==0) continue; 60 #endif 61 #if defined(PETSC_HAVE_MKL_CPARDISO) 62 if (j==1) stype = MATSOLVERMKL_CPARDISO; 63 #else 64 if (j==1) continue; 65 #endif 66 if (issbaij) { 67 ierr = MatGetFactor(A,stype,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr); 68 ierr = MatCholeskyFactorSymbolic(F,A,NULL,NULL);CHKERRQ(ierr); 69 ierr = MatCholeskyFactorNumeric(F,A,NULL);CHKERRQ(ierr); 70 } else { 71 ierr = MatGetFactor(A,stype,MAT_FACTOR_LU,&F);CHKERRQ(ierr); 72 ierr = MatLUFactorSymbolic(F,A,NULL,NULL,NULL);CHKERRQ(ierr); 73 ierr = MatLUFactorNumeric(F,A,NULL);CHKERRQ(ierr); 74 } 75 for (i=0; i<10; i++) { 76 ierr = VecSetRandom(b,rdm);CHKERRQ(ierr); 77 ierr = MatSolve(F,b,y);CHKERRQ(ierr); 78 /* Check the error */ 79 ierr = MatMult(A,y,x);CHKERRQ(ierr); 80 ierr = VecAXPY(x,-1.0,b);CHKERRQ(ierr); 81 ierr = VecNorm(x,NORM_2,&norm2);CHKERRQ(ierr); 82 if (norm2>tol) { 83 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error:MatSolve(), norm2: %g\n",(double)norm2);CHKERRQ(ierr); 84 } 85 } 86 ierr = MatDestroy(&F);CHKERRQ(ierr); 87 } 88 ierr = VecDestroy(&x);CHKERRQ(ierr); 89 ierr = VecDestroy(&y);CHKERRQ(ierr); 90 ierr = VecDestroy(&b);CHKERRQ(ierr); 91 ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr); 92 #endif 93 ierr = MatDestroy(&A);CHKERRQ(ierr); 94 PetscFunctionReturn(0); 95 } 96 97 int main(int argc,char *argv[]) 98 { 99 PetscErrorCode ierr; 100 MPI_Comm comm; 101 PetscMPIInt size; 102 103 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 104 comm = PETSC_COMM_WORLD; 105 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 106 if (size != 2) SETERRQ(comm,PETSC_ERR_USER,"This example must be run with exactly two processes"); 107 ierr = Assemble(comm,2,MATMPIBAIJ);CHKERRQ(ierr); 108 ierr = Assemble(comm,2,MATMPISBAIJ);CHKERRQ(ierr); 109 ierr = Assemble(comm,1,MATMPIBAIJ);CHKERRQ(ierr); 110 ierr = Assemble(comm,1,MATMPISBAIJ);CHKERRQ(ierr); 111 ierr = PetscFinalize(); 112 return ierr; 113 } 114 115 116 /*TEST 117 118 test: 119 nsize: 2 120 args: -mat_ignore_lower_triangular -vecscatter_type sf 121 filter: sed -e "s~mem [0-9]*~mem~g" 122 123 TEST*/ 124