Actual source code: packm.c

      2: #include <../src/dm/impls/composite/packimpl.h>

      4: static PetscErrorCode DMCreateMatrix_Composite_Nest(DM dm,Mat *J)
      5: {
      6:   const DM_Composite           *com = (DM_Composite*)dm->data;
      7:   const struct DMCompositeLink *rlink,*clink;
      8:   PetscErrorCode               ierr;
      9:   IS                           *isg;
     10:   Mat                          *submats;
     11:   PetscInt                     i,j,n;

     14:   n = com->nDM;                 /* Total number of entries */

     16:   /* Explicit index sets are not required for MatCreateNest, but getting them here allows MatNest to do consistency
     17:    * checking and allows ISEqual to compare by identity instead of by contents. */
     18:   DMCompositeGetGlobalISs(dm,&isg);

     20:   /* Get submatrices */
     21:   PetscMalloc1(n*n,&submats);
     22:   for (i=0,rlink=com->next; rlink; i++,rlink=rlink->next) {
     23:     for (j=0,clink=com->next; clink; j++,clink=clink->next) {
     24:       Mat sub = NULL;
     25:       if (i == j) {
     26:         DMCreateMatrix(rlink->dm,&sub);
     27:       } else if (com->FormCoupleLocations) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot manage off-diagonal parts yet");
     28:       submats[i*n+j] = sub;
     29:     }
     30:   }

     32:   MatCreateNest(PetscObjectComm((PetscObject)dm),n,isg,n,isg,submats,J);

     34:   /* Disown references */
     35:   for (i=0; i<n; i++) {ISDestroy(&isg[i]);}
     36:   PetscFree(isg);

     38:   for (i=0; i<n*n; i++) {
     39:     if (submats[i]) {MatDestroy(&submats[i]);}
     40:   }
     41:   PetscFree(submats);
     42:   return(0);
     43: }

     45: static PetscErrorCode DMCreateMatrix_Composite_AIJ(DM dm,Mat *J)
     46: {
     47:   PetscErrorCode         ierr;
     48:   DM_Composite           *com = (DM_Composite*)dm->data;
     49:   struct DMCompositeLink *next;
     50:   PetscInt               m,*dnz,*onz,i,j,mA;
     51:   Mat                    Atmp;
     52:   PetscMPIInt            rank;
     53:   PetscBool              dense = PETSC_FALSE;

     56:   /* use global vector to determine layout needed for matrix */
     57:   m = com->n;

     59:   MatCreate(PetscObjectComm((PetscObject)dm),J);
     60:   MatSetSizes(*J,m,m,PETSC_DETERMINE,PETSC_DETERMINE);
     61:   MatSetType(*J,dm->mattype);

     63:   /*
     64:      Extremely inefficient but will compute entire Jacobian for testing
     65:   */
     66:   PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dmcomposite_dense_jacobian",&dense,NULL);
     67:   if (dense) {
     68:     PetscInt    rstart,rend,*indices;
     69:     PetscScalar *values;

     71:     mA   = com->N;
     72:     MatMPIAIJSetPreallocation(*J,mA,NULL,mA-m,NULL);
     73:     MatSeqAIJSetPreallocation(*J,mA,NULL);

     75:     MatGetOwnershipRange(*J,&rstart,&rend);
     76:     PetscMalloc2(mA,&values,mA,&indices);
     77:     PetscArrayzero(values,mA);
     78:     for (i=0; i<mA; i++) indices[i] = i;
     79:     for (i=rstart; i<rend; i++) {
     80:       MatSetValues(*J,1,&i,mA,indices,values,INSERT_VALUES);
     81:     }
     82:     PetscFree2(values,indices);
     83:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
     84:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
     85:     return(0);
     86:   }

     88:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
     89:   MatPreallocateInitialize(PetscObjectComm((PetscObject)dm),m,m,dnz,onz);
     90:   /* loop over packed objects, handling one at at time */
     91:   next = com->next;
     92:   while (next) {
     93:     PetscInt       nc,rstart,*ccols,maxnc;
     94:     const PetscInt *cols,*rstarts;
     95:     PetscMPIInt    proc;

     97:     DMCreateMatrix(next->dm,&Atmp);
     98:     MatGetOwnershipRange(Atmp,&rstart,NULL);
     99:     MatGetOwnershipRanges(Atmp,&rstarts);
    100:     MatGetLocalSize(Atmp,&mA,NULL);

    102:     maxnc = 0;
    103:     for (i=0; i<mA; i++) {
    104:       MatGetRow(Atmp,rstart+i,&nc,NULL,NULL);
    105:       maxnc = PetscMax(nc,maxnc);
    106:       MatRestoreRow(Atmp,rstart+i,&nc,NULL,NULL);
    107:     }
    108:     PetscMalloc1(maxnc,&ccols);
    109:     for (i=0; i<mA; i++) {
    110:       MatGetRow(Atmp,rstart+i,&nc,&cols,NULL);
    111:       /* remap the columns taking into how much they are shifted on each process */
    112:       for (j=0; j<nc; j++) {
    113:         proc = 0;
    114:         while (cols[j] >= rstarts[proc+1]) proc++;
    115:         ccols[j] = cols[j] + next->grstarts[proc] - rstarts[proc];
    116:       }
    117:       MatPreallocateSet(com->rstart+next->rstart+i,nc,ccols,dnz,onz);
    118:       MatRestoreRow(Atmp,rstart+i,&nc,&cols,NULL);
    119:     }
    120:     PetscFree(ccols);
    121:     MatDestroy(&Atmp);
    122:     next = next->next;
    123:   }
    124:   if (com->FormCoupleLocations) {
    125:     (*com->FormCoupleLocations)(dm,NULL,dnz,onz,__rstart,__nrows,__start,__end);
    126:   }
    127:   MatMPIAIJSetPreallocation(*J,0,dnz,0,onz);
    128:   MatSeqAIJSetPreallocation(*J,0,dnz);
    129:   MatPreallocateFinalize(dnz,onz);

    131:   if (dm->prealloc_only) return(0);

    133:   next = com->next;
    134:   while (next) {
    135:     PetscInt          nc,rstart,row,maxnc,*ccols;
    136:     const PetscInt    *cols,*rstarts;
    137:     const PetscScalar *values;
    138:     PetscMPIInt       proc;

    140:     DMCreateMatrix(next->dm,&Atmp);
    141:     MatGetOwnershipRange(Atmp,&rstart,NULL);
    142:     MatGetOwnershipRanges(Atmp,&rstarts);
    143:     MatGetLocalSize(Atmp,&mA,NULL);
    144:     maxnc = 0;
    145:     for (i=0; i<mA; i++) {
    146:       MatGetRow(Atmp,rstart+i,&nc,NULL,NULL);
    147:       maxnc = PetscMax(nc,maxnc);
    148:       MatRestoreRow(Atmp,rstart+i,&nc,NULL,NULL);
    149:     }
    150:     PetscMalloc1(maxnc,&ccols);
    151:     for (i=0; i<mA; i++) {
    152:       MatGetRow(Atmp,rstart+i,&nc,(const PetscInt**)&cols,&values);
    153:       for (j=0; j<nc; j++) {
    154:         proc = 0;
    155:         while (cols[j] >= rstarts[proc+1]) proc++;
    156:         ccols[j] = cols[j] + next->grstarts[proc] - rstarts[proc];
    157:       }
    158:       row  = com->rstart+next->rstart+i;
    159:       MatSetValues(*J,1,&row,nc,ccols,values,INSERT_VALUES);
    160:       MatRestoreRow(Atmp,rstart+i,&nc,(const PetscInt**)&cols,&values);
    161:     }
    162:     PetscFree(ccols);
    163:     MatDestroy(&Atmp);
    164:     next = next->next;
    165:   }
    166:   if (com->FormCoupleLocations) {
    167:     PetscInt __rstart;
    168:     MatGetOwnershipRange(*J,&__rstart,NULL);
    169:     (*com->FormCoupleLocations)(dm,*J,NULL,NULL,__rstart,0,0,0);
    170:   }
    171:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
    172:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
    173:   return(0);
    174: }

    176: PetscErrorCode DMCreateMatrix_Composite(DM dm,Mat *J)
    177: {
    178:   PetscErrorCode         ierr;
    179:   PetscBool              usenest;
    180:   ISLocalToGlobalMapping ltogmap;

    183:   DMSetFromOptions(dm);
    184:   DMSetUp(dm);
    185:   PetscStrcmp(dm->mattype,MATNEST,&usenest);
    186:   if (usenest) {
    187:     DMCreateMatrix_Composite_Nest(dm,J);
    188:   } else {
    189:     DMCreateMatrix_Composite_AIJ(dm,J);
    190:   }

    192:   DMGetLocalToGlobalMapping(dm,&ltogmap);
    193:   MatSetLocalToGlobalMapping(*J,ltogmap,ltogmap);
    194:   MatSetDM(*J,dm);
    195:   return(0);
    196: }