1 2 const char help[] = "Test MatCreateLocalRef()\n\n"; 3 4 #include <petscmat.h> 5 6 static PetscErrorCode GetLocalRef(Mat A,IS isrow,IS iscol,Mat *B) 7 { 8 PetscErrorCode ierr; 9 IS istmp; 10 11 PetscFunctionBegin; 12 ierr = PetscPrintf(PETSC_COMM_WORLD,"Extracting LocalRef with isrow:\n");CHKERRQ(ierr); 13 ierr = ISOnComm(isrow,PETSC_COMM_WORLD,PETSC_COPY_VALUES,&istmp);CHKERRQ(ierr); 14 ierr = ISView(istmp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 15 ierr = ISDestroy(&istmp);CHKERRQ(ierr); 16 ierr = PetscPrintf(PETSC_COMM_WORLD,"Extracting LocalRef with iscol (only rank=0 shown):\n");CHKERRQ(ierr); 17 ierr = ISOnComm(iscol,PETSC_COMM_WORLD,PETSC_COPY_VALUES,&istmp);CHKERRQ(ierr); 18 ierr = ISView(istmp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 19 ierr = ISDestroy(&istmp);CHKERRQ(ierr); 20 21 ierr = MatCreateLocalRef(A,isrow,iscol,B);CHKERRQ(ierr); 22 PetscFunctionReturn(0); 23 } 24 25 int main(int argc,char *argv[]) 26 { 27 PetscErrorCode ierr; 28 MPI_Comm comm; 29 Mat J,B; 30 PetscInt i,j,k,l,rstart,rend,m,n,top_bs,row_bs,col_bs,nlocblocks,*idx,nrowblocks,ncolblocks,*ridx,*cidx,*icol,*irow; 31 PetscScalar *vals; 32 ISLocalToGlobalMapping brmap; 33 IS is0,is1; 34 PetscBool diag,blocked; 35 36 ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr; 37 comm = PETSC_COMM_WORLD; 38 39 ierr = PetscOptionsBegin(comm,NULL,"LocalRef Test Options",NULL);CHKERRQ(ierr); 40 { 41 top_bs = 2; row_bs = 2; col_bs = 2; diag = PETSC_FALSE; blocked = PETSC_FALSE; 42 ierr = PetscOptionsInt("-top_bs","Block size of top-level matrix",0,top_bs,&top_bs,NULL);CHKERRQ(ierr); 43 ierr = PetscOptionsInt("-row_bs","Block size of row map",0,row_bs,&row_bs,NULL);CHKERRQ(ierr); 44 ierr = PetscOptionsInt("-col_bs","Block size of col map",0,col_bs,&col_bs,NULL);CHKERRQ(ierr); 45 ierr = PetscOptionsBool("-diag","Extract a diagonal black",0,diag,&diag,NULL);CHKERRQ(ierr); 46 ierr = PetscOptionsBool("-blocked","Use block insertion",0,blocked,&blocked,NULL);CHKERRQ(ierr); 47 } 48 ierr = PetscOptionsEnd();CHKERRQ(ierr); 49 50 ierr = MatCreate(comm,&J);CHKERRQ(ierr); 51 ierr = MatSetSizes(J,6,6,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 52 ierr = MatSetBlockSize(J,top_bs);CHKERRQ(ierr); 53 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 54 ierr = MatSeqBAIJSetPreallocation(J,top_bs,PETSC_DECIDE,0);CHKERRQ(ierr); 55 ierr = MatMPIBAIJSetPreallocation(J,top_bs,PETSC_DECIDE,0,PETSC_DECIDE,0);CHKERRQ(ierr); 56 ierr = MatSetUp(J);CHKERRQ(ierr); 57 ierr = MatGetSize(J,&m,&n);CHKERRQ(ierr); 58 ierr = MatGetOwnershipRange(J,&rstart,&rend);CHKERRQ(ierr); 59 60 nlocblocks = (rend-rstart)/top_bs + 2; 61 ierr = PetscMalloc1(nlocblocks,&idx);CHKERRQ(ierr); 62 for (i=0; i<nlocblocks; i++) { 63 idx[i] = (rstart/top_bs + i - 1 + m/top_bs) % (m/top_bs); 64 } 65 ierr = ISLocalToGlobalMappingCreate(comm,top_bs,nlocblocks,idx,PETSC_OWN_POINTER,&brmap);CHKERRQ(ierr); 66 ierr = PetscPrintf(comm,"Block ISLocalToGlobalMapping:\n");CHKERRQ(ierr); 67 ierr = ISLocalToGlobalMappingView(brmap,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 68 69 ierr = MatSetLocalToGlobalMapping(J,brmap,brmap);CHKERRQ(ierr); 70 ierr = ISLocalToGlobalMappingDestroy(&brmap);CHKERRQ(ierr); 71 72 /* Create index sets for local submatrix */ 73 nrowblocks = (rend-rstart)/row_bs; 74 ncolblocks = (rend-rstart)/col_bs; 75 ierr = PetscMalloc2(nrowblocks,&ridx,ncolblocks,&cidx);CHKERRQ(ierr); 76 for (i=0; i<nrowblocks; i++) ridx[i] = i + ((i > nrowblocks/2) ^ !rstart); 77 for (i=0; i<ncolblocks; i++) cidx[i] = i + ((i < ncolblocks/2) ^ !!rstart); 78 ierr = ISCreateBlock(PETSC_COMM_SELF,row_bs,nrowblocks,ridx,PETSC_COPY_VALUES,&is0);CHKERRQ(ierr); 79 ierr = ISCreateBlock(PETSC_COMM_SELF,col_bs,ncolblocks,cidx,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr); 80 ierr = PetscFree2(ridx,cidx);CHKERRQ(ierr); 81 82 if (diag) { 83 ierr = ISDestroy(&is1);CHKERRQ(ierr); 84 ierr = PetscObjectReference((PetscObject)is0);CHKERRQ(ierr); 85 is1 = is0; 86 ncolblocks = nrowblocks; 87 } 88 89 ierr = GetLocalRef(J,is0,is1,&B);CHKERRQ(ierr); 90 91 ierr = MatZeroEntries(J);CHKERRQ(ierr); 92 93 ierr = PetscMalloc3(row_bs,&irow,col_bs,&icol,row_bs*col_bs,&vals);CHKERRQ(ierr); 94 for (i=0; i<nrowblocks; i++) { 95 for (j=0; j<ncolblocks; j++) { 96 for (k=0; k<row_bs; k++) { 97 for (l=0; l<col_bs; l++) { 98 vals[k*col_bs+l] = i * 100000 + j * 1000 + k * 10 + l; 99 } 100 irow[k] = i*row_bs+k; 101 } 102 for (l=0; l<col_bs; l++) icol[l] = j*col_bs+l; 103 if (blocked) { 104 ierr = MatSetValuesBlockedLocal(B,1,&i,1,&j,vals,ADD_VALUES);CHKERRQ(ierr); 105 } else { 106 ierr = MatSetValuesLocal(B,row_bs,irow,col_bs,icol,vals,ADD_VALUES);CHKERRQ(ierr); 107 } 108 } 109 } 110 ierr = PetscFree3(irow,icol,vals);CHKERRQ(ierr); 111 112 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 113 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 114 115 ierr = MatView(J,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 116 117 ierr = ISDestroy(&is0);CHKERRQ(ierr); 118 ierr = ISDestroy(&is1);CHKERRQ(ierr); 119 ierr = MatDestroy(&B);CHKERRQ(ierr); 120 ierr = MatDestroy(&J);CHKERRQ(ierr); 121 ierr = PetscFinalize(); 122 return ierr; 123 } 124 125 /*TEST 126 127 test: 128 nsize: 2 129 filter: grep -v "type: mpi" 130 args: -blocked {{0 1}} -mat_type {{aij baij}} 131 132 TEST*/ 133