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