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