Actual source code: matis.c

      1: #define PETSCMAT_DLL

      3: /*
      4:     Creates a matrix class for using the Neumann-Neumann type preconditioners.
      5:    This stores the matrices in globally unassembled form. Each processor
      6:    assembles only its local Neumann problem and the parallel matrix vector
      7:    product is handled "implicitly".

      9:      We provide:
     10:          MatMult()

     12:     Currently this allows for only one subdomain per processor.

     14: */

     16:  #include src/mat/impls/is/matis.h

     20: PetscErrorCode MatDestroy_IS(Mat A)
     21: {
     23:   Mat_IS         *b = (Mat_IS*)A->data;

     26:   if (b->A) {
     27:     MatDestroy(b->A);
     28:   }
     29:   if (b->ctx) {
     30:     VecScatterDestroy(b->ctx);
     31:   }
     32:   if (b->x) {
     33:     VecDestroy(b->x);
     34:   }
     35:   if (b->y) {
     36:     VecDestroy(b->y);
     37:   }
     38:   if (b->mapping) {
     39:     ISLocalToGlobalMappingDestroy(b->mapping);
     40:   }
     41:   PetscFree(b);
     42:   PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C","",PETSC_NULL);
     43:   return(0);
     44: }

     48: PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y)
     49: {
     51:   Mat_IS         *is = (Mat_IS*)A->data;
     52:   PetscScalar    zero = 0.0;

     55:   /*  scatter the global vector x into the local work vector */
     56:   VecScatterBegin(x,is->x,INSERT_VALUES,SCATTER_FORWARD,is->ctx);
     57:   VecScatterEnd(x,is->x,INSERT_VALUES,SCATTER_FORWARD,is->ctx);

     59:   /* multiply the local matrix */
     60:   MatMult(is->A,is->x,is->y);

     62:   /* scatter product back into global memory */
     63:   VecSet(y,zero);
     64:   VecScatterBegin(is->y,y,ADD_VALUES,SCATTER_REVERSE,is->ctx);
     65:   VecScatterEnd(is->y,y,ADD_VALUES,SCATTER_REVERSE,is->ctx);

     67:   return(0);
     68: }

     72: PetscErrorCode MatView_IS(Mat A,PetscViewer viewer)
     73: {
     74:   Mat_IS         *a = (Mat_IS*)A->data;
     76:   PetscViewer    sviewer;

     79:   PetscViewerGetSingleton(viewer,&sviewer);
     80:   MatView(a->A,sviewer);
     81:   PetscViewerRestoreSingleton(viewer,&sviewer);
     82:   return(0);
     83: }

     87: PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping mapping)
     88: {
     90:   PetscInt       n;
     91:   Mat_IS         *is = (Mat_IS*)A->data;
     92:   IS             from,to;
     93:   Vec            global;

     96:   is->mapping = mapping;
     97:   PetscObjectReference((PetscObject)mapping);

     99:   /* Create the local matrix A */
    100:   ISLocalToGlobalMappingGetSize(mapping,&n);
    101:   MatCreate(PETSC_COMM_SELF,&is->A);
    102:   MatSetSizes(is->A,n,n,n,n);
    103:   PetscObjectSetOptionsPrefix((PetscObject)is->A,"is");
    104:   MatSetFromOptions(is->A);

    106:   /* Create the local work vectors */
    107:   VecCreateSeq(PETSC_COMM_SELF,n,&is->x);
    108:   VecDuplicate(is->x,&is->y);

    110:   /* setup the global to local scatter */
    111:   ISCreateStride(PETSC_COMM_SELF,n,0,1,&to);
    112:   ISLocalToGlobalMappingApplyIS(mapping,to,&from);
    113:   VecCreateMPI(A->comm,A->n,A->N,&global);
    114:   VecScatterCreate(global,from,is->x,to,&is->ctx);
    115:   VecDestroy(global);
    116:   ISDestroy(to);
    117:   ISDestroy(from);
    118:   return(0);
    119: }

    124: PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv)
    125: {
    127:   Mat_IS         *is = (Mat_IS*)A->data;

    130:   MatSetValues(is->A,m,rows,n,cols,values,addv);
    131:   return(0);
    132: }

    136: PetscErrorCode MatZeroRowsLocal_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag)
    137: {
    138:   Mat_IS         *is = (Mat_IS*)A->data;
    140:   PetscInt       i;
    141:   PetscScalar    *array;

    144:   {
    145:     /*
    146:        Set up is->x as a "counting vector". This is in order to MatMult_IS
    147:        work properly in the interface nodes.
    148:     */
    149:     Vec         counter;
    150:     PetscScalar one=1.0, zero=0.0;
    151:     VecCreateMPI(A->comm,A->n,A->N,&counter);
    152:     VecSet(counter,zero);
    153:     VecSet(is->x,one);
    154:     VecScatterBegin(is->x,counter,ADD_VALUES,SCATTER_REVERSE,is->ctx);
    155:     VecScatterEnd  (is->x,counter,ADD_VALUES,SCATTER_REVERSE,is->ctx);
    156:     VecScatterBegin(counter,is->x,INSERT_VALUES,SCATTER_FORWARD,is->ctx);
    157:     VecScatterEnd  (counter,is->x,INSERT_VALUES,SCATTER_FORWARD,is->ctx);
    158:     VecDestroy(counter);
    159:   }
    160:   if (!n) {
    161:     is->pure_neumann = PETSC_TRUE;
    162:   } else {
    163:     is->pure_neumann = PETSC_FALSE;
    164:     VecGetArray(is->x,&array);
    165:     MatZeroRows(is->A,n,rows,diag);
    166:     for (i=0; i<n; i++) {
    167:       MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);
    168:     }
    169:     MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);
    170:     MatAssemblyEnd  (is->A,MAT_FINAL_ASSEMBLY);
    171:     VecRestoreArray(is->x,&array);
    172:   }

    174:   return(0);
    175: }

    179: PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type)
    180: {
    181:   Mat_IS         *is = (Mat_IS*)A->data;

    185:   MatAssemblyBegin(is->A,type);
    186:   return(0);
    187: }

    191: PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type)
    192: {
    193:   Mat_IS         *is = (Mat_IS*)A->data;

    197:   MatAssemblyEnd(is->A,type);
    198:   return(0);
    199: }

    204: PetscErrorCode PETSCMAT_DLLEXPORT MatISGetLocalMat_IS(Mat mat,Mat *local)
    205: {
    206:   Mat_IS *is = (Mat_IS *)mat->data;
    207:
    209:   *local = is->A;
    210:   return(0);
    211: }

    216: /*@
    217:     MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix.

    219:   Input Parameter:
    220: .  mat - the matrix

    222:   Output Parameter:
    223: .  local - the local matrix usually MATSEQAIJ

    225:   Level: advanced

    227:   Notes:
    228:     This can be called if you have precomputed the nonzero structure of the
    229:   matrix and want to provide it to the inner matrix object to improve the performance
    230:   of the MatSetValues() operation.

    232: .seealso: MATIS
    233: @*/
    234: PetscErrorCode PETSCMAT_DLLEXPORT MatISGetLocalMat(Mat mat,Mat *local)
    235: {
    236:   PetscErrorCode ierr,(*f)(Mat,Mat *);

    241:   PetscObjectQueryFunction((PetscObject)mat,"MatISGetLocalMat_C",(void (**)(void))&f);
    242:   if (f) {
    243:     (*f)(mat,local);
    244:   } else {
    245:     local = 0;
    246:   }
    247:   return(0);
    248: }

    252: PetscErrorCode MatZeroEntries_IS(Mat A)
    253: {
    254:   Mat_IS         *a = (Mat_IS*)A->data;

    258:   MatZeroEntries(a->A);
    259:   return(0);
    260: }

    264: PetscErrorCode MatSetOption_IS(Mat A,MatOption op)
    265: {
    266:   Mat_IS         *a = (Mat_IS*)A->data;

    270:   MatSetOption(a->A,op);
    271:   return(0);
    272: }

    274: /*MC
    275:    MATIS - MATIS = "is" - A matrix type to be used for using the Neumann-Neumann type preconditioners.
    276:    This stores the matrices in globally unassembled form. Each processor
    277:    assembles only its local Neumann problem and the parallel matrix vector
    278:    product is handled "implicitly".

    280:    Operations Provided:
    281: +  MatMult()
    282: .  MatZeroEntries()
    283: .  MatSetOption()
    284: .  MatZeroRowsLocal()
    285: .  MatSetValuesLocal()
    286: -  MatSetLocalToGlobalMapping()

    288:    Options Database Keys:
    289: . -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions()

    291:    Notes: Options prefix for the inner matrix are given by -is_mat_xxx
    292:
    293:           You must call MatSetLocalToGlobalMapping() before using this matrix type.

    295:           You can do matrix preallocation on the local matrix after you obtain it with
    296:           MatISGetLocalMat()

    298:   Level: advanced

    300: .seealso: PC, MatISGetLocalMat(), MatSetLocalToGlobalMapping()

    302: M*/

    307: PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_IS(Mat A)
    308: {
    310:   Mat_IS         *b;

    313:   PetscNew(Mat_IS,&b);
    314:   A->data             = (void*)b;
    315:   PetscMemzero(A->ops,sizeof(struct _MatOps));
    316:   A->factor           = 0;
    317:   A->mapping          = 0;

    319:   A->ops->mult                    = MatMult_IS;
    320:   A->ops->destroy                 = MatDestroy_IS;
    321:   A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS;
    322:   A->ops->setvalueslocal          = MatSetValuesLocal_IS;
    323:   A->ops->zerorowslocal           = MatZeroRowsLocal_IS;
    324:   A->ops->assemblybegin           = MatAssemblyBegin_IS;
    325:   A->ops->assemblyend             = MatAssemblyEnd_IS;
    326:   A->ops->view                    = MatView_IS;
    327:   A->ops->zeroentries             = MatZeroEntries_IS;
    328:   A->ops->setoption               = MatSetOption_IS;

    330:   PetscSplitOwnership(A->comm,&A->m,&A->M);
    331:   PetscSplitOwnership(A->comm,&A->n,&A->N);
    332:   MPI_Scan(&A->m,&b->rend,1,MPIU_INT,MPI_SUM,A->comm);
    333:   b->rstart = b->rend - A->m;

    335:   b->A          = 0;
    336:   b->ctx        = 0;
    337:   b->x          = 0;
    338:   b->y          = 0;
    339:   PetscObjectComposeFunctionDynamic((PetscObject)A,"MatISGetLocalMat_C","MatISGetLocalMat_IS",MatISGetLocalMat_IS);

    341:   return(0);
    342: }