Actual source code: matrart.c

      2: /*
      3:   Defines projective product routines where A is a SeqAIJ matrix
      4:           C = R * A * R^T
      5: */

      7: #include <../src/mat/impls/aij/seq/aij.h>
      8: #include <../src/mat/utils/freespace.h>
      9: #include <../src/mat/impls/dense/seq/dense.h>

     11: PetscErrorCode MatDestroy_SeqAIJ_RARt(void *data)
     12: {
     14:   Mat_RARt       *rart = (Mat_RARt*)data;

     17:   MatTransposeColoringDestroy(&rart->matcoloring);
     18:   MatDestroy(&rart->Rt);
     19:   MatDestroy(&rart->RARt);
     20:   MatDestroy(&rart->ARt);
     21:   PetscFree(rart->work);
     22:   if (rart->destroy) {
     23:     (*rart->destroy)(rart->data);
     24:   }
     25:   PetscFree(rart);
     26:   return(0);
     27: }

     29: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,PetscReal fill,Mat C)
     30: {
     31:   PetscErrorCode       ierr;
     32:   Mat                  P;
     33:   PetscInt             *rti,*rtj;
     34:   Mat_RARt             *rart;
     35:   MatColoring          coloring;
     36:   MatTransposeColoring matcoloring;
     37:   ISColoring           iscoloring;
     38:   Mat                  Rt_dense,RARt_dense;

     41:   MatCheckProduct(C,4);
     42:   if (C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data not empty");
     43:   /* create symbolic P=Rt */
     44:   MatGetSymbolicTranspose_SeqAIJ(R,&rti,&rtj);
     45:   MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,R->cmap->n,R->rmap->n,rti,rtj,NULL,&P);

     47:   /* get symbolic C=Pt*A*P */
     48:   MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy(A,P,fill,C);
     49:   MatSetBlockSizes(C,PetscAbs(R->rmap->bs),PetscAbs(R->rmap->bs));
     50:   C->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart;

     52:   /* create a supporting struct */
     53:   PetscNew(&rart);
     54:   C->product->data    = rart;
     55:   C->product->destroy = MatDestroy_SeqAIJ_RARt;

     57:   /* ------ Use coloring ---------- */
     58:   /* inode causes memory problem */
     59:   MatSetOption(C,MAT_USE_INODES,PETSC_FALSE);

     61:   /* Create MatTransposeColoring from symbolic C=R*A*R^T */
     62:   MatColoringCreate(C,&coloring);
     63:   MatColoringSetDistance(coloring,2);
     64:   MatColoringSetType(coloring,MATCOLORINGSL);
     65:   MatColoringSetFromOptions(coloring);
     66:   MatColoringApply(coloring,&iscoloring);
     67:   MatColoringDestroy(&coloring);
     68:   MatTransposeColoringCreate(C,iscoloring,&matcoloring);

     70:   rart->matcoloring = matcoloring;
     71:   ISColoringDestroy(&iscoloring);

     73:   /* Create Rt_dense */
     74:   MatCreate(PETSC_COMM_SELF,&Rt_dense);
     75:   MatSetSizes(Rt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);
     76:   MatSetType(Rt_dense,MATSEQDENSE);
     77:   MatSeqDenseSetPreallocation(Rt_dense,NULL);

     79:   Rt_dense->assembled = PETSC_TRUE;
     80:   rart->Rt            = Rt_dense;

     82:   /* Create RARt_dense = R*A*Rt_dense */
     83:   MatCreate(PETSC_COMM_SELF,&RARt_dense);
     84:   MatSetSizes(RARt_dense,C->rmap->n,matcoloring->ncolors,C->rmap->n,matcoloring->ncolors);
     85:   MatSetType(RARt_dense,MATSEQDENSE);
     86:   MatSeqDenseSetPreallocation(RARt_dense,NULL);

     88:   rart->RARt = RARt_dense;

     90:   /* Allocate work array to store columns of A*R^T used in MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense() */
     91:   PetscMalloc1(A->rmap->n*4,&rart->work);

     93:   /* clean up */
     94:   MatRestoreSymbolicTranspose_SeqAIJ(R,&rti,&rtj);
     95:   MatDestroy(&P);

     97: #if defined(PETSC_USE_INFO)
     98:   {
     99:     Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;
    100:     PetscReal density = (PetscReal)(c->nz)/(RARt_dense->rmap->n*RARt_dense->cmap->n);
    101:     PetscInfo(C,"C=R*(A*Rt) via coloring C - use sparse-dense inner products\n");
    102:     PetscInfo6(C,"RARt_den %D %D; Rt %D %D (RARt->nz %D)/(m*ncolors)=%g\n",RARt_dense->rmap->n,RARt_dense->cmap->n,R->cmap->n,R->rmap->n,c->nz,density);
    103:   }
    104: #endif
    105:   return(0);
    106: }

    108: /*
    109:  RAB = R * A * B, R and A in seqaij format, B in dense format;
    110: */
    111: PetscErrorCode MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(Mat R,Mat A,Mat B,Mat RAB,PetscScalar *work)
    112: {
    113:   Mat_SeqAIJ        *a=(Mat_SeqAIJ*)A->data,*r=(Mat_SeqAIJ*)R->data;
    114:   PetscErrorCode    ierr;
    115:   PetscScalar       r1,r2,r3,r4;
    116:   const PetscScalar *b,*b1,*b2,*b3,*b4;
    117:   MatScalar         *aa,*ra;
    118:   PetscInt          cn =B->cmap->n,bm=B->rmap->n,col,i,j,n,*ai=a->i,*aj,am=A->rmap->n;
    119:   PetscInt          am2=2*am,am3=3*am,bm4=4*bm;
    120:   PetscScalar       *d,*c,*c2,*c3,*c4;
    121:   PetscInt          *rj,rm=R->rmap->n,dm=RAB->rmap->n,dn=RAB->cmap->n;
    122:   PetscInt         rm2=2*rm,rm3=3*rm,colrm;

    125:   if (!dm || !dn) return(0);
    126:   if (bm != A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number columns in A %D not equal rows in B %D\n",A->cmap->n,bm);
    127:   if (am != R->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number columns in R %D not equal rows in A %D\n",R->cmap->n,am);
    128:   if (R->rmap->n != RAB->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows in RAB %D not equal rows in R %D\n",RAB->rmap->n,R->rmap->n);
    129:   if (B->cmap->n != RAB->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number columns in RAB %D not equal columns in B %D\n",RAB->cmap->n,B->cmap->n);

    131:   { /*
    132:      This approach is not as good as original ones (will be removed later), but it reveals that
    133:      AB_den=A*B takes almost all execution time in R*A*B for src/ksp/ksp/tutorials/ex56.c
    134:      */
    135:     PetscBool via_matmatmult=PETSC_FALSE;
    136:     PetscOptionsGetBool(NULL,NULL,"-matrart_via_matmatmult",&via_matmatmult,NULL);
    137:     if (via_matmatmult) {
    138:       Mat AB_den = NULL;
    139:       MatCreate(PetscObjectComm((PetscObject)A),&AB_den);
    140:       MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,0.0,AB_den);
    141:       MatMatMultNumeric_SeqAIJ_SeqDense(A,B,AB_den);
    142:       MatMatMultNumeric_SeqAIJ_SeqDense(R,AB_den,RAB);
    143:       MatDestroy(&AB_den);
    144:       return(0);
    145:     }
    146:   }

    148:   MatDenseGetArrayRead(B,&b);
    149:   MatDenseGetArray(RAB,&d);
    150:   b1   = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;
    151:   c    = work; c2 = c + am; c3 = c2 + am; c4 = c3 + am;
    152:   for (col=0; col<cn-4; col += 4) {  /* over columns of C */
    153:     for (i=0; i<am; i++) {        /* over rows of A in those columns */
    154:       r1 = r2 = r3 = r4 = 0.0;
    155:       n  = ai[i+1] - ai[i];
    156:       aj = a->j + ai[i];
    157:       aa = a->a + ai[i];
    158:       for (j=0; j<n; j++) {
    159:         r1 += (*aa)*b1[*aj];
    160:         r2 += (*aa)*b2[*aj];
    161:         r3 += (*aa)*b3[*aj];
    162:         r4 += (*aa++)*b4[*aj++];
    163:       }
    164:       c[i]       = r1;
    165:       c[am  + i] = r2;
    166:       c[am2 + i] = r3;
    167:       c[am3 + i] = r4;
    168:     }
    169:     b1 += bm4;
    170:     b2 += bm4;
    171:     b3 += bm4;
    172:     b4 += bm4;

    174:     /* RAB[:,col] = R*C[:,col] */
    175:     colrm = col*rm;
    176:     for (i=0; i<rm; i++) {        /* over rows of R in those columns */
    177:       r1 = r2 = r3 = r4 = 0.0;
    178:       n  = r->i[i+1] - r->i[i];
    179:       rj = r->j + r->i[i];
    180:       ra = r->a + r->i[i];
    181:       for (j=0; j<n; j++) {
    182:         r1 += (*ra)*c[*rj];
    183:         r2 += (*ra)*c2[*rj];
    184:         r3 += (*ra)*c3[*rj];
    185:         r4 += (*ra++)*c4[*rj++];
    186:       }
    187:       d[colrm + i]       = r1;
    188:       d[colrm + rm + i]  = r2;
    189:       d[colrm + rm2 + i] = r3;
    190:       d[colrm + rm3 + i] = r4;
    191:     }
    192:   }
    193:   for (; col<cn; col++) {     /* over extra columns of C */
    194:     for (i=0; i<am; i++) {  /* over rows of A in those columns */
    195:       r1 = 0.0;
    196:       n  = a->i[i+1] - a->i[i];
    197:       aj = a->j + a->i[i];
    198:       aa = a->a + a->i[i];
    199:       for (j=0; j<n; j++) {
    200:         r1 += (*aa++)*b1[*aj++];
    201:       }
    202:       c[i] = r1;
    203:     }
    204:     b1 += bm;

    206:     for (i=0; i<rm; i++) {  /* over rows of R in those columns */
    207:       r1 = 0.0;
    208:       n  = r->i[i+1] - r->i[i];
    209:       rj = r->j + r->i[i];
    210:       ra = r->a + r->i[i];
    211:       for (j=0; j<n; j++) {
    212:         r1 += (*ra++)*c[*rj++];
    213:       }
    214:       d[col*rm + i] = r1;
    215:     }
    216:   }
    217:   PetscLogFlops(cn*2.0*(a->nz + r->nz));

    219:   MatDenseRestoreArrayRead(B,&b);
    220:   MatDenseRestoreArray(RAB,&d);
    221:   MatAssemblyBegin(RAB,MAT_FINAL_ASSEMBLY);
    222:   MatAssemblyEnd(RAB,MAT_FINAL_ASSEMBLY);
    223:   return(0);
    224: }

    226: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_colorrart(Mat A,Mat R,Mat C)
    227: {
    228:   PetscErrorCode       ierr;
    229:   Mat_RARt             *rart;
    230:   MatTransposeColoring matcoloring;
    231:   Mat                  Rt,RARt;

    234:   MatCheckProduct(C,3);
    235:   if (!C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data empty");
    236:   rart = (Mat_RARt*)C->product->data;

    238:   /* Get dense Rt by Apply MatTransposeColoring to R */
    239:   matcoloring = rart->matcoloring;
    240:   Rt          = rart->Rt;
    241:   MatTransColoringApplySpToDen(matcoloring,R,Rt);

    243:   /* Get dense RARt = R*A*Rt -- dominates! */
    244:   RARt = rart->RARt;
    245:   MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqDense(R,A,Rt,RARt,rart->work);

    247:   /* Recover C from C_dense */
    248:   MatTransColoringApplyDenToSp(matcoloring,RARt,C);
    249:   return(0);
    250: }

    252: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,PetscReal fill,Mat C)
    253: {
    255:   Mat            ARt;
    256:   Mat_RARt       *rart;
    257:   char           *alg;

    260:   MatCheckProduct(C,4);
    261:   if (C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data not empty");
    262:   /* create symbolic ARt = A*R^T  */
    263:   MatProductCreate(A,R,NULL,&ARt);
    264:   MatProductSetType(ARt,MATPRODUCT_ABt);
    265:   MatProductSetAlgorithm(ARt,"sorted");
    266:   MatProductSetFill(ARt,fill);
    267:   MatProductSetFromOptions(ARt);
    268:   MatProductSymbolic(ARt);

    270:   /* compute symbolic C = R*ARt */
    271:   /* set algorithm for C = R*ARt */
    272:   PetscStrallocpy(C->product->alg,&alg);
    273:   MatProductSetAlgorithm(C,"sorted");
    274:   MatMatMultSymbolic_SeqAIJ_SeqAIJ(R,ARt,fill,C);
    275:   /* resume original algorithm for C */
    276:   MatProductSetAlgorithm(C,alg);
    277:   PetscFree(alg);

    279:   C->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult;

    281:   PetscNew(&rart);
    282:   rart->ARt = ARt;
    283:   C->product->data    = rart;
    284:   C->product->destroy = MatDestroy_SeqAIJ_RARt;
    285:   PetscInfo(C,"Use ARt=A*R^T, C=R*ARt via MatMatTransposeMult(). Coloring can be applied to A*R^T.\n");
    286:   return(0);
    287: }

    289: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ_matmattransposemult(Mat A,Mat R,Mat C)
    290: {
    292:   Mat_RARt       *rart;

    295:   MatCheckProduct(C,3);
    296:   if (!C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data empty");
    297:   rart = (Mat_RARt*)C->product->data;
    298:   MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,R,rart->ARt); /* dominate! */
    299:   MatMatMultNumeric_SeqAIJ_SeqAIJ(R,rart->ARt,C);
    300:   return(0);
    301: }

    303: PetscErrorCode MatRARtSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat R,PetscReal fill,Mat C)
    304: {
    306:   Mat            Rt;
    307:   Mat_RARt       *rart;

    310:   MatCheckProduct(C,4);
    311:   if (C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data not empty");
    312:   MatTranspose_SeqAIJ(R,MAT_INITIAL_MATRIX,&Rt);
    313:   MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ(R,A,Rt,fill,C);

    315:   PetscNew(&rart);
    316:   rart->data = C->product->data;
    317:   rart->destroy = C->product->destroy;
    318:   rart->Rt = Rt;
    319:   C->product->data    = rart;
    320:   C->product->destroy = MatDestroy_SeqAIJ_RARt;
    321:   C->ops->rartnumeric = MatRARtNumeric_SeqAIJ_SeqAIJ;
    322:   PetscInfo(C,"Use Rt=R^T and C=R*A*Rt via MatMatMatMult() to avoid sparse inner products\n");
    323:   return(0);
    324: }

    326: PetscErrorCode MatRARtNumeric_SeqAIJ_SeqAIJ(Mat A,Mat R,Mat C)
    327: {
    329:   Mat_RARt       *rart;

    332:   MatCheckProduct(C,3);
    333:   if (!C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data empty");
    334:   rart = (Mat_RARt*)C->product->data;
    335:   MatTranspose_SeqAIJ(R,MAT_REUSE_MATRIX,&rart->Rt);
    336:   /* MatMatMatMultSymbolic used a different data */
    337:   C->product->data = rart->data;
    338:   MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ(R,A,rart->Rt,C);
    339:   C->product->data = rart;
    340:   return(0);
    341: }

    343: PetscErrorCode MatRARt_SeqAIJ_SeqAIJ(Mat A,Mat R,MatReuse scall,PetscReal fill,Mat *C)
    344: {
    346:   const char     *algTypes[3] = {"matmatmatmult","matmattransposemult","coloring_rart"};
    347:   PetscInt       alg=0; /* set default algorithm */

    350:   if (scall == MAT_INITIAL_MATRIX) {
    351:     PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatRARt","Mat");
    352:     PetscOptionsEList("-matrart_via","Algorithmic approach","MatRARt",algTypes,3,algTypes[0],&alg,NULL);
    353:     PetscOptionsEnd();

    355:     PetscLogEventBegin(MAT_RARtSymbolic,A,R,0,0);
    356:     MatCreate(PETSC_COMM_SELF,C);
    357:     switch (alg) {
    358:     case 1:
    359:       /* via matmattransposemult: ARt=A*R^T, C=R*ARt - matrix coloring can be applied to A*R^T */
    360:       MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(A,R,fill,*C);
    361:       break;
    362:     case 2:
    363:       /* via coloring_rart: apply coloring C = R*A*R^T                          */
    364:       MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(A,R,fill,*C);
    365:       break;
    366:     default:
    367:       /* via matmatmatmult: Rt=R^T, C=R*A*Rt - avoid inefficient sparse inner products */
    368:       MatRARtSymbolic_SeqAIJ_SeqAIJ(A,R,fill,*C);
    369:       break;
    370:     }
    371:     PetscLogEventEnd(MAT_RARtSymbolic,A,R,0,0);
    372:   }

    374:   PetscLogEventBegin(MAT_RARtNumeric,A,R,0,0);
    375:   ((*C)->ops->rartnumeric)(A,R,*C);
    376:   PetscLogEventEnd(MAT_RARtNumeric,A,R,0,0);
    377:   return(0);
    378: }

    380: /* ------------------------------------------------------------- */
    381: PetscErrorCode MatProductSymbolic_RARt_SeqAIJ_SeqAIJ(Mat C)
    382: {
    383:   PetscErrorCode      ierr;
    384:   Mat_Product         *product = C->product;
    385:   Mat                 A=product->A,R=product->B;
    386:   MatProductAlgorithm alg=product->alg;
    387:   PetscReal           fill=product->fill;
    388:   PetscBool           flg;

    391:   PetscStrcmp(alg,"r*a*rt",&flg);
    392:   if (flg) {
    393:     MatRARtSymbolic_SeqAIJ_SeqAIJ(A,R,fill,C);
    394:     goto next;
    395:   }

    397:   PetscStrcmp(alg,"r*art",&flg);
    398:   if (flg) {
    399:     MatRARtSymbolic_SeqAIJ_SeqAIJ_matmattransposemult(A,R,fill,C);
    400:     goto next;
    401:   }

    403:   PetscStrcmp(alg,"coloring_rart",&flg);
    404:   if (flg) {
    405:     MatRARtSymbolic_SeqAIJ_SeqAIJ_colorrart(A,R,fill,C);
    406:     goto next;
    407:   }

    409:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatProductAlgorithm is not supported");

    411: next:
    412:   C->ops->productnumeric = MatProductNumeric_RARt;
    413:   return(0);
    414: }