Actual source code: ex139.c
petsc-3.14.6 2021-03-30
2: const char help[] = "Test MatCreateLocalRef()\n\n";
4: #include <petscmat.h>
6: static PetscErrorCode GetLocalRef(Mat A,IS isrow,IS iscol,Mat *B)
7: {
9: IS istmp;
12: PetscPrintf(PETSC_COMM_WORLD,"Extracting LocalRef with isrow:\n");
13: ISOnComm(isrow,PETSC_COMM_WORLD,PETSC_COPY_VALUES,&istmp);
14: ISView(istmp,PETSC_VIEWER_STDOUT_WORLD);
15: ISDestroy(&istmp);
16: PetscPrintf(PETSC_COMM_WORLD,"Extracting LocalRef with iscol (only rank=0 shown):\n");
17: ISOnComm(iscol,PETSC_COMM_WORLD,PETSC_COPY_VALUES,&istmp);
18: ISView(istmp,PETSC_VIEWER_STDOUT_WORLD);
19: ISDestroy(&istmp);
21: MatCreateLocalRef(A,isrow,iscol,B);
22: return(0);
23: }
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;
36: PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
37: comm = PETSC_COMM_WORLD;
39: PetscOptionsBegin(comm,NULL,"LocalRef Test Options",NULL);
40: {
41: top_bs = 2; row_bs = 2; col_bs = 2; diag = PETSC_FALSE; blocked = PETSC_FALSE;
42: PetscOptionsInt("-top_bs","Block size of top-level matrix",0,top_bs,&top_bs,NULL);
43: PetscOptionsInt("-row_bs","Block size of row map",0,row_bs,&row_bs,NULL);
44: PetscOptionsInt("-col_bs","Block size of col map",0,col_bs,&col_bs,NULL);
45: PetscOptionsBool("-diag","Extract a diagonal black",0,diag,&diag,NULL);
46: PetscOptionsBool("-blocked","Use block insertion",0,blocked,&blocked,NULL);
47: }
48: PetscOptionsEnd();
50: MatCreate(comm,&J);
51: MatSetSizes(J,6,6,PETSC_DETERMINE,PETSC_DETERMINE);
52: MatSetBlockSize(J,top_bs);
53: MatSetFromOptions(J);
54: MatSeqBAIJSetPreallocation(J,top_bs,PETSC_DECIDE,0);
55: MatMPIBAIJSetPreallocation(J,top_bs,PETSC_DECIDE,0,PETSC_DECIDE,0);
56: MatSetUp(J);
57: MatGetSize(J,&m,&n);
58: MatGetOwnershipRange(J,&rstart,&rend);
60: nlocblocks = (rend-rstart)/top_bs + 2;
61: PetscMalloc1(nlocblocks,&idx);
62: for (i=0; i<nlocblocks; i++) {
63: idx[i] = (rstart/top_bs + i - 1 + m/top_bs) % (m/top_bs);
64: }
65: ISLocalToGlobalMappingCreate(comm,top_bs,nlocblocks,idx,PETSC_OWN_POINTER,&brmap);
66: PetscPrintf(comm,"Block ISLocalToGlobalMapping:\n");
67: ISLocalToGlobalMappingView(brmap,PETSC_VIEWER_STDOUT_WORLD);
69: MatSetLocalToGlobalMapping(J,brmap,brmap);
70: ISLocalToGlobalMappingDestroy(&brmap);
72: /* Create index sets for local submatrix */
73: nrowblocks = (rend-rstart)/row_bs;
74: ncolblocks = (rend-rstart)/col_bs;
75: PetscMalloc2(nrowblocks,&ridx,ncolblocks,&cidx);
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: ISCreateBlock(PETSC_COMM_SELF,row_bs,nrowblocks,ridx,PETSC_COPY_VALUES,&is0);
79: ISCreateBlock(PETSC_COMM_SELF,col_bs,ncolblocks,cidx,PETSC_COPY_VALUES,&is1);
80: PetscFree2(ridx,cidx);
82: if (diag) {
83: ISDestroy(&is1);
84: PetscObjectReference((PetscObject)is0);
85: is1 = is0;
86: ncolblocks = nrowblocks;
87: }
89: GetLocalRef(J,is0,is1,&B);
91: MatZeroEntries(J);
93: PetscMalloc3(row_bs,&irow,col_bs,&icol,row_bs*col_bs,&vals);
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: MatSetValuesBlockedLocal(B,1,&i,1,&j,vals,ADD_VALUES);
105: } else {
106: MatSetValuesLocal(B,row_bs,irow,col_bs,icol,vals,ADD_VALUES);
107: }
108: }
109: }
110: PetscFree3(irow,icol,vals);
112: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
113: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
115: MatView(J,PETSC_VIEWER_STDOUT_WORLD);
117: ISDestroy(&is0);
118: ISDestroy(&is1);
119: MatDestroy(&B);
120: MatDestroy(&J);
121: PetscFinalize();
122: return ierr;
123: }
125: /*TEST
127: test:
128: nsize: 2
129: filter: grep -v "type: mpi"
130: args: -blocked {{0 1}} -mat_type {{aij baij}}
132: TEST*/