Actual source code: ml.c

      1: #define PETSCKSP_DLL

      3: /*
      4:    Provides an interface to the ML 3.0 smoothed Aggregation
      5: */
     6:  #include src/ksp/pc/pcimpl.h
     7:  #include src/ksp/pc/impls/mg/mgimpl.h
     8:  #include src/mat/impls/aij/seq/aij.h
     9:  #include src/mat/impls/aij/mpi/mpiaij.h

     12: #include <math.h>
     13: #include "ml_include.h"

     16: /* The context (data structure) at each grid level */
     17: typedef struct {
     18:   Vec        x,b,r;           /* global vectors */
     19:   Mat        A,P,R;
     20:   KSP        ksp;
     21: } GridCtx;

     23: /* The context used to input PETSc matrix into ML at fine grid */
     24: typedef struct {
     25:   Mat          A,Aloc;
     26:   Vec          x,y;
     27:   ML_Operator  *mlmat;
     28:   PetscScalar  *pwork; /* tmp array used by PetscML_comm() */
     29: } FineGridCtx;

     31: /* The context associates a ML matrix with a PETSc shell matrix */
     32: typedef struct {
     33:   Mat          A;       /* PETSc shell matrix associated with mlmat */
     34:   ML_Operator  *mlmat;  /* ML matrix assorciated with A */
     35:   Vec          y;
     36: } Mat_MLShell;

     38: /* Private context for the ML preconditioner */
     39: typedef struct {
     40:   ML           *ml_object;
     41:   ML_Aggregate *agg_object;
     42:   GridCtx      *gridctx;
     43:   FineGridCtx  *PetscMLdata;
     44:   PetscInt     fine_level,MaxNlevels,MaxCoarseSize,CoarsenScheme;
     45:   PetscReal    Threshold,DampingFactor;
     46:   PetscTruth   SpectralNormScheme_Anorm;
     47:   PetscMPIInt  size;
     48:   PetscErrorCode (*PCSetUp)(PC);
     49:   PetscErrorCode (*PCDestroy)(PC);
     50: } PC_ML;

     53:    int allocated_space,int columns[],double values[],int row_lengths[]);

     64: /* -------------------------------------------------------------------------- */
     65: /*
     66:    PCSetUp_ML - Prepares for the use of the ML preconditioner
     67:                     by setting data structures and options.

     69:    Input Parameter:
     70: .  pc - the preconditioner context

     72:    Application Interface Routine: PCSetUp()

     74:    Notes:
     75:    The interface routine PCSetUp() is not usually called directly by
     76:    the user, but instead is called by PCApply() if necessary.
     77: */
     81: PetscErrorCode PCSetUp_ML(PC pc)
     82: {
     83:   PetscErrorCode       ierr;
     84:   PetscMPIInt          size;
     85:   FineGridCtx          *PetscMLdata;
     86:   ML                   *ml_object;
     87:   ML_Aggregate         *agg_object;
     88:   ML_Operator          *mlmat;
     89:   PetscInt             nlocal_allcols,Nlevels,mllevel,level,level1,m,fine_level;
     90:   Mat                  A,Aloc;
     91:   GridCtx              *gridctx;
     92:   PC_ML                *pc_ml=PETSC_NULL;
     93:   PetscObjectContainer container;

     96:   PetscObjectQuery((PetscObject)pc,"PC_ML",(PetscObject *)&container);
     97:   if (container) {
     98:     PetscObjectContainerGetPointer(container,(void **)&pc_ml);
     99:   } else {
    100:     SETERRQ(PETSC_ERR_ARG_NULL,"Container does not exit");
    101:   }
    102:
    103:   /* setup special features of PCML */
    104:   /*--------------------------------*/
    105:   /* covert A to Aloc to be used by ML at fine grid */
    106:   A = pc->pmat;
    107:   MPI_Comm_size(A->comm,&size);
    108:   pc_ml->size = size;
    109:   if (size > 1){
    110:     MatConvert_MPIAIJ_ML(A,PETSC_NULL,MAT_INITIAL_MATRIX,&Aloc);
    111:   } else {
    112:     Aloc = A;
    113:   }

    115:   /* create and initialize struct 'PetscMLdata' */
    116:   PetscNew(FineGridCtx,&PetscMLdata);
    117:   PetscMLdata->A    = A;
    118:   PetscMLdata->Aloc = Aloc;
    119:   PetscMalloc((Aloc->n+1)*sizeof(PetscScalar),&PetscMLdata->pwork);
    120:   pc_ml->PetscMLdata = PetscMLdata;

    122:   VecCreate(PETSC_COMM_SELF,&PetscMLdata->x);
    123:   if (size == 1){
    124:     VecSetSizes(PetscMLdata->x,A->n,A->n);
    125:   } else {
    126:     VecSetSizes(PetscMLdata->x,Aloc->n,Aloc->n);
    127:   }
    128:   VecSetType(PetscMLdata->x,VECSEQ);

    130:   VecCreate(PETSC_COMM_SELF,&PetscMLdata->y);
    131:   VecSetSizes(PetscMLdata->y,A->m,PETSC_DECIDE);
    132:   VecSetType(PetscMLdata->y,VECSEQ);
    133:
    134:   /* create ML discretization matrix at fine grid */
    135:   MatGetSize(Aloc,&m,&nlocal_allcols);
    136:   ML_Create(&ml_object,pc_ml->MaxNlevels);
    137:   ML_Init_Amatrix(ml_object,0,m,m,PetscMLdata);
    138:   ML_Set_Amatrix_Getrow(ml_object,0,PetscML_getrow,PetscML_comm,nlocal_allcols);
    139:   ML_Set_Amatrix_Matvec(ml_object,0,PetscML_matvec);

    141:   /* aggregation */
    142:   ML_Aggregate_Create(&agg_object);
    143:   ML_Aggregate_Set_MaxCoarseSize(agg_object,pc_ml->MaxCoarseSize);
    144:   /* set options */
    145:   switch (pc_ml->CoarsenScheme) {
    146:   case 1:
    147:     ML_Aggregate_Set_CoarsenScheme_Coupled(agg_object);break;
    148:   case 2:
    149:     ML_Aggregate_Set_CoarsenScheme_MIS(agg_object);break;
    150:   case 3:
    151:     ML_Aggregate_Set_CoarsenScheme_METIS(agg_object);break;
    152:   }
    153:   ML_Aggregate_Set_Threshold(agg_object,pc_ml->Threshold);
    154:   ML_Aggregate_Set_DampingFactor(agg_object,pc_ml->DampingFactor);
    155:   if (pc_ml->SpectralNormScheme_Anorm){
    156:     ML_Aggregate_Set_SpectralNormScheme_Anorm(agg_object);
    157:   }
    158:
    159:   Nlevels = ML_Gen_MGHierarchy_UsingAggregation(ml_object,0,ML_INCREASING,agg_object);
    160:   if (Nlevels<=0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Nlevels %d must > 0",Nlevels);
    161:   PCMGSetLevels(pc,Nlevels,PETSC_NULL);
    162:   PCSetFromOptions_MG(pc); /* should be called in PCSetFromOptions_ML(), but cannot be called prior to PCMGSetLevels() */
    163:   pc_ml->ml_object  = ml_object;
    164:   pc_ml->agg_object = agg_object;

    166:   PetscMalloc(Nlevels*sizeof(GridCtx),&gridctx);
    167:   fine_level = Nlevels - 1;
    168:   pc_ml->gridctx = gridctx;
    169:   pc_ml->fine_level = fine_level;

    171:   /* wrap ML matrices by PETSc shell matrices at coarsened grids.
    172:      Level 0 is the finest grid for ML, but coarsest for PETSc! */
    173:   gridctx[fine_level].A = A;
    174:   level = fine_level - 1;
    175:   if (size == 1){ /* convert ML P, R and A into seqaij format */
    176:     for (mllevel=1; mllevel<Nlevels; mllevel++){
    177:       mlmat  = &(ml_object->Pmat[mllevel]);
    178:       MatWrapML_SeqAIJ(mlmat,&gridctx[level].P);
    179:       mlmat  = &(ml_object->Amat[mllevel]);
    180:       MatWrapML_SeqAIJ(mlmat,&gridctx[level].A);
    181:       mlmat  = &(ml_object->Rmat[mllevel-1]);
    182:       MatWrapML_SeqAIJ(mlmat,&gridctx[level].R);
    183:       level--;
    184:     }
    185:   } else { /* convert ML P and R into shell format, ML A into mpiaij format */
    186:     for (mllevel=1; mllevel<Nlevels; mllevel++){
    187:       mlmat  = &(ml_object->Pmat[mllevel]);
    188:       MatWrapML_SHELL(mlmat,&gridctx[level].P);
    189:       mlmat  = &(ml_object->Rmat[mllevel-1]);
    190:       MatWrapML_SHELL(mlmat,&gridctx[level].R);
    191:       mlmat  = &(ml_object->Amat[mllevel]);
    192:       MatWrapML_MPIAIJ(mlmat,&gridctx[level].A);
    193:       level--;
    194:     }
    195:   }

    197:   /* create coarse level and the interpolation between the levels */
    198:   for (level=0; level<fine_level; level++){
    199:     VecCreate(gridctx[level].A->comm,&gridctx[level].x);
    200:     VecSetSizes(gridctx[level].x,gridctx[level].A->n,PETSC_DECIDE);
    201:     VecSetType(gridctx[level].x,VECMPI);
    202:     PCMGSetX(pc,level,gridctx[level].x);
    203:
    204:     VecCreate(gridctx[level].A->comm,&gridctx[level].b);
    205:     VecSetSizes(gridctx[level].b,gridctx[level].A->m,PETSC_DECIDE);
    206:     VecSetType(gridctx[level].b,VECMPI);
    207:     PCMGSetRhs(pc,level,gridctx[level].b);
    208:
    209:     level1 = level + 1;
    210:     VecCreate(gridctx[level1].A->comm,&gridctx[level1].r);
    211:     VecSetSizes(gridctx[level1].r,gridctx[level1].A->m,PETSC_DECIDE);
    212:     VecSetType(gridctx[level1].r,VECMPI);
    213:     PCMGSetR(pc,level1,gridctx[level1].r);

    215:     PCMGSetInterpolate(pc,level1,gridctx[level].P);
    216:     PCMGSetRestriction(pc,level1,gridctx[level].R);

    218:     if (level == 0){
    219:       PCMGGetCoarseSolve(pc,&gridctx[level].ksp);
    220:     } else {
    221:       PCMGGetSmoother(pc,level,&gridctx[level].ksp);
    222:       PCMGSetResidual(pc,level,PCMGDefaultResidual,gridctx[level].A);
    223:     }
    224:     KSPSetOperators(gridctx[level].ksp,gridctx[level].A,gridctx[level].A,DIFFERENT_NONZERO_PATTERN);
    225:   }
    226:   PCMGGetSmoother(pc,fine_level,&gridctx[fine_level].ksp);
    227:   PCMGSetResidual(pc,fine_level,PCMGDefaultResidual,gridctx[fine_level].A);
    228:   KSPSetOperators(gridctx[fine_level].ksp,gridctx[level].A,gridctx[fine_level].A,DIFFERENT_NONZERO_PATTERN);
    229:   KSPSetOptionsPrefix(gridctx[fine_level].ksp,"mg_fine_");
    230:
    231:   /* now call PCSetUp_MG()         */
    232:   /*--------------------------------*/
    233:   (*pc_ml->PCSetUp)(pc);
    234:   return(0);
    235: }

    239: PetscErrorCode PetscObjectContainerDestroy_PC_ML(void *ptr)
    240: {
    241:   PetscErrorCode       ierr;
    242:   PC_ML                *pc_ml = (PC_ML*)ptr;
    243:   PetscInt             level;

    246:   if (pc_ml->size > 1){MatDestroy(pc_ml->PetscMLdata->Aloc);}
    247:   ML_Aggregate_Destroy(&pc_ml->agg_object);
    248:   ML_Destroy(&pc_ml->ml_object);

    250:   PetscFree(pc_ml->PetscMLdata->pwork);
    251:   if (pc_ml->PetscMLdata->x){VecDestroy(pc_ml->PetscMLdata->x);}
    252:   if (pc_ml->PetscMLdata->y){VecDestroy(pc_ml->PetscMLdata->y);}
    253:   PetscFree(pc_ml->PetscMLdata);

    255:   for (level=0; level<pc_ml->fine_level; level++){
    256:     if (pc_ml->gridctx[level].A){MatDestroy(pc_ml->gridctx[level].A);}
    257:     if (pc_ml->gridctx[level].P){MatDestroy(pc_ml->gridctx[level].P);}
    258:     if (pc_ml->gridctx[level].R){MatDestroy(pc_ml->gridctx[level].R);}
    259:     if (pc_ml->gridctx[level].x){VecDestroy(pc_ml->gridctx[level].x);}
    260:     if (pc_ml->gridctx[level].b){VecDestroy(pc_ml->gridctx[level].b);}
    261:     if (pc_ml->gridctx[level+1].r){VecDestroy(pc_ml->gridctx[level+1].r);}
    262:   }
    263:   PetscFree(pc_ml->gridctx);
    264:   PetscFree(pc_ml);
    265:   return(0);
    266: }
    267: /* -------------------------------------------------------------------------- */
    268: /*
    269:    PCDestroy_ML - Destroys the private context for the ML preconditioner
    270:    that was created with PCCreate_ML().

    272:    Input Parameter:
    273: .  pc - the preconditioner context

    275:    Application Interface Routine: PCDestroy()
    276: */
    279: PetscErrorCode PCDestroy_ML(PC pc)
    280: {
    281:   PetscErrorCode       ierr;
    282:   PC_ML                *pc_ml=PETSC_NULL;
    283:   PetscObjectContainer container;

    286:   PetscObjectQuery((PetscObject)pc,"PC_ML",(PetscObject *)&container);
    287:   if (container) {
    288:     PetscObjectContainerGetPointer(container,(void **)&pc_ml);
    289:     pc->ops->destroy = pc_ml->PCDestroy;
    290:   } else {
    291:     SETERRQ(PETSC_ERR_ARG_NULL,"Container does not exit");
    292:   }
    293:   /* detach pc and PC_ML and dereference container */
    294:   PetscObjectCompose((PetscObject)pc,"PC_ML",0);
    295:   (*pc->ops->destroy)(pc);

    297:   PetscObjectContainerDestroy(container);
    298:   return(0);
    299: }

    303: PetscErrorCode PCSetFromOptions_ML(PC pc)
    304: {
    305:   PetscErrorCode       ierr;
    306:   PetscInt             indx,m,PrintLevel,MaxNlevels,MaxCoarseSize;
    307:   PetscReal            Threshold,DampingFactor;
    308:   PetscTruth           flg;
    309:   const char           *scheme[] = {"Uncoupled","Coupled","MIS","METIS"};
    310:   PC_ML                *pc_ml=PETSC_NULL;
    311:   PetscObjectContainer container;
    312:   PCMGType             mgtype;

    315:   PetscObjectQuery((PetscObject)pc,"PC_ML",(PetscObject *)&container);
    316:   if (container) {
    317:     PetscObjectContainerGetPointer(container,(void **)&pc_ml);
    318:   } else {
    319:     SETERRQ(PETSC_ERR_ARG_NULL,"Container does not exit");
    320:   }

    322:   /* inherited MG options */
    323:   PetscOptionsHead("Multigrid options(inherited)");
    324:     PetscOptionsInt("-pc_mg_cycles","1 for V cycle, 2 for W-cycle","MGSetCycles",1,&m,&flg);
    325:     PetscOptionsInt("-pc_mg_smoothup","Number of post-smoothing steps","MGSetNumberSmoothUp",1,&m,&flg);
    326:     PetscOptionsInt("-pc_mg_smoothdown","Number of pre-smoothing steps","MGSetNumberSmoothDown",1,&m,&flg);
    327:     PetscOptionsEnum("-pc_mg_type","Multigrid type","PCMGSetType",PCMGTypes,(PetscEnum)PC_MG_MULTIPLICATIVE,(PetscEnum*)&mgtype,&flg);
    328:   PetscOptionsTail();

    330:   /* ML options */
    331:   PetscOptionsHead("ML options");
    332:   /* set defaults */
    333:   PrintLevel    = 0;
    334:   MaxNlevels    = 10;
    335:   MaxCoarseSize = 1;
    336:   indx          = 0;
    337:   Threshold     = 0.0;
    338:   DampingFactor = 4.0/3.0;
    339:
    340:   PetscOptionsInt("-pc_ml_PrintLevel","Print level","ML_Set_PrintLevel",PrintLevel,&PrintLevel,PETSC_NULL);
    341:   ML_Set_PrintLevel(PrintLevel);

    343:   PetscOptionsInt("-pc_ml_maxNlevels","Maximum number of levels","None",MaxNlevels,&MaxNlevels,PETSC_NULL);
    344:   pc_ml->MaxNlevels = MaxNlevels;

    346:   PetscOptionsInt("-pc_ml_maxCoarseSize","Maximum coarsest mesh size","ML_Aggregate_Set_MaxCoarseSize",MaxCoarseSize,&MaxCoarseSize,PETSC_NULL);
    347:   pc_ml->MaxCoarseSize = MaxCoarseSize;

    349:   PetscOptionsEList("-pc_ml_CoarsenScheme","Aggregate Coarsen Scheme","ML_Aggregate_Set_CoarsenScheme_*",scheme,4,scheme[0],&indx,PETSC_NULL);
    350:   pc_ml->CoarsenScheme = indx;

    352:   PetscOptionsReal("-pc_ml_DampingFactor","P damping factor","ML_Aggregate_Set_DampingFactor",DampingFactor,&DampingFactor,PETSC_NULL);
    353:   pc_ml->DampingFactor = DampingFactor;
    354:
    355:   PetscOptionsReal("-pc_ml_Threshold","Smoother drop tol","ML_Aggregate_Set_Threshold",Threshold,&Threshold,PETSC_NULL);
    356:   pc_ml->Threshold = Threshold;

    358:   PetscOptionsTruth("-pc_ml_SpectralNormScheme_Anorm","Method used for estimating spectral radius","ML_Aggregate_Set_SpectralNormScheme_Anorm",PETSC_FALSE,&pc_ml->SpectralNormScheme_Anorm,PETSC_FALSE);
    359:
    360:   PetscOptionsTail();
    361:   return(0);
    362: }

    364: /* -------------------------------------------------------------------------- */
    365: /*
    366:    PCCreate_ML - Creates a ML preconditioner context, PC_ML,
    367:    and sets this as the private data within the generic preconditioning
    368:    context, PC, that was created within PCCreate().

    370:    Input Parameter:
    371: .  pc - the preconditioner context

    373:    Application Interface Routine: PCCreate()
    374: */

    376: /*MC
    377:      PCML - Use geometric multigrid preconditioning. This preconditioner requires you provide
    378:        fine grid discretization matrix. The coarser grid matrices and restriction/interpolation
    379:        operators are computed by ML, with the matrices coverted to PETSc matrices in aij format
    380:        and the restriction/interpolation operators wrapped as PETSc shell matrices.

    382:    Options Database Key:
    383:    Multigrid options(inherited)
    384: +  -pc_mg_cycles <1>: 1 for V cycle, 2 for W-cycle (MGSetCycles)
    385: .  -pc_mg_smoothup <1>: Number of post-smoothing steps (MGSetNumberSmoothUp)
    386: .  -pc_mg_smoothdown <1>: Number of pre-smoothing steps (MGSetNumberSmoothDown)
    387: -  -pc_mg_type <multiplicative> (one of) additive multiplicative full cascade kascade
    388:
    389:    ML options
    390: +  -pc_ml_PrintLevel <0>: Print level (ML_Set_PrintLevel)
    391: .  -pc_ml_maxNlevels <10>: Maximum number of levels (None)
    392: .  -pc_ml_maxCoarseSize <1>: Maximum coarsest mesh size (ML_Aggregate_Set_MaxCoarseSize)
    393: .  -pc_ml_CoarsenScheme <Uncoupled> (one of) Uncoupled Coupled MIS METIS
    394: .  -pc_ml_DampingFactor <1.33333>: P damping factor (ML_Aggregate_Set_DampingFactor)
    395: .  -pc_ml_Threshold <0>: Smoother drop tol (ML_Aggregate_Set_Threshold)
    396: -  -pc_ml_SpectralNormScheme_Anorm: <false> Method used for estimating spectral radius (ML_Aggregate_Set_SpectralNormScheme_Anorm)

    398:    Level: intermediate

    400:   Concepts: multigrid
    401:
    402: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType,
    403:            PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), MPSetCycles(), PCMGSetNumberSmoothDown(),
    404:            PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(),
    405:            PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(),
    406:            PCMGSetCyclesOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR()
    407: M*/

    412: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_ML(PC pc)
    413: {
    414:   PetscErrorCode       ierr;
    415:   PC_ML                *pc_ml;
    416:   PetscObjectContainer container;

    419:   /* initialize pc as PCMG */
    420:   PCSetType(pc,PCMG); /* calls PCCreate_MG() and MGCreate_Private() */

    422:   /* create a supporting struct and attach it to pc */
    423:   PetscNew(PC_ML,&pc_ml);
    424:   PetscObjectContainerCreate(PETSC_COMM_SELF,&container);
    425:   PetscObjectContainerSetPointer(container,pc_ml);
    426:   PetscObjectContainerSetUserDestroy(container,PetscObjectContainerDestroy_PC_ML);
    427:   PetscObjectCompose((PetscObject)pc,"PC_ML",(PetscObject)container);
    428:
    429:   pc_ml->PCSetUp   = pc->ops->setup;
    430:   pc_ml->PCDestroy = pc->ops->destroy;

    432:   /* overwrite the pointers of PCMG by the functions of PCML */
    433:   pc->ops->setfromoptions = PCSetFromOptions_ML;
    434:   pc->ops->setup          = PCSetUp_ML;
    435:   pc->ops->destroy        = PCDestroy_ML;
    436:   return(0);
    437: }

    440: int PetscML_getrow(ML_Operator *ML_data, int N_requested_rows, int requested_rows[],
    441:    int allocated_space, int columns[], double values[], int row_lengths[])
    442: {
    444:   Mat            Aloc;
    445:   Mat_SeqAIJ     *a;
    446:   PetscInt       m,i,j,k=0,row,*aj;
    447:   PetscScalar    *aa;
    448:   FineGridCtx    *ml=(FineGridCtx*)ML_Get_MyGetrowData(ML_data);

    450:   Aloc = ml->Aloc;
    451:   a    = (Mat_SeqAIJ*)Aloc->data;
    452:   MatGetSize(Aloc,&m,PETSC_NULL);

    454:   for (i = 0; i<N_requested_rows; i++) {
    455:     row   = requested_rows[i];
    456:     row_lengths[i] = a->ilen[row];
    457:     if (allocated_space < k+row_lengths[i]) return(0);
    458:     if ( (row >= 0) || (row <= (m-1)) ) {
    459:       aj = a->j + a->i[row];
    460:       aa = a->a + a->i[row];
    461:       for (j=0; j<row_lengths[i]; j++){
    462:         columns[k]  = aj[j];
    463:         values[k++] = aa[j];
    464:       }
    465:     }
    466:   }
    467:   return(1);
    468: }

    470: int PetscML_matvec(ML_Operator *ML_data,int in_length,double p[],int out_length,double ap[])
    471: {
    473:   FineGridCtx    *ml=(FineGridCtx*)ML_Get_MyMatvecData(ML_data);
    474:   Mat            A=ml->A, Aloc=ml->Aloc;
    475:   PetscMPIInt    size;
    476:   PetscScalar    *pwork=ml->pwork;
    477:   PetscInt       i;

    479:   MPI_Comm_size(A->comm,&size);
    480:   if (size == 1){
    481:     VecPlaceArray(ml->x,p);
    482:   } else {
    483:     for (i=0; i<in_length; i++) pwork[i] = p[i];
    484:     PetscML_comm(pwork,ml);
    485:     VecPlaceArray(ml->x,pwork);
    486:   }
    487:   VecPlaceArray(ml->y,ap);
    488:   MatMult(Aloc,ml->x,ml->y);
    489:   return 0;
    490: }

    492: int PetscML_comm(double p[],void *ML_data)
    493: {
    495:   FineGridCtx    *ml=(FineGridCtx*)ML_data;
    496:   Mat            A=ml->A;
    497:   Mat_MPIAIJ     *a = (Mat_MPIAIJ*)A->data;
    498:   PetscMPIInt    size;
    499:   PetscInt       i,in_length=A->m,out_length=ml->Aloc->n;
    500:   PetscScalar    *array;

    502:   MPI_Comm_size(A->comm,&size);
    503:   if (size == 1) return 0;
    504:
    505:   VecPlaceArray(ml->y,p);
    506:   VecScatterBegin(ml->y,a->lvec,INSERT_VALUES,SCATTER_FORWARD,a->Mvctx);
    507:   VecScatterEnd(ml->y,a->lvec,INSERT_VALUES,SCATTER_FORWARD,a->Mvctx);
    508:   VecGetArray(a->lvec,&array);
    509:   for (i=in_length; i<out_length; i++){
    510:     p[i] = array[i-in_length];
    511:   }
    512:   return 0;
    513: }
    516: PetscErrorCode MatMult_ML(Mat A,Vec x,Vec y)
    517: {
    518:   PetscErrorCode   ierr;
    519:   Mat_MLShell      *shell;
    520:   PetscScalar      *xarray,*yarray;
    521:   PetscInt         x_length,y_length;
    522:
    524:   MatShellGetContext(A,(void *)&shell);
    525:   VecGetArray(x,&xarray);
    526:   VecGetArray(y,&yarray);
    527:   x_length = shell->mlmat->invec_leng;
    528:   y_length = shell->mlmat->outvec_leng;

    530:   ML_Operator_Apply(shell->mlmat,x_length,xarray,y_length,yarray);

    532:   VecRestoreArray(x,&xarray);
    533:   VecRestoreArray(y,&yarray);
    534:   return(0);
    535: }
    536: /* MatMultAdd_ML -  Compute y = w + A*x */
    539: PetscErrorCode MatMultAdd_ML(Mat A,Vec x,Vec w,Vec y)
    540: {
    541:   PetscErrorCode    ierr;
    542:   Mat_MLShell       *shell;
    543:   PetscScalar       *xarray,*yarray;
    544:   const PetscScalar one=1.0;
    545:   PetscInt          x_length,y_length;
    546:
    548:   MatShellGetContext(A,(void *)&shell);
    549:   VecGetArray(x,&xarray);
    550:   VecGetArray(y,&yarray);

    552:   x_length = shell->mlmat->invec_leng;
    553:   y_length = shell->mlmat->outvec_leng;

    555:   ML_Operator_Apply(shell->mlmat,x_length,xarray,y_length,yarray);

    557:   VecRestoreArray(x,&xarray);
    558:   VecRestoreArray(y,&yarray);
    559:   VecAXPY(y,one,w);

    561:   return(0);
    562: }

    564: /* newtype is ignored because "ml" is not listed under Petsc MatType yet */
    567: PetscErrorCode MatConvert_MPIAIJ_ML(Mat A,MatType newtype,MatReuse scall,Mat *Aloc)
    568: {
    569:   PetscErrorCode  ierr;
    570:   Mat_MPIAIJ      *mpimat=(Mat_MPIAIJ*)A->data;
    571:   Mat_SeqAIJ      *mat,*a=(Mat_SeqAIJ*)(mpimat->A)->data,*b=(Mat_SeqAIJ*)(mpimat->B)->data;
    572:   PetscInt        *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j;
    573:   PetscScalar     *aa=a->a,*ba=b->a,*ca;
    574:   PetscInt        am=A->m,an=A->n,i,j,k;
    575:   PetscInt        *ci,*cj,ncols;

    578:   if (am != an) SETERRQ2(PETSC_ERR_ARG_WRONG,"A must have a square diagonal portion, am: %d != an: %d",am,an);

    580:   if (scall == MAT_INITIAL_MATRIX){
    581:     PetscMalloc((1+am)*sizeof(PetscInt),&ci);
    582:     ci[0] = 0;
    583:     for (i=0; i<am; i++){
    584:       ci[i+1] = ci[i] + (ai[i+1] - ai[i]) + (bi[i+1] - bi[i]);
    585:     }
    586:     PetscMalloc((1+ci[am])*sizeof(PetscInt),&cj);
    587:     PetscMalloc((1+ci[am])*sizeof(PetscScalar),&ca);

    589:     k = 0;
    590:     for (i=0; i<am; i++){
    591:       /* diagonal portion of A */
    592:       ncols = ai[i+1] - ai[i];
    593:       for (j=0; j<ncols; j++) {
    594:         cj[k]   = *aj++;
    595:         ca[k++] = *aa++;
    596:       }
    597:       /* off-diagonal portion of A */
    598:       ncols = bi[i+1] - bi[i];
    599:       for (j=0; j<ncols; j++) {
    600:         cj[k]   = an + (*bj); bj++;
    601:         ca[k++] = *ba++;
    602:       }
    603:     }
    604:     if (k != ci[am]) SETERRQ2(PETSC_ERR_ARG_WRONG,"k: %d != ci[am]: %d",k,ci[am]);

    606:     /* put together the new matrix */
    607:     an = mpimat->A->n+mpimat->B->n;
    608:     MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,am,an,ci,cj,ca,Aloc);

    610:     /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
    611:     /* Since these are PETSc arrays, change flags to free them as necessary. */
    612:     mat = (Mat_SeqAIJ*)(*Aloc)->data;
    613:     mat->freedata = PETSC_TRUE;
    614:     mat->nonew    = 0;
    615:   } else if (scall == MAT_REUSE_MATRIX){
    616:     mat=(Mat_SeqAIJ*)(*Aloc)->data;
    617:     ci = mat->i; cj = mat->j; ca = mat->a;
    618:     for (i=0; i<am; i++) {
    619:       /* diagonal portion of A */
    620:       ncols = ai[i+1] - ai[i];
    621:       for (j=0; j<ncols; j++) *ca++ = *aa++;
    622:       /* off-diagonal portion of A */
    623:       ncols = bi[i+1] - bi[i];
    624:       for (j=0; j<ncols; j++) *ca++ = *ba++;
    625:     }
    626:   } else {
    627:     SETERRQ1(PETSC_ERR_ARG_WRONG,"Invalid MatReuse %d",(int)scall);
    628:   }
    629:   return(0);
    630: }
    634: PetscErrorCode MatDestroy_ML(Mat A)
    635: {
    637:   Mat_MLShell    *shell;

    640:   MatShellGetContext(A,(void *)&shell);
    641:   VecDestroy(shell->y);
    642:   PetscFree(shell);
    643:   MatDestroy_Shell(A);
    644:   return(0);
    645: }

    649: PetscErrorCode MatWrapML_SeqAIJ(ML_Operator *mlmat,Mat *newmat)
    650: {
    651:   struct ML_CSR_MSRdata *matdata = (struct ML_CSR_MSRdata *)mlmat->data;
    652:   PetscErrorCode        ierr;
    653:   PetscInt              m=mlmat->outvec_leng,n=mlmat->invec_leng,*nnz,nz_max;
    654:   PetscInt              *ml_cols=matdata->columns,*aj,i,j,k;
    655:   PetscScalar           *ml_vals=matdata->values,*aa;
    656:
    658:   if ( mlmat->getrow == NULL) SETERRQ(PETSC_ERR_ARG_NULL,"mlmat->getrow = NULL");
    659:   if (m != n){ /* ML Pmat and Rmat are in CSR format. Pass array pointers into SeqAIJ matrix */
    660:     MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,m,n,matdata->rowptr,ml_cols,ml_vals,newmat);
    661:     return(0);
    662:   }

    664:   /* ML Amat is in MSR format. Copy its data into SeqAIJ matrix */
    665:   MatCreate(PETSC_COMM_SELF,newmat);
    666:   MatSetSizes(*newmat,m,n,PETSC_DECIDE,PETSC_DECIDE);
    667:   MatSetType(*newmat,MATSEQAIJ);
    668:   PetscMalloc((m+1)*sizeof(PetscInt),&nnz);

    670:   nz_max = 0;
    671:   for (i=0; i<m; i++) {
    672:     nnz[i] = ml_cols[i+1] - ml_cols[i] + 1;
    673:     if (nnz[i] > nz_max) nz_max = nnz[i];
    674:   }
    675:   MatSeqAIJSetPreallocation(*newmat,0,nnz);
    676:   MatSetOption(*newmat,MAT_COLUMNS_SORTED); /* check! */

    678:   nz_max++;
    679:   PetscMalloc(nz_max*(sizeof(PetscInt)+sizeof(PetscScalar)),&aj);
    680:   aa = (PetscScalar*)(aj + nz_max);

    682:   for (i=0; i<m; i++){
    683:     k = 0;
    684:     /* diagonal entry */
    685:     aj[k] = i; aa[k++] = ml_vals[i];
    686:     /* off diagonal entries */
    687:     for (j=ml_cols[i]; j<ml_cols[i+1]; j++){
    688:       aj[k] = ml_cols[j]; aa[k++] = ml_vals[j];
    689:     }
    690:     /* sort aj and aa */
    691:     PetscSortIntWithScalarArray(nnz[i],aj,aa);
    692:     MatSetValues(*newmat,1,&i,nnz[i],aj,aa,INSERT_VALUES);
    693:   }
    694:   MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
    695:   MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
    696:   PetscFree(aj);
    697:   PetscFree(nnz);
    698:   return(0);
    699: }

    703: PetscErrorCode MatWrapML_SHELL(ML_Operator *mlmat,Mat *newmat)
    704: {
    706:   PetscInt       m,n;
    707:   ML_Comm        *MLcomm;
    708:   Mat_MLShell    *shellctx;

    711:   m = mlmat->outvec_leng;
    712:   n = mlmat->invec_leng;
    713:   if (!m || !n){
    714:     newmat = PETSC_NULL;
    715:   } else {
    716:     MLcomm = mlmat->comm;
    717:     PetscNew(Mat_MLShell,&shellctx);
    718:     MatCreateShell(MLcomm->USR_comm,m,n,PETSC_DETERMINE,PETSC_DETERMINE,shellctx,newmat);
    719:     MatShellSetOperation(*newmat,MATOP_MULT,(void(*)(void))MatMult_ML);
    720:     MatShellSetOperation(*newmat,MATOP_MULT_ADD,(void(*)(void))MatMultAdd_ML);
    721:     shellctx->A         = *newmat;
    722:     shellctx->mlmat     = mlmat;
    723:     VecCreate(PETSC_COMM_WORLD,&shellctx->y);
    724:     VecSetSizes(shellctx->y,m,PETSC_DECIDE);
    725:     VecSetFromOptions(shellctx->y);
    726:     (*newmat)->ops->destroy = MatDestroy_ML;
    727:   }
    728:   return(0);
    729: }

    733: PetscErrorCode MatWrapML_MPIAIJ(ML_Operator *mlmat,Mat *newmat)
    734: {
    735:   struct ML_CSR_MSRdata *matdata = (struct ML_CSR_MSRdata *)mlmat->data;
    736:   PetscInt              *ml_cols=matdata->columns,*aj;
    737:   PetscScalar           *ml_vals=matdata->values,*aa;
    738:   PetscErrorCode        ierr;
    739:   PetscInt              i,j,k,*gordering;
    740:   PetscInt              m=mlmat->outvec_leng,n,*nnzA,*nnzB,*nnz,nz_max,row;
    741:   Mat                   A;

    744:   if (mlmat->getrow == NULL) SETERRQ(PETSC_ERR_ARG_NULL,"mlmat->getrow = NULL");
    745:   n = mlmat->invec_leng;
    746:   if (m != n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"m %d must equal to n %d",m,n);

    748:   MatCreate(mlmat->comm->USR_comm,&A);
    749:   MatSetSizes(A,m,n,PETSC_DECIDE,PETSC_DECIDE);
    750:   MatSetType(A,MATMPIAIJ);
    751:   PetscMalloc3(m,PetscInt,&nnzA,m,PetscInt,&nnzB,m,PetscInt,&nnz);
    752:
    753:   nz_max = 0;
    754:   for (i=0; i<m; i++){
    755:     nnz[i] = ml_cols[i+1] - ml_cols[i] + 1;
    756:     if (nz_max < nnz[i]) nz_max = nnz[i];
    757:     nnzA[i] = 1; /* diag */
    758:     for (j=ml_cols[i]; j<ml_cols[i+1]; j++){
    759:       if (ml_cols[j] < m) nnzA[i]++;
    760:     }
    761:     nnzB[i] = nnz[i] - nnzA[i];
    762:   }
    763:   MatMPIAIJSetPreallocation(A,0,nnzA,0,nnzB);

    765:   /* insert mat values -- remap row and column indices */
    766:   nz_max++;
    767:   PetscMalloc(nz_max*(sizeof(PetscInt)+sizeof(PetscScalar)),&aj);
    768:   aa = (PetscScalar*)(aj + nz_max);
    769:   ML_build_global_numbering(mlmat,mlmat->comm,&gordering);
    770:   for (i=0; i<m; i++){
    771:     row = gordering[i];
    772:     k = 0;
    773:     /* diagonal entry */
    774:     aj[k] = row; aa[k++] = ml_vals[i];
    775:     /* off diagonal entries */
    776:     for (j=ml_cols[i]; j<ml_cols[i+1]; j++){
    777:       aj[k] = gordering[ml_cols[j]]; aa[k++] = ml_vals[j];
    778:     }
    779:     MatSetValues(A,1,&row,nnz[i],aj,aa,INSERT_VALUES);
    780:   }
    781:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    782:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    783:   *newmat = A;

    785:   PetscFree3(nnzA,nnzB,nnz);
    786:   PetscFree(aj);
    787:   return(0);
    788: }