Actual source code: ex139.c


  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: {
  8:   IS             istmp;

 10:   PetscPrintf(PETSC_COMM_WORLD,"Extracting LocalRef with isrow:\n");
 11:   ISOnComm(isrow,PETSC_COMM_WORLD,PETSC_COPY_VALUES,&istmp);
 12:   ISView(istmp,PETSC_VIEWER_STDOUT_WORLD);
 13:   ISDestroy(&istmp);
 14:   PetscPrintf(PETSC_COMM_WORLD,"Extracting LocalRef with iscol (only rank=0 shown):\n");
 15:   ISOnComm(iscol,PETSC_COMM_WORLD,PETSC_COPY_VALUES,&istmp);
 16:   ISView(istmp,PETSC_VIEWER_STDOUT_WORLD);
 17:   ISDestroy(&istmp);

 19:   MatCreateLocalRef(A,isrow,iscol,B);
 20:   return 0;
 21: }

 23: int main(int argc,char *argv[])
 24: {
 25:   PetscErrorCode         ierr;
 26:   MPI_Comm               comm;
 27:   Mat                    J,B;
 28:   PetscInt               i,j,k,l,rstart,rend,m,n,top_bs,row_bs,col_bs,nlocblocks,*idx,nrowblocks,ncolblocks,*ridx,*cidx,*icol,*irow;
 29:   PetscScalar            *vals;
 30:   ISLocalToGlobalMapping brmap;
 31:   IS                     is0,is1;
 32:   PetscBool              diag,blocked;

 34:   PetscInitialize(&argc,&argv,0,help);
 35:   comm = PETSC_COMM_WORLD;

 37:   PetscOptionsBegin(comm,NULL,"LocalRef Test Options",NULL);
 38:   {
 39:     top_bs = 2; row_bs = 2; col_bs = 2; diag = PETSC_FALSE; blocked = PETSC_FALSE;
 40:     PetscOptionsInt("-top_bs","Block size of top-level matrix",0,top_bs,&top_bs,NULL);
 41:     PetscOptionsInt("-row_bs","Block size of row map",0,row_bs,&row_bs,NULL);
 42:     PetscOptionsInt("-col_bs","Block size of col map",0,col_bs,&col_bs,NULL);
 43:     PetscOptionsBool("-diag","Extract a diagonal black",0,diag,&diag,NULL);
 44:     PetscOptionsBool("-blocked","Use block insertion",0,blocked,&blocked,NULL);
 45:   }
 46:   PetscOptionsEnd();

 48:   MatCreate(comm,&J);
 49:   MatSetSizes(J,6,6,PETSC_DETERMINE,PETSC_DETERMINE);
 50:   MatSetBlockSize(J,top_bs);
 51:   MatSetFromOptions(J);
 52:   MatSeqBAIJSetPreallocation(J,top_bs,PETSC_DECIDE,0);
 53:   MatMPIBAIJSetPreallocation(J,top_bs,PETSC_DECIDE,0,PETSC_DECIDE,0);
 54:   MatSetUp(J);
 55:   MatGetSize(J,&m,&n);
 56:   MatGetOwnershipRange(J,&rstart,&rend);

 58:   nlocblocks = (rend-rstart)/top_bs + 2;
 59:   PetscMalloc1(nlocblocks,&idx);
 60:   for (i=0; i<nlocblocks; i++) {
 61:     idx[i] = (rstart/top_bs + i - 1 + m/top_bs) % (m/top_bs);
 62:   }
 63:   ISLocalToGlobalMappingCreate(comm,top_bs,nlocblocks,idx,PETSC_OWN_POINTER,&brmap);
 64:   PetscPrintf(comm,"Block ISLocalToGlobalMapping:\n");
 65:   ISLocalToGlobalMappingView(brmap,PETSC_VIEWER_STDOUT_WORLD);

 67:   MatSetLocalToGlobalMapping(J,brmap,brmap);
 68:   ISLocalToGlobalMappingDestroy(&brmap);

 70:   /* Create index sets for local submatrix */
 71:   nrowblocks = (rend-rstart)/row_bs;
 72:   ncolblocks = (rend-rstart)/col_bs;
 73:   PetscMalloc2(nrowblocks,&ridx,ncolblocks,&cidx);
 74:   for (i=0; i<nrowblocks; i++) ridx[i] = i + ((i > nrowblocks/2) ^ !rstart);
 75:   for (i=0; i<ncolblocks; i++) cidx[i] = i + ((i < ncolblocks/2) ^ !!rstart);
 76:   ISCreateBlock(PETSC_COMM_SELF,row_bs,nrowblocks,ridx,PETSC_COPY_VALUES,&is0);
 77:   ISCreateBlock(PETSC_COMM_SELF,col_bs,ncolblocks,cidx,PETSC_COPY_VALUES,&is1);
 78:   PetscFree2(ridx,cidx);

 80:   if (diag) {
 81:     ISDestroy(&is1);
 82:     PetscObjectReference((PetscObject)is0);
 83:     is1        = is0;
 84:     ncolblocks = nrowblocks;
 85:   }

 87:   GetLocalRef(J,is0,is1,&B);

 89:   MatZeroEntries(J);

 91:   PetscMalloc3(row_bs,&irow,col_bs,&icol,row_bs*col_bs,&vals);
 92:   for (i=0; i<nrowblocks; i++) {
 93:     for (j=0; j<ncolblocks; j++) {
 94:       for (k=0; k<row_bs; k++) {
 95:         for (l=0; l<col_bs; l++) {
 96:           vals[k*col_bs+l] = i * 100000 + j * 1000 + k * 10 + l;
 97:         }
 98:         irow[k] = i*row_bs+k;
 99:       }
100:       for (l=0; l<col_bs; l++) icol[l] = j*col_bs+l;
101:       if (blocked) {
102:         MatSetValuesBlockedLocal(B,1,&i,1,&j,vals,ADD_VALUES);
103:       } else {
104:         MatSetValuesLocal(B,row_bs,irow,col_bs,icol,vals,ADD_VALUES);
105:       }
106:     }
107:   }
108:   PetscFree3(irow,icol,vals);

110:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
111:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);

113:   MatView(J,PETSC_VIEWER_STDOUT_WORLD);

115:   ISDestroy(&is0);
116:   ISDestroy(&is1);
117:   MatDestroy(&B);
118:   MatDestroy(&J);
119:   PetscFinalize();
120:   return 0;
121: }

123: /*TEST

125:    test:
126:       nsize: 2
127:       filter: grep -v "type: mpi"
128:       args: -blocked {{0 1}} -mat_type {{aij baij}}

130: TEST*/