Actual source code: kaij.c

      1: /*
      2:   Defines the basic matrix operations for the KAIJ  matrix storage format.
      3:   This format is used to evaluate matrices of the form:

      5:     [I \otimes S + A \otimes T]

      7:   where
      8:     S is a dense (p \times q) matrix
      9:     T is a dense (p \times q) matrix
     10:     A is an AIJ  (n \times n) matrix
     11:     I is the identity matrix

     13:   The resulting matrix is (np \times nq)

     15:   We provide:
     16:      MatMult()
     17:      MatMultAdd()
     18:      MatInvertBlockDiagonal()
     19:   and
     20:      MatCreateKAIJ(Mat,PetscInt,PetscInt,const PetscScalar[],const PetscScalar[],Mat*)

     22:   This single directory handles both the sequential and parallel codes
     23: */

     25: #include <../src/mat/impls/kaij/kaij.h>
     26: #include <../src/mat/utils/freespace.h>
     27: #include <petsc/private/vecimpl.h>

     29: /*@
     30:   MatKAIJGetAIJ - Get the `MATAIJ` matrix describing the blockwise action of the `MATKAIJ` matrix

     32:   Not Collective, but if the `MATKAIJ` matrix is parallel, the `MATAIJ` matrix is also parallel

     34:   Input Parameter:
     35: . A - the `MATKAIJ` matrix

     37:   Output Parameter:
     38: . B - the `MATAIJ` matrix

     40:   Level: advanced

     42:   Note:
     43:   The reference count on the `MATAIJ` matrix is not increased so you should not destroy it.

     45: .seealso: [](ch_matrices), `Mat`, `MatCreateKAIJ()`, `MATKAIJ`, `MATAIJ`
     46: @*/
     47: PetscErrorCode MatKAIJGetAIJ(Mat A, Mat *B)
     48: {
     49:   PetscBool ismpikaij, isseqkaij;

     51:   PetscFunctionBegin;
     52:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIKAIJ, &ismpikaij));
     53:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQKAIJ, &isseqkaij));
     54:   if (ismpikaij) {
     55:     Mat_MPIKAIJ *b = (Mat_MPIKAIJ *)A->data;

     57:     *B = b->A;
     58:   } else if (isseqkaij) {
     59:     Mat_SeqKAIJ *b = (Mat_SeqKAIJ *)A->data;

     61:     *B = b->AIJ;
     62:   } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix passed in is not of type KAIJ");
     63:   PetscFunctionReturn(PETSC_SUCCESS);
     64: }

     66: /*@C
     67:   MatKAIJGetS - Get the `S` matrix describing the shift action of the `MATKAIJ` matrix

     69:   Not Collective; the entire `S` is stored and returned independently on all processes.

     71:   Input Parameter:
     72: . A - the `MATKAIJ` matrix

     74:   Output Parameters:
     75: + m - the number of rows in `S`
     76: . n - the number of columns in `S`
     77: - S - the S matrix, in form of a scalar array in column-major format

     79:   Level: advanced

     81:   Note:
     82:   All output parameters are optional (pass `NULL` if not desired)

     84: .seealso: [](ch_matrices), `Mat`, `MATKAIJ`, `MatCreateKAIJ()`, `MatGetBlockSizes()`
     85: @*/
     86: PetscErrorCode MatKAIJGetS(Mat A, PetscInt *m, PetscInt *n, PetscScalar *S[])
     87: {
     88:   Mat_SeqKAIJ *b = (Mat_SeqKAIJ *)A->data;

     90:   PetscFunctionBegin;
     91:   if (m) *m = b->p;
     92:   if (n) *n = b->q;
     93:   if (S) *S = b->S;
     94:   PetscFunctionReturn(PETSC_SUCCESS);
     95: }

     97: /*@C
     98:   MatKAIJGetSRead - Get a read-only pointer to the `S` matrix describing the shift action of the `MATKAIJ` matrix

    100:   Not Collective; the entire `S` is stored and returned independently on all processes.

    102:   Input Parameter:
    103: . A - the `MATKAIJ` matrix

    105:   Output Parameters:
    106: + m - the number of rows in `S`
    107: . n - the number of columns in `S`
    108: - S - the S matrix, in form of a scalar array in column-major format

    110:   Level: advanced

    112:   Note:
    113:   All output parameters are optional (pass `NULL` if not desired)

    115: .seealso: [](ch_matrices), `Mat`, `MATKAIJ`, `MatCreateKAIJ()`, `MatGetBlockSizes()`
    116: @*/
    117: PetscErrorCode MatKAIJGetSRead(Mat A, PetscInt *m, PetscInt *n, const PetscScalar *S[])
    118: {
    119:   Mat_SeqKAIJ *b = (Mat_SeqKAIJ *)A->data;

    121:   PetscFunctionBegin;
    122:   if (m) *m = b->p;
    123:   if (n) *n = b->q;
    124:   if (S) *S = b->S;
    125:   PetscFunctionReturn(PETSC_SUCCESS);
    126: }

    128: /*@C
    129:   MatKAIJRestoreS - Restore array obtained with `MatKAIJGetS()`

    131:   Not Collective

    133:   Input Parameters:
    134: + A - the `MATKAIJ` matrix
    135: - S - location of pointer to array obtained with `MatKAIJGetS()`

    137:   Level: advanced

    139:   Note:
    140:   This routine zeros the array pointer to prevent accidental reuse after it has been restored.
    141:   If `NULL` is passed, it will not attempt to zero the array pointer.

    143: .seealso: [](ch_matrices), `Mat`, `MATKAIJ`, `MatKAIJGetS()`, `MatKAIJGetSRead()`, `MatKAIJRestoreSRead()`
    144: @*/
    145: PetscErrorCode MatKAIJRestoreS(Mat A, PetscScalar *S[])
    146: {
    147:   PetscFunctionBegin;
    148:   if (S) *S = NULL;
    149:   PetscCall(PetscObjectStateIncrease((PetscObject)A));
    150:   PetscFunctionReturn(PETSC_SUCCESS);
    151: }

    153: /*@C
    154:   MatKAIJRestoreSRead - Restore array obtained with `MatKAIJGetSRead()`

    156:   Not Collective

    158:   Input Parameters:
    159: + A - the `MATKAIJ` matrix
    160: - S - location of pointer to array obtained with `MatKAIJGetS()`

    162:   Level: advanced

    164:   Note:
    165:   This routine zeros the array pointer to prevent accidental reuse after it has been restored.
    166:   If `NULL` is passed, it will not attempt to zero the array pointer.

    168: .seealso: [](ch_matrices), `Mat`, `MATKAIJ`, `MatKAIJGetS()`, `MatKAIJGetSRead()`
    169: @*/
    170: PetscErrorCode MatKAIJRestoreSRead(Mat A, const PetscScalar *S[])
    171: {
    172:   PetscFunctionBegin;
    173:   if (S) *S = NULL;
    174:   PetscFunctionReturn(PETSC_SUCCESS);
    175: }

    177: /*@C
    178:   MatKAIJGetT - Get the transformation matrix `T` associated with the `MATKAIJ` matrix

    180:   Not Collective; the entire `T` is stored and returned independently on all processes

    182:   Input Parameter:
    183: . A - the `MATKAIJ` matrix

    185:   Output Parameters:
    186: + m - the number of rows in `T`
    187: . n - the number of columns in `T`
    188: - T - the T matrix, in form of a scalar array in column-major format

    190:   Level: advanced

    192:   Note:
    193:   All output parameters are optional (pass `NULL` if not desired)

    195: .seealso: [](ch_matrices), `Mat`, `MATKAIJ`, `MatCreateKAIJ()`, `MatGetBlockSizes()`
    196: @*/
    197: PetscErrorCode MatKAIJGetT(Mat A, PetscInt *m, PetscInt *n, PetscScalar *T[])
    198: {
    199:   Mat_SeqKAIJ *b = (Mat_SeqKAIJ *)A->data;

    201:   PetscFunctionBegin;
    202:   if (m) *m = b->p;
    203:   if (n) *n = b->q;
    204:   if (T) *T = b->T;
    205:   PetscFunctionReturn(PETSC_SUCCESS);
    206: }

    208: /*@C
    209:   MatKAIJGetTRead - Get a read-only pointer to the transformation matrix `T` associated with the `MATKAIJ` matrix

    211:   Not Collective; the entire `T` is stored and returned independently on all processes

    213:   Input Parameter:
    214: . A - the `MATKAIJ` matrix

    216:   Output Parameters:
    217: + m - the number of rows in `T`
    218: . n - the number of columns in `T`
    219: - T - the T matrix, in form of a scalar array in column-major format

    221:   Level: advanced

    223:   Note:
    224:   All output parameters are optional (pass `NULL` if not desired)

    226: .seealso: [](ch_matrices), `Mat`, `MATKAIJ`, `MatCreateKAIJ()`, `MatGetBlockSizes()`
    227: @*/
    228: PetscErrorCode MatKAIJGetTRead(Mat A, PetscInt *m, PetscInt *n, const PetscScalar *T[])
    229: {
    230:   Mat_SeqKAIJ *b = (Mat_SeqKAIJ *)A->data;

    232:   PetscFunctionBegin;
    233:   if (m) *m = b->p;
    234:   if (n) *n = b->q;
    235:   if (T) *T = b->T;
    236:   PetscFunctionReturn(PETSC_SUCCESS);
    237: }

    239: /*@C
    240:   MatKAIJRestoreT - Restore array obtained with `MatKAIJGetT()`

    242:   Not Collective

    244:   Input Parameters:
    245: + A - the `MATKAIJ` matrix
    246: - T - location of pointer to array obtained with `MatKAIJGetS()`

    248:   Level: advanced

    250:   Note:
    251:   This routine zeros the array pointer to prevent accidental reuse after it has been restored.
    252:   If `NULL` is passed, it will not attempt to zero the array pointer.

    254: .seealso: [](ch_matrices), `Mat`, `MATKAIJ`, `MatKAIJGetT()`, `MatKAIJGetTRead()`, `MatKAIJRestoreTRead()`
    255: @*/
    256: PetscErrorCode MatKAIJRestoreT(Mat A, PetscScalar *T[])
    257: {
    258:   PetscFunctionBegin;
    259:   if (T) *T = NULL;
    260:   PetscCall(PetscObjectStateIncrease((PetscObject)A));
    261:   PetscFunctionReturn(PETSC_SUCCESS);
    262: }

    264: /*@C
    265:   MatKAIJRestoreTRead - Restore array obtained with `MatKAIJGetTRead()`

    267:   Not Collective

    269:   Input Parameters:
    270: + A - the `MATKAIJ` matrix
    271: - T - location of pointer to array obtained with `MatKAIJGetS()`

    273:   Level: advanced

    275:   Note:
    276:   This routine zeros the array pointer to prevent accidental reuse after it has been restored.
    277:   If `NULL` is passed, it will not attempt to zero the array pointer.

    279: .seealso: [](ch_matrices), `Mat`, `MATKAIJ`, `MatKAIJGetT()`, `MatKAIJGetTRead()`
    280: @*/
    281: PetscErrorCode MatKAIJRestoreTRead(Mat A, const PetscScalar *T[])
    282: {
    283:   PetscFunctionBegin;
    284:   if (T) *T = NULL;
    285:   PetscFunctionReturn(PETSC_SUCCESS);
    286: }

    288: /*@
    289:   MatKAIJSetAIJ - Set the `MATAIJ` matrix describing the blockwise action of the `MATKAIJ` matrix

    291:   Logically Collective; if the `MATAIJ` matrix is parallel, the `MATKAIJ` matrix is also parallel

    293:   Input Parameters:
    294: + A - the `MATKAIJ` matrix
    295: - B - the `MATAIJ` matrix

    297:   Level: advanced

    299:   Notes:
    300:   This function increases the reference count on the `MATAIJ` matrix, so the user is free to destroy the matrix if it is not needed.

    302:   Changes to the entries of the `MATAIJ` matrix will immediately affect the `MATKAIJ` matrix.

    304: .seealso: [](ch_matrices), `Mat`, `MATKAIJ`, `MatKAIJGetAIJ()`, `MatKAIJSetS()`, `MatKAIJSetT()`
    305: @*/
    306: PetscErrorCode MatKAIJSetAIJ(Mat A, Mat B)
    307: {
    308:   PetscMPIInt size;
    309:   PetscBool   flg;

    311:   PetscFunctionBegin;
    312:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
    313:   if (size == 1) {
    314:     PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSEQAIJ, &flg));
    315:     PetscCheck(flg, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "MatKAIJSetAIJ() with MATSEQKAIJ does not support %s as the AIJ mat", ((PetscObject)B)->type_name);
    316:     Mat_SeqKAIJ *a = (Mat_SeqKAIJ *)A->data;
    317:     a->AIJ         = B;
    318:   } else {
    319:     Mat_MPIKAIJ *a = (Mat_MPIKAIJ *)A->data;
    320:     a->A           = B;
    321:   }
    322:   PetscCall(PetscObjectReference((PetscObject)B));
    323:   PetscFunctionReturn(PETSC_SUCCESS);
    324: }

    326: /*@
    327:   MatKAIJSetS - Set the `S` matrix describing the shift action of the `MATKAIJ` matrix

    329:   Logically Collective; the entire `S` is stored independently on all processes.

    331:   Input Parameters:
    332: + A - the `MATKAIJ` matrix
    333: . p - the number of rows in `S`
    334: . q - the number of columns in `S`
    335: - S - the S matrix, in form of a scalar array in column-major format

    337:   Level: advanced

    339:   Notes:
    340:   The dimensions `p` and `q` must match those of the transformation matrix `T` associated with the `MATKAIJ` matrix.

    342:   The `S` matrix is copied, so the user can destroy this array.

    344: .seealso: [](ch_matrices), `Mat`, `MATKAIJ`, `MatKAIJGetS()`, `MatKAIJSetT()`, `MatKAIJSetAIJ()`
    345: @*/
    346: PetscErrorCode MatKAIJSetS(Mat A, PetscInt p, PetscInt q, const PetscScalar S[])
    347: {
    348:   Mat_SeqKAIJ *a = (Mat_SeqKAIJ *)A->data;

    350:   PetscFunctionBegin;
    351:   PetscCall(PetscFree(a->S));
    352:   if (S) {
    353:     PetscCall(PetscMalloc1(p * q, &a->S));
    354:     PetscCall(PetscArraycpy(a->S, S, p * q));
    355:   } else a->S = NULL;

    357:   a->p = p;
    358:   a->q = q;
    359:   PetscFunctionReturn(PETSC_SUCCESS);
    360: }

    362: /*@
    363:   MatKAIJGetScaledIdentity - Check if both `S` and `T` are scaled identities.

    365:   Logically Collective.

    367:   Input Parameter:
    368: . A - the `MATKAIJ` matrix

    370:   Output Parameter:
    371: . identity - the Boolean value

    373:   Level: advanced

    375: .seealso: [](ch_matrices), `Mat`, `MATKAIJ`, `MatKAIJGetS()`, `MatKAIJGetT()`
    376: @*/
    377: PetscErrorCode MatKAIJGetScaledIdentity(Mat A, PetscBool *identity)
    378: {
    379:   Mat_SeqKAIJ *a = (Mat_SeqKAIJ *)A->data;
    380:   PetscInt     i, j;

    382:   PetscFunctionBegin;
    383:   if (a->p != a->q) {
    384:     *identity = PETSC_FALSE;
    385:     PetscFunctionReturn(PETSC_SUCCESS);
    386:   } else *identity = PETSC_TRUE;
    387:   if (!a->isTI || a->S) {
    388:     for (i = 0; i < a->p && *identity; i++) {
    389:       for (j = 0; j < a->p && *identity; j++) {
    390:         if (i != j) {
    391:           if (a->S && PetscAbsScalar(a->S[i + j * a->p]) > PETSC_SMALL) *identity = PETSC_FALSE;
    392:           if (a->T && PetscAbsScalar(a->T[i + j * a->p]) > PETSC_SMALL) *identity = PETSC_FALSE;
    393:         } else {
    394:           if (a->S && PetscAbsScalar(a->S[i * (a->p + 1)] - a->S[0]) > PETSC_SMALL) *identity = PETSC_FALSE;
    395:           if (a->T && PetscAbsScalar(a->T[i * (a->p + 1)] - a->T[0]) > PETSC_SMALL) *identity = PETSC_FALSE;
    396:         }
    397:       }
    398:     }
    399:   }
    400:   PetscFunctionReturn(PETSC_SUCCESS);
    401: }

    403: /*@
    404:   MatKAIJSetT - Set the transformation matrix `T` associated with the `MATKAIJ` matrix

    406:   Logically Collective; the entire `T` is stored independently on all processes.

    408:   Input Parameters:
    409: + A - the `MATKAIJ` matrix
    410: . p - the number of rows in `S`
    411: . q - the number of columns in `S`
    412: - T - the `T` matrix, in form of a scalar array in column-major format

    414:   Level: advanced

    416:   Notes:
    417:   The dimensions `p` and `q` must match those of the shift matrix `S` associated with the `MATKAIJ` matrix.

    419:   The `T` matrix is copied, so the user can destroy this array.

    421: .seealso: [](ch_matrices), `Mat`, `MATKAIJ`, `MatKAIJGetT()`, `MatKAIJSetS()`, `MatKAIJSetAIJ()`
    422: @*/
    423: PetscErrorCode MatKAIJSetT(Mat A, PetscInt p, PetscInt q, const PetscScalar T[])
    424: {
    425:   PetscInt     i, j;
    426:   Mat_SeqKAIJ *a    = (Mat_SeqKAIJ *)A->data;
    427:   PetscBool    isTI = PETSC_FALSE;

    429:   PetscFunctionBegin;
    430:   /* check if T is an identity matrix */
    431:   if (T && (p == q)) {
    432:     isTI = PETSC_TRUE;
    433:     for (i = 0; i < p; i++) {
    434:       for (j = 0; j < q; j++) {
    435:         if (i == j) {
    436:           /* diagonal term must be 1 */
    437:           if (T[i + j * p] != 1.0) isTI = PETSC_FALSE;
    438:         } else {
    439:           /* off-diagonal term must be 0 */
    440:           if (T[i + j * p] != 0.0) isTI = PETSC_FALSE;
    441:         }
    442:       }
    443:     }
    444:   }
    445:   a->isTI = isTI;

    447:   PetscCall(PetscFree(a->T));
    448:   if (T && (!isTI)) {
    449:     PetscCall(PetscMalloc1(p * q, &a->T));
    450:     PetscCall(PetscArraycpy(a->T, T, p * q));
    451:   } else a->T = NULL;

    453:   a->p = p;
    454:   a->q = q;
    455:   PetscFunctionReturn(PETSC_SUCCESS);
    456: }

    458: static PetscErrorCode MatDestroy_SeqKAIJ(Mat A)
    459: {
    460:   Mat_SeqKAIJ *b = (Mat_SeqKAIJ *)A->data;

    462:   PetscFunctionBegin;
    463:   PetscCall(MatDestroy(&b->AIJ));
    464:   PetscCall(PetscFree(b->S));
    465:   PetscCall(PetscFree(b->T));
    466:   PetscCall(PetscFree(b->ibdiag));
    467:   PetscCall(PetscFree5(b->sor.w, b->sor.y, b->sor.work, b->sor.t, b->sor.arr));
    468:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqkaij_seqaij_C", NULL));
    469:   PetscCall(PetscFree(A->data));
    470:   PetscFunctionReturn(PETSC_SUCCESS);
    471: }

    473: static PetscErrorCode MatKAIJ_build_AIJ_OAIJ(Mat A)
    474: {
    475:   Mat_MPIKAIJ     *a;
    476:   Mat_MPIAIJ      *mpiaij;
    477:   PetscScalar     *T;
    478:   PetscInt         i, j;
    479:   PetscObjectState state;

    481:   PetscFunctionBegin;
    482:   a      = (Mat_MPIKAIJ *)A->data;
    483:   mpiaij = (Mat_MPIAIJ *)a->A->data;

    485:   PetscCall(PetscObjectStateGet((PetscObject)a->A, &state));
    486:   if (state == a->state) {
    487:     /* The existing AIJ and KAIJ members are up-to-date, so simply exit. */
    488:     PetscFunctionReturn(PETSC_SUCCESS);
    489:   } else {
    490:     PetscCall(MatDestroy(&a->AIJ));
    491:     PetscCall(MatDestroy(&a->OAIJ));
    492:     if (a->isTI) {
    493:       /* If the transformation matrix associated with the parallel matrix A is the identity matrix, then a->T will be NULL.
    494:        * In this case, if we pass a->T directly to the MatCreateKAIJ() calls to create the sequential submatrices, the routine will
    495:        * not be able to tell that transformation matrix should be set to the identity; thus we create a temporary identity matrix
    496:        * to pass in. */
    497:       PetscCall(PetscMalloc1(a->p * a->q, &T));
    498:       for (i = 0; i < a->p; i++) {
    499:         for (j = 0; j < a->q; j++) {
    500:           if (i == j) T[i + j * a->p] = 1.0;
    501:           else T[i + j * a->p] = 0.0;
    502:         }
    503:       }
    504:     } else T = a->T;
    505:     PetscCall(MatCreateKAIJ(mpiaij->A, a->p, a->q, a->S, T, &a->AIJ));
    506:     PetscCall(MatCreateKAIJ(mpiaij->B, a->p, a->q, NULL, T, &a->OAIJ));
    507:     if (a->isTI) PetscCall(PetscFree(T));
    508:     a->state = state;
    509:   }
    510:   PetscFunctionReturn(PETSC_SUCCESS);
    511: }

    513: static PetscErrorCode MatSetUp_KAIJ(Mat A)
    514: {
    515:   PetscInt     n;
    516:   PetscMPIInt  size;
    517:   Mat_SeqKAIJ *seqkaij = (Mat_SeqKAIJ *)A->data;

    519:   PetscFunctionBegin;
    520:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
    521:   if (size == 1) {
    522:     PetscCall(MatSetSizes(A, seqkaij->p * seqkaij->AIJ->rmap->n, seqkaij->q * seqkaij->AIJ->cmap->n, seqkaij->p * seqkaij->AIJ->rmap->N, seqkaij->q * seqkaij->AIJ->cmap->N));
    523:     PetscCall(PetscLayoutSetBlockSize(A->rmap, seqkaij->p));
    524:     PetscCall(PetscLayoutSetBlockSize(A->cmap, seqkaij->q));
    525:     PetscCall(PetscLayoutSetUp(A->rmap));
    526:     PetscCall(PetscLayoutSetUp(A->cmap));
    527:   } else {
    528:     Mat_MPIKAIJ *a;
    529:     Mat_MPIAIJ  *mpiaij;
    530:     IS           from, to;
    531:     Vec          gvec;

    533:     a      = (Mat_MPIKAIJ *)A->data;
    534:     mpiaij = (Mat_MPIAIJ *)a->A->data;
    535:     PetscCall(MatSetSizes(A, a->p * a->A->rmap->n, a->q * a->A->cmap->n, a->p * a->A->rmap->N, a->q * a->A->cmap->N));
    536:     PetscCall(PetscLayoutSetBlockSize(A->rmap, seqkaij->p));
    537:     PetscCall(PetscLayoutSetBlockSize(A->cmap, seqkaij->q));
    538:     PetscCall(PetscLayoutSetUp(A->rmap));
    539:     PetscCall(PetscLayoutSetUp(A->cmap));

    541:     PetscCall(MatKAIJ_build_AIJ_OAIJ(A));

    543:     PetscCall(VecGetSize(mpiaij->lvec, &n));
    544:     PetscCall(VecCreate(PETSC_COMM_SELF, &a->w));
    545:     PetscCall(VecSetSizes(a->w, n * a->q, n * a->q));
    546:     PetscCall(VecSetBlockSize(a->w, a->q));
    547:     PetscCall(VecSetType(a->w, VECSEQ));

    549:     /* create two temporary Index sets for build scatter gather */
    550:     PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)a->A), a->q, n, mpiaij->garray, PETSC_COPY_VALUES, &from));
    551:     PetscCall(ISCreateStride(PETSC_COMM_SELF, n * a->q, 0, 1, &to));

    553:     /* create temporary global vector to generate scatter context */
    554:     PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)a->A), a->q, a->q * a->A->cmap->n, a->q * a->A->cmap->N, NULL, &gvec));

    556:     /* generate the scatter context */
    557:     PetscCall(VecScatterCreate(gvec, from, a->w, to, &a->ctx));

    559:     PetscCall(ISDestroy(&from));
    560:     PetscCall(ISDestroy(&to));
    561:     PetscCall(VecDestroy(&gvec));
    562:   }

    564:   A->assembled = PETSC_TRUE;
    565:   PetscFunctionReturn(PETSC_SUCCESS);
    566: }

    568: static PetscErrorCode MatView_KAIJ(Mat A, PetscViewer viewer)
    569: {
    570:   PetscViewerFormat format;
    571:   Mat_SeqKAIJ      *a = (Mat_SeqKAIJ *)A->data;
    572:   Mat               B;
    573:   PetscInt          i;
    574:   PetscBool         ismpikaij;

    576:   PetscFunctionBegin;
    577:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIKAIJ, &ismpikaij));
    578:   PetscCall(PetscViewerGetFormat(viewer, &format));
    579:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL || format == PETSC_VIEWER_ASCII_IMPL) {
    580:     PetscCall(PetscViewerASCIIPrintf(viewer, "S and T have %" PetscInt_FMT " rows and %" PetscInt_FMT " columns\n", a->p, a->q));

    582:     /* Print appropriate details for S. */
    583:     if (!a->S) {
    584:       PetscCall(PetscViewerASCIIPrintf(viewer, "S is NULL\n"));
    585:     } else if (format == PETSC_VIEWER_ASCII_IMPL) {
    586:       PetscCall(PetscViewerASCIIPrintf(viewer, "Entries of S are "));
    587:       for (i = 0; i < (a->p * a->q); i++) {
    588: #if defined(PETSC_USE_COMPLEX)
    589:         PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e %18.16e ", (double)PetscRealPart(a->S[i]), (double)PetscImaginaryPart(a->S[i])));
    590: #else
    591:         PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e ", (double)PetscRealPart(a->S[i])));
    592: #endif
    593:       }
    594:       PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
    595:     }

    597:     /* Print appropriate details for T. */
    598:     if (a->isTI) {
    599:       PetscCall(PetscViewerASCIIPrintf(viewer, "T is the identity matrix\n"));
    600:     } else if (!a->T) {
    601:       PetscCall(PetscViewerASCIIPrintf(viewer, "T is NULL\n"));
    602:     } else if (format == PETSC_VIEWER_ASCII_IMPL) {
    603:       PetscCall(PetscViewerASCIIPrintf(viewer, "Entries of T are "));
    604:       for (i = 0; i < (a->p * a->q); i++) {
    605: #if defined(PETSC_USE_COMPLEX)
    606:         PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e %18.16e ", (double)PetscRealPart(a->T[i]), (double)PetscImaginaryPart(a->T[i])));
    607: #else
    608:         PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e ", (double)PetscRealPart(a->T[i])));
    609: #endif
    610:       }
    611:       PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
    612:     }

    614:     /* Now print details for the AIJ matrix, using the AIJ viewer. */
    615:     PetscCall(PetscViewerASCIIPrintf(viewer, "Now viewing the associated AIJ matrix:\n"));
    616:     if (ismpikaij) {
    617:       Mat_MPIKAIJ *b = (Mat_MPIKAIJ *)A->data;
    618:       PetscCall(MatView(b->A, viewer));
    619:     } else {
    620:       PetscCall(MatView(a->AIJ, viewer));
    621:     }

    623:   } else {
    624:     /* For all other matrix viewer output formats, simply convert to an AIJ matrix and call MatView() on that. */
    625:     PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B));
    626:     PetscCall(MatView(B, viewer));
    627:     PetscCall(MatDestroy(&B));
    628:   }
    629:   PetscFunctionReturn(PETSC_SUCCESS);
    630: }

    632: static PetscErrorCode MatDestroy_MPIKAIJ(Mat A)
    633: {
    634:   Mat_MPIKAIJ *b = (Mat_MPIKAIJ *)A->data;

    636:   PetscFunctionBegin;
    637:   PetscCall(MatDestroy(&b->AIJ));
    638:   PetscCall(MatDestroy(&b->OAIJ));
    639:   PetscCall(MatDestroy(&b->A));
    640:   PetscCall(VecScatterDestroy(&b->ctx));
    641:   PetscCall(VecDestroy(&b->w));
    642:   PetscCall(PetscFree(b->S));
    643:   PetscCall(PetscFree(b->T));
    644:   PetscCall(PetscFree(b->ibdiag));
    645:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatGetDiagonalBlock_C", NULL));
    646:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_mpikaij_mpiaij_C", NULL));
    647:   PetscCall(PetscFree(A->data));
    648:   PetscFunctionReturn(PETSC_SUCCESS);
    649: }

    651: /* zz = yy + Axx */
    652: static PetscErrorCode MatMultAdd_SeqKAIJ(Mat A, Vec xx, Vec yy, Vec zz)
    653: {
    654:   Mat_SeqKAIJ       *b = (Mat_SeqKAIJ *)A->data;
    655:   Mat_SeqAIJ        *a = (Mat_SeqAIJ *)b->AIJ->data;
    656:   const PetscScalar *s = b->S, *t = b->T;
    657:   const PetscScalar *x, *v, *bx;
    658:   PetscScalar       *y, *sums;
    659:   const PetscInt     m = b->AIJ->rmap->n, *idx, *ii;
    660:   PetscInt           n, i, jrow, j, l, p = b->p, q = b->q, k;

    662:   PetscFunctionBegin;
    663:   if (!yy) PetscCall(VecSet(zz, 0.0));
    664:   else PetscCall(VecCopy(yy, zz));
    665:   if ((!s) && (!t) && (!b->isTI)) PetscFunctionReturn(PETSC_SUCCESS);

    667:   PetscCall(VecGetArrayRead(xx, &x));
    668:   PetscCall(VecGetArray(zz, &y));
    669:   idx = a->j;
    670:   v   = a->a;
    671:   ii  = a->i;

    673:   if (b->isTI) {
    674:     for (i = 0; i < m; i++) {
    675:       jrow = ii[i];
    676:       n    = ii[i + 1] - jrow;
    677:       sums = y + p * i;
    678:       for (j = 0; j < n; j++) {
    679:         for (k = 0; k < p; k++) sums[k] += v[jrow + j] * x[q * idx[jrow + j] + k];
    680:       }
    681:     }
    682:     PetscCall(PetscLogFlops(3.0 * (a->nz) * p));
    683:   } else if (t) {
    684:     for (i = 0; i < m; i++) {
    685:       jrow = ii[i];
    686:       n    = ii[i + 1] - jrow;
    687:       sums = y + p * i;
    688:       for (j = 0; j < n; j++) {
    689:         for (k = 0; k < p; k++) {
    690:           for (l = 0; l < q; l++) sums[k] += v[jrow + j] * t[k + l * p] * x[q * idx[jrow + j] + l];
    691:         }
    692:       }
    693:     }
    694:     /* The flop count below assumes that v[jrow+j] is hoisted out (which an optimizing compiler is likely to do),
    695:      * and also that T part is hoisted outside this loop (in exchange for temporary storage) as (A \otimes I) (I \otimes T),
    696:      * so that this multiply doesn't have to be redone for each matrix entry, but just once per column. The latter
    697:      * transformation is much less likely to be applied, but we nonetheless count the minimum flops required. */
    698:     PetscCall(PetscLogFlops((2.0 * p * q - p) * m + 2.0 * p * a->nz));
    699:   }
    700:   if (s) {
    701:     for (i = 0; i < m; i++) {
    702:       sums = y + p * i;
    703:       bx   = x + q * i;
    704:       if (i < b->AIJ->cmap->n) {
    705:         for (j = 0; j < q; j++) {
    706:           for (k = 0; k < p; k++) sums[k] += s[k + j * p] * bx[j];
    707:         }
    708:       }
    709:     }
    710:     PetscCall(PetscLogFlops(2.0 * m * p * q));
    711:   }

    713:   PetscCall(VecRestoreArrayRead(xx, &x));
    714:   PetscCall(VecRestoreArray(zz, &y));
    715:   PetscFunctionReturn(PETSC_SUCCESS);
    716: }

    718: static PetscErrorCode MatMult_SeqKAIJ(Mat A, Vec xx, Vec yy)
    719: {
    720:   PetscFunctionBegin;
    721:   PetscCall(MatMultAdd_SeqKAIJ(A, xx, NULL, yy));
    722:   PetscFunctionReturn(PETSC_SUCCESS);
    723: }

    725: #include <petsc/private/kernels/blockinvert.h>

    727: static PetscErrorCode MatInvertBlockDiagonal_SeqKAIJ(Mat A, const PetscScalar **values)
    728: {
    729:   Mat_SeqKAIJ       *b = (Mat_SeqKAIJ *)A->data;
    730:   Mat_SeqAIJ        *a = (Mat_SeqAIJ *)b->AIJ->data;
    731:   const PetscScalar *S = b->S;
    732:   const PetscScalar *T = b->T;
    733:   const PetscScalar *v = a->a;
    734:   const PetscInt     p = b->p, q = b->q, m = b->AIJ->rmap->n, *idx = a->j, *ii = a->i;
    735:   PetscInt           i, j, *v_pivots, dof, dof2;
    736:   PetscScalar       *diag, aval, *v_work;

    738:   PetscFunctionBegin;
    739:   PetscCheck(p == q, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MATKAIJ: Block size must be square to calculate inverse.");
    740:   PetscCheck(S || T || b->isTI, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MATKAIJ: Cannot invert a zero matrix.");

    742:   dof  = p;
    743:   dof2 = dof * dof;

    745:   if (b->ibdiagvalid) {
    746:     if (values) *values = b->ibdiag;
    747:     PetscFunctionReturn(PETSC_SUCCESS);
    748:   }
    749:   if (!b->ibdiag) PetscCall(PetscMalloc1(dof2 * m, &b->ibdiag));
    750:   if (values) *values = b->ibdiag;
    751:   diag = b->ibdiag;

    753:   PetscCall(PetscMalloc2(dof, &v_work, dof, &v_pivots));
    754:   for (i = 0; i < m; i++) {
    755:     if (S) {
    756:       PetscCall(PetscArraycpy(diag, S, dof2));
    757:     } else {
    758:       PetscCall(PetscArrayzero(diag, dof2));
    759:     }
    760:     if (b->isTI) {
    761:       aval = 0;
    762:       for (j = ii[i]; j < ii[i + 1]; j++)
    763:         if (idx[j] == i) aval = v[j];
    764:       for (j = 0; j < dof; j++) diag[j + dof * j] += aval;
    765:     } else if (T) {
    766:       aval = 0;
    767:       for (j = ii[i]; j < ii[i + 1]; j++)
    768:         if (idx[j] == i) aval = v[j];
    769:       for (j = 0; j < dof2; j++) diag[j] += aval * T[j];
    770:     }
    771:     PetscCall(PetscKernel_A_gets_inverse_A(dof, diag, v_pivots, v_work, PETSC_FALSE, NULL));
    772:     diag += dof2;
    773:   }
    774:   PetscCall(PetscFree2(v_work, v_pivots));

    776:   b->ibdiagvalid = PETSC_TRUE;
    777:   PetscFunctionReturn(PETSC_SUCCESS);
    778: }

    780: static PetscErrorCode MatGetDiagonalBlock_MPIKAIJ(Mat A, Mat *B)
    781: {
    782:   Mat_MPIKAIJ *kaij = (Mat_MPIKAIJ *)A->data;

    784:   PetscFunctionBegin;
    785:   *B = kaij->AIJ;
    786:   PetscFunctionReturn(PETSC_SUCCESS);
    787: }

    789: static PetscErrorCode MatConvert_KAIJ_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
    790: {
    791:   Mat_SeqKAIJ    *a = (Mat_SeqKAIJ *)A->data;
    792:   Mat             AIJ, OAIJ, B;
    793:   PetscInt       *d_nnz, *o_nnz = NULL, nz, i, j, m, d;
    794:   const PetscInt  p = a->p, q = a->q;
    795:   PetscBool       ismpikaij, diagDense;
    796:   const PetscInt *aijdiag;

    798:   PetscFunctionBegin;
    799:   if (reuse != MAT_REUSE_MATRIX) {
    800:     PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIKAIJ, &ismpikaij));
    801:     if (ismpikaij) {
    802:       Mat_MPIKAIJ *b = (Mat_MPIKAIJ *)A->data;
    803:       AIJ            = ((Mat_SeqKAIJ *)b->AIJ->data)->AIJ;
    804:       OAIJ           = ((Mat_SeqKAIJ *)b->OAIJ->data)->AIJ;
    805:     } else {
    806:       AIJ  = a->AIJ;
    807:       OAIJ = NULL;
    808:     }
    809:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
    810:     PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
    811:     PetscCall(MatSetType(B, MATAIJ));
    812:     PetscCall(MatGetSize(AIJ, &m, NULL));

    814:     /*
    815:       Determine the first row of AIJ with missing diagonal and assume all successive rows also have a missing diagonal
    816:     */
    817:     PetscCall(MatGetDiagonalMarkers_SeqAIJ(AIJ, &aijdiag, &diagDense));
    818:     if (diagDense || !a->S) d = m;
    819:     else {
    820:       Mat_SeqAIJ *aij = (Mat_SeqAIJ *)AIJ->data;

    822:       for (d = 0; d < m; d++) {
    823:         if (aijdiag[d] == aij->i[d + 1]) break;
    824:       }
    825:     }
    826:     PetscCall(PetscMalloc1(m * p, &d_nnz));
    827:     for (i = 0; i < m; ++i) {
    828:       PetscCall(MatGetRow_SeqAIJ(AIJ, i, &nz, NULL, NULL));
    829:       for (j = 0; j < p; ++j) d_nnz[i * p + j] = nz * q + (i >= d) * q;
    830:       PetscCall(MatRestoreRow_SeqAIJ(AIJ, i, &nz, NULL, NULL));
    831:     }
    832:     if (OAIJ) {
    833:       PetscCall(PetscMalloc1(m * p, &o_nnz));
    834:       for (i = 0; i < m; ++i) {
    835:         PetscCall(MatGetRow_SeqAIJ(OAIJ, i, &nz, NULL, NULL));
    836:         for (j = 0; j < p; ++j) o_nnz[i * p + j] = nz * q;
    837:         PetscCall(MatRestoreRow_SeqAIJ(OAIJ, i, &nz, NULL, NULL));
    838:       }
    839:       PetscCall(MatMPIAIJSetPreallocation(B, 0, d_nnz, 0, o_nnz));
    840:     } else {
    841:       PetscCall(MatSeqAIJSetPreallocation(B, 0, d_nnz));
    842:     }
    843:     PetscCall(PetscFree(d_nnz));
    844:     PetscCall(PetscFree(o_nnz));
    845:   } else B = *newmat;
    846:   PetscCall(MatConvert_Basic(A, newtype, MAT_REUSE_MATRIX, &B));
    847:   if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
    848:   else *newmat = B;
    849:   PetscFunctionReturn(PETSC_SUCCESS);
    850: }

    852: static PetscErrorCode MatSOR_SeqKAIJ(Mat A, Vec bb, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec xx)
    853: {
    854:   Mat_SeqKAIJ       *kaij = (Mat_SeqKAIJ *)A->data;
    855:   Mat_SeqAIJ        *a    = (Mat_SeqAIJ *)kaij->AIJ->data;
    856:   const PetscScalar *aa = a->a, *T = kaij->T, *v;
    857:   const PetscInt     m = kaij->AIJ->rmap->n, *ai = a->i, *aj = a->j, p = kaij->p, q = kaij->q, *diag, *vi;
    858:   const PetscScalar *b, *xb, *idiag;
    859:   PetscScalar       *x, *work, *workt, *w, *y, *arr, *t, *arrt;
    860:   PetscInt           i, j, k, i2, bs, bs2, nz;

    862:   PetscFunctionBegin;
    863:   its = its * lits;
    864:   PetscCheck(!(flag & SOR_EISENSTAT), PETSC_COMM_SELF, PETSC_ERR_SUP, "No support yet for Eisenstat");
    865:   PetscCheck(its > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Relaxation requires global its %" PetscInt_FMT " and local its %" PetscInt_FMT " both positive", its, lits);
    866:   PetscCheck(!fshift, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for diagonal shift");
    867:   PetscCheck(!(flag & SOR_APPLY_UPPER) && !(flag & SOR_APPLY_LOWER), PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for applying upper or lower triangular parts");
    868:   PetscCheck(p == q, PETSC_COMM_SELF, PETSC_ERR_SUP, "MatSOR for KAIJ: No support for non-square dense blocks");
    869:   bs  = p;
    870:   bs2 = bs * bs;

    872:   if (!m) PetscFunctionReturn(PETSC_SUCCESS);

    874:   if (!kaij->ibdiagvalid) PetscCall(MatInvertBlockDiagonal_SeqKAIJ(A, NULL));
    875:   idiag = kaij->ibdiag;
    876:   PetscCall(MatGetDiagonalMarkers_SeqAIJ(kaij->AIJ, &diag, NULL));

    878:   if (!kaij->sor.setup) {
    879:     PetscCall(PetscMalloc5(bs, &kaij->sor.w, bs, &kaij->sor.y, m * bs, &kaij->sor.work, m * bs, &kaij->sor.t, m * bs2, &kaij->sor.arr));
    880:     kaij->sor.setup = PETSC_TRUE;
    881:   }
    882:   y    = kaij->sor.y;
    883:   w    = kaij->sor.w;
    884:   work = kaij->sor.work;
    885:   t    = kaij->sor.t;
    886:   arr  = kaij->sor.arr;

    888:   PetscCall(VecGetArray(xx, &x));
    889:   PetscCall(VecGetArrayRead(bb, &b));

    891:   if (flag & SOR_ZERO_INITIAL_GUESS) {
    892:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
    893:       PetscKernel_w_gets_Ar_times_v(bs, bs, b, idiag, x); /* x[0:bs] <- D^{-1} b[0:bs] */
    894:       PetscCall(PetscArraycpy(t, b, bs));
    895:       i2 = bs;
    896:       idiag += bs2;
    897:       for (i = 1; i < m; i++) {
    898:         v  = aa + ai[i];
    899:         vi = aj + ai[i];
    900:         nz = diag[i] - ai[i];

    902:         if (T) { /* b - T (Arow * x) */
    903:           PetscCall(PetscArrayzero(w, bs));
    904:           for (j = 0; j < nz; j++) {
    905:             for (k = 0; k < bs; k++) w[k] -= v[j] * x[vi[j] * bs + k];
    906:           }
    907:           PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs, w, T, &t[i2]);
    908:           for (k = 0; k < bs; k++) t[i2 + k] += b[i2 + k];
    909:         } else if (kaij->isTI) {
    910:           PetscCall(PetscArraycpy(t + i2, b + i2, bs));
    911:           for (j = 0; j < nz; j++) {
    912:             for (k = 0; k < bs; k++) t[i2 + k] -= v[j] * x[vi[j] * bs + k];
    913:           }
    914:         } else {
    915:           PetscCall(PetscArraycpy(t + i2, b + i2, bs));
    916:         }

    918:         PetscKernel_w_gets_Ar_times_v(bs, bs, t + i2, idiag, y);
    919:         for (j = 0; j < bs; j++) x[i2 + j] = omega * y[j];

    921:         idiag += bs2;
    922:         i2 += bs;
    923:       }
    924:       /* for logging purposes assume number of nonzero in lower half is 1/2 of total */
    925:       PetscCall(PetscLogFlops(1.0 * bs2 * a->nz));
    926:       xb = t;
    927:     } else xb = b;
    928:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
    929:       idiag = kaij->ibdiag + bs2 * (m - 1);
    930:       i2    = bs * (m - 1);
    931:       PetscCall(PetscArraycpy(w, xb + i2, bs));
    932:       PetscKernel_w_gets_Ar_times_v(bs, bs, w, idiag, x + i2);
    933:       i2 -= bs;
    934:       idiag -= bs2;
    935:       for (i = m - 2; i >= 0; i--) {
    936:         v  = aa + diag[i] + 1;
    937:         vi = aj + diag[i] + 1;
    938:         nz = ai[i + 1] - diag[i] - 1;

    940:         if (T) { /* FIXME: This branch untested */
    941:           PetscCall(PetscArraycpy(w, xb + i2, bs));
    942:           /* copy all rows of x that are needed into contiguous space */
    943:           workt = work;
    944:           for (j = 0; j < nz; j++) {
    945:             PetscCall(PetscArraycpy(workt, x + bs * (*vi++), bs));
    946:             workt += bs;
    947:           }
    948:           arrt = arr;
    949:           for (j = 0; j < nz; j++) {
    950:             PetscCall(PetscArraycpy(arrt, T, bs2));
    951:             for (k = 0; k < bs2; k++) arrt[k] *= v[j];
    952:             arrt += bs2;
    953:           }
    954:           PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, arr, work);
    955:         } else if (kaij->isTI) {
    956:           PetscCall(PetscArraycpy(w, t + i2, bs));
    957:           for (j = 0; j < nz; j++) {
    958:             for (k = 0; k < bs; k++) w[k] -= v[j] * x[vi[j] * bs + k];
    959:           }
    960:         }

    962:         PetscKernel_w_gets_Ar_times_v(bs, bs, w, idiag, y); /* RHS incorrect for omega != 1.0 */
    963:         for (j = 0; j < bs; j++) x[i2 + j] = (1.0 - omega) * x[i2 + j] + omega * y[j];

    965:         idiag -= bs2;
    966:         i2 -= bs;
    967:       }
    968:       PetscCall(PetscLogFlops(1.0 * bs2 * (a->nz)));
    969:     }
    970:     its--;
    971:   }
    972:   while (its--) { /* FIXME: This branch not updated */
    973:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
    974:       i2    = 0;
    975:       idiag = kaij->ibdiag;
    976:       for (i = 0; i < m; i++) {
    977:         PetscCall(PetscArraycpy(w, b + i2, bs));

    979:         v     = aa + ai[i];
    980:         vi    = aj + ai[i];
    981:         nz    = diag[i] - ai[i];
    982:         workt = work;
    983:         for (j = 0; j < nz; j++) {
    984:           PetscCall(PetscArraycpy(workt, x + bs * (*vi++), bs));
    985:           workt += bs;
    986:         }
    987:         arrt = arr;
    988:         if (T) {
    989:           for (j = 0; j < nz; j++) {
    990:             PetscCall(PetscArraycpy(arrt, T, bs2));
    991:             for (k = 0; k < bs2; k++) arrt[k] *= v[j];
    992:             arrt += bs2;
    993:           }
    994:           PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, arr, work);
    995:         } else if (kaij->isTI) {
    996:           for (j = 0; j < nz; j++) {
    997:             PetscCall(PetscArrayzero(arrt, bs2));
    998:             for (k = 0; k < bs; k++) arrt[k + bs * k] = v[j];
    999:             arrt += bs2;
    1000:           }
    1001:           PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, arr, work);
    1002:         }
    1003:         PetscCall(PetscArraycpy(t + i2, w, bs));

    1005:         v     = aa + diag[i] + 1;
    1006:         vi    = aj + diag[i] + 1;
    1007:         nz    = ai[i + 1] - diag[i] - 1;
    1008:         workt = work;
    1009:         for (j = 0; j < nz; j++) {
    1010:           PetscCall(PetscArraycpy(workt, x + bs * (*vi++), bs));
    1011:           workt += bs;
    1012:         }
    1013:         arrt = arr;
    1014:         if (T) {
    1015:           for (j = 0; j < nz; j++) {
    1016:             PetscCall(PetscArraycpy(arrt, T, bs2));
    1017:             for (k = 0; k < bs2; k++) arrt[k] *= v[j];
    1018:             arrt += bs2;
    1019:           }
    1020:           PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, arr, work);
    1021:         } else if (kaij->isTI) {
    1022:           for (j = 0; j < nz; j++) {
    1023:             PetscCall(PetscArrayzero(arrt, bs2));
    1024:             for (k = 0; k < bs; k++) arrt[k + bs * k] = v[j];
    1025:             arrt += bs2;
    1026:           }
    1027:           PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, arr, work);
    1028:         }

    1030:         PetscKernel_w_gets_Ar_times_v(bs, bs, w, idiag, y);
    1031:         for (j = 0; j < bs; j++) *(x + i2 + j) = (1.0 - omega) * *(x + i2 + j) + omega * *(y + j);

    1033:         idiag += bs2;
    1034:         i2 += bs;
    1035:       }
    1036:       xb = t;
    1037:     } else xb = b;
    1038:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
    1039:       idiag = kaij->ibdiag + bs2 * (m - 1);
    1040:       i2    = bs * (m - 1);
    1041:       if (xb == b) {
    1042:         for (i = m - 1; i >= 0; i--) {
    1043:           PetscCall(PetscArraycpy(w, b + i2, bs));

    1045:           v     = aa + ai[i];
    1046:           vi    = aj + ai[i];
    1047:           nz    = diag[i] - ai[i];
    1048:           workt = work;
    1049:           for (j = 0; j < nz; j++) {
    1050:             PetscCall(PetscArraycpy(workt, x + bs * (*vi++), bs));
    1051:             workt += bs;
    1052:           }
    1053:           arrt = arr;
    1054:           if (T) {
    1055:             for (j = 0; j < nz; j++) {
    1056:               PetscCall(PetscArraycpy(arrt, T, bs2));
    1057:               for (k = 0; k < bs2; k++) arrt[k] *= v[j];
    1058:               arrt += bs2;
    1059:             }
    1060:             PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, arr, work);
    1061:           } else if (kaij->isTI) {
    1062:             for (j = 0; j < nz; j++) {
    1063:               PetscCall(PetscArrayzero(arrt, bs2));
    1064:               for (k = 0; k < bs; k++) arrt[k + bs * k] = v[j];
    1065:               arrt += bs2;
    1066:             }
    1067:             PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, arr, work);
    1068:           }

    1070:           v     = aa + diag[i] + 1;
    1071:           vi    = aj + diag[i] + 1;
    1072:           nz    = ai[i + 1] - diag[i] - 1;
    1073:           workt = work;
    1074:           for (j = 0; j < nz; j++) {
    1075:             PetscCall(PetscArraycpy(workt, x + bs * (*vi++), bs));
    1076:             workt += bs;
    1077:           }
    1078:           arrt = arr;
    1079:           if (T) {
    1080:             for (j = 0; j < nz; j++) {
    1081:               PetscCall(PetscArraycpy(arrt, T, bs2));
    1082:               for (k = 0; k < bs2; k++) arrt[k] *= v[j];
    1083:               arrt += bs2;
    1084:             }
    1085:             PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, arr, work);
    1086:           } else if (kaij->isTI) {
    1087:             for (j = 0; j < nz; j++) {
    1088:               PetscCall(PetscArrayzero(arrt, bs2));
    1089:               for (k = 0; k < bs; k++) arrt[k + bs * k] = v[j];
    1090:               arrt += bs2;
    1091:             }
    1092:             PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, arr, work);
    1093:           }

    1095:           PetscKernel_w_gets_Ar_times_v(bs, bs, w, idiag, y);
    1096:           for (j = 0; j < bs; j++) *(x + i2 + j) = (1.0 - omega) * *(x + i2 + j) + omega * *(y + j);
    1097:         }
    1098:       } else {
    1099:         for (i = m - 1; i >= 0; i--) {
    1100:           PetscCall(PetscArraycpy(w, xb + i2, bs));
    1101:           v     = aa + diag[i] + 1;
    1102:           vi    = aj + diag[i] + 1;
    1103:           nz    = ai[i + 1] - diag[i] - 1;
    1104:           workt = work;
    1105:           for (j = 0; j < nz; j++) {
    1106:             PetscCall(PetscArraycpy(workt, x + bs * (*vi++), bs));
    1107:             workt += bs;
    1108:           }
    1109:           arrt = arr;
    1110:           if (T) {
    1111:             for (j = 0; j < nz; j++) {
    1112:               PetscCall(PetscArraycpy(arrt, T, bs2));
    1113:               for (k = 0; k < bs2; k++) arrt[k] *= v[j];
    1114:               arrt += bs2;
    1115:             }
    1116:             PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, arr, work);
    1117:           } else if (kaij->isTI) {
    1118:             for (j = 0; j < nz; j++) {
    1119:               PetscCall(PetscArrayzero(arrt, bs2));
    1120:               for (k = 0; k < bs; k++) arrt[k + bs * k] = v[j];
    1121:               arrt += bs2;
    1122:             }
    1123:             PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, arr, work);
    1124:           }
    1125:           PetscKernel_w_gets_Ar_times_v(bs, bs, w, idiag, y);
    1126:           for (j = 0; j < bs; j++) *(x + i2 + j) = (1.0 - omega) * *(x + i2 + j) + omega * *(y + j);
    1127:         }
    1128:       }
    1129:       PetscCall(PetscLogFlops(1.0 * bs2 * (a->nz)));
    1130:     }
    1131:   }

    1133:   PetscCall(VecRestoreArray(xx, &x));
    1134:   PetscCall(VecRestoreArrayRead(bb, &b));
    1135:   PetscFunctionReturn(PETSC_SUCCESS);
    1136: }

    1138: /*===================================================================================*/

    1140: static PetscErrorCode MatMultAdd_MPIKAIJ(Mat A, Vec xx, Vec yy, Vec zz)
    1141: {
    1142:   Mat_MPIKAIJ *b = (Mat_MPIKAIJ *)A->data;

    1144:   PetscFunctionBegin;
    1145:   if (!yy) PetscCall(VecSet(zz, 0.0));
    1146:   else PetscCall(VecCopy(yy, zz));
    1147:   PetscCall(MatKAIJ_build_AIJ_OAIJ(A)); /* Ensure b->AIJ and b->OAIJ are up to date. */
    1148:   /* start the scatter */
    1149:   PetscCall(VecScatterBegin(b->ctx, xx, b->w, INSERT_VALUES, SCATTER_FORWARD));
    1150:   PetscCall((*b->AIJ->ops->multadd)(b->AIJ, xx, zz, zz));
    1151:   PetscCall(VecScatterEnd(b->ctx, xx, b->w, INSERT_VALUES, SCATTER_FORWARD));
    1152:   PetscCall((*b->OAIJ->ops->multadd)(b->OAIJ, b->w, zz, zz));
    1153:   PetscFunctionReturn(PETSC_SUCCESS);
    1154: }

    1156: static PetscErrorCode MatMult_MPIKAIJ(Mat A, Vec xx, Vec yy)
    1157: {
    1158:   PetscFunctionBegin;
    1159:   PetscCall(MatMultAdd_MPIKAIJ(A, xx, NULL, yy));
    1160:   PetscFunctionReturn(PETSC_SUCCESS);
    1161: }

    1163: static PetscErrorCode MatInvertBlockDiagonal_MPIKAIJ(Mat A, const PetscScalar **values)
    1164: {
    1165:   Mat_MPIKAIJ *b = (Mat_MPIKAIJ *)A->data;

    1167:   PetscFunctionBegin;
    1168:   PetscCall(MatKAIJ_build_AIJ_OAIJ(A)); /* Ensure b->AIJ is up to date. */
    1169:   PetscCall((*b->AIJ->ops->invertblockdiagonal)(b->AIJ, values));
    1170:   PetscFunctionReturn(PETSC_SUCCESS);
    1171: }

    1173: static PetscErrorCode MatGetRow_SeqKAIJ(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **values)
    1174: {
    1175:   Mat_SeqKAIJ *b    = (Mat_SeqKAIJ *)A->data;
    1176:   PetscBool    diag = PETSC_FALSE;
    1177:   PetscInt     nzaij, nz, *colsaij, *idx, i, j, p = b->p, q = b->q, r = row / p, s = row % p, c;
    1178:   PetscScalar *vaij, *v, *S = b->S, *T = b->T;

    1180:   PetscFunctionBegin;
    1181:   PetscCheck(!b->getrowactive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Already active");
    1182:   b->getrowactive = PETSC_TRUE;
    1183:   PetscCheck(row >= 0 && row < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " out of range", row);

    1185:   if ((!S) && (!T) && (!b->isTI)) {
    1186:     if (ncols) *ncols = 0;
    1187:     if (cols) *cols = NULL;
    1188:     if (values) *values = NULL;
    1189:     PetscFunctionReturn(PETSC_SUCCESS);
    1190:   }

    1192:   if (T || b->isTI) {
    1193:     PetscCall(MatGetRow_SeqAIJ(b->AIJ, r, &nzaij, &colsaij, &vaij));
    1194:     c = nzaij;
    1195:     for (i = 0; i < nzaij; i++) {
    1196:       /* check if this row contains a diagonal entry */
    1197:       if (colsaij[i] == r) {
    1198:         diag = PETSC_TRUE;
    1199:         c    = i;
    1200:       }
    1201:     }
    1202:   } else nzaij = c = 0;

    1204:   /* calculate size of row */
    1205:   nz = 0;
    1206:   if (S) nz += q;
    1207:   if (T || b->isTI) nz += (diag && S ? (nzaij - 1) * q : nzaij * q);

    1209:   if (cols || values) {
    1210:     PetscCall(PetscMalloc2(nz, &idx, nz, &v));
    1211:     for (i = 0; i < q; i++) {
    1212:       /* We need to initialize the v[i] to zero to handle the case in which T is NULL (not the identity matrix). */
    1213:       v[i] = 0.0;
    1214:     }
    1215:     if (b->isTI) {
    1216:       for (i = 0; i < nzaij; i++) {
    1217:         for (j = 0; j < q; j++) {
    1218:           idx[i * q + j] = colsaij[i] * q + j;
    1219:           v[i * q + j]   = (j == s ? vaij[i] : 0);
    1220:         }
    1221:       }
    1222:     } else if (T) {
    1223:       for (i = 0; i < nzaij; i++) {
    1224:         for (j = 0; j < q; j++) {
    1225:           idx[i * q + j] = colsaij[i] * q + j;
    1226:           v[i * q + j]   = vaij[i] * T[s + j * p];
    1227:         }
    1228:       }
    1229:     }
    1230:     if (S) {
    1231:       for (j = 0; j < q; j++) {
    1232:         idx[c * q + j] = r * q + j;
    1233:         v[c * q + j] += S[s + j * p];
    1234:       }
    1235:     }
    1236:   }

    1238:   if (ncols) *ncols = nz;
    1239:   if (cols) *cols = idx;
    1240:   if (values) *values = v;
    1241:   PetscFunctionReturn(PETSC_SUCCESS);
    1242: }

    1244: static PetscErrorCode MatRestoreRow_SeqKAIJ(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
    1245: {
    1246:   PetscFunctionBegin;
    1247:   PetscCall(PetscFree2(*idx, *v));
    1248:   ((Mat_SeqKAIJ *)A->data)->getrowactive = PETSC_FALSE;
    1249:   PetscFunctionReturn(PETSC_SUCCESS);
    1250: }

    1252: static PetscErrorCode MatGetRow_MPIKAIJ(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **values)
    1253: {
    1254:   Mat_MPIKAIJ   *b    = (Mat_MPIKAIJ *)A->data;
    1255:   Mat            AIJ  = b->A;
    1256:   PetscBool      diag = PETSC_FALSE;
    1257:   Mat            MatAIJ, MatOAIJ;
    1258:   const PetscInt rstart = A->rmap->rstart, rend = A->rmap->rend, p = b->p, q = b->q, *garray;
    1259:   PetscInt       nz, *idx, ncolsaij = 0, ncolsoaij = 0, *colsaij, *colsoaij, r, s, c, i, j, lrow;
    1260:   PetscScalar   *v, *vals, *ovals, *S = b->S, *T = b->T;

    1262:   PetscFunctionBegin;
    1263:   PetscCall(MatKAIJ_build_AIJ_OAIJ(A)); /* Ensure b->AIJ and b->OAIJ are up to date. */
    1264:   MatAIJ  = ((Mat_SeqKAIJ *)b->AIJ->data)->AIJ;
    1265:   MatOAIJ = ((Mat_SeqKAIJ *)b->OAIJ->data)->AIJ;
    1266:   PetscCheck(!b->getrowactive, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Already active");
    1267:   b->getrowactive = PETSC_TRUE;
    1268:   PetscCheck(row >= rstart && row < rend, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Only local rows");
    1269:   lrow = row - rstart;

    1271:   if ((!S) && (!T) && (!b->isTI)) {
    1272:     if (ncols) *ncols = 0;
    1273:     if (cols) *cols = NULL;
    1274:     if (values) *values = NULL;
    1275:     PetscFunctionReturn(PETSC_SUCCESS);
    1276:   }

    1278:   r = lrow / p;
    1279:   s = lrow % p;

    1281:   if (T || b->isTI) {
    1282:     PetscCall(MatMPIAIJGetSeqAIJ(AIJ, NULL, NULL, &garray));
    1283:     PetscCall(MatGetRow_SeqAIJ(MatAIJ, lrow / p, &ncolsaij, &colsaij, &vals));
    1284:     PetscCall(MatGetRow_SeqAIJ(MatOAIJ, lrow / p, &ncolsoaij, &colsoaij, &ovals));

    1286:     c = ncolsaij + ncolsoaij;
    1287:     for (i = 0; i < ncolsaij; i++) {
    1288:       /* check if this row contains a diagonal entry */
    1289:       if (colsaij[i] == r) {
    1290:         diag = PETSC_TRUE;
    1291:         c    = i;
    1292:       }
    1293:     }
    1294:   } else c = 0;

    1296:   /* calculate size of row */
    1297:   nz = 0;
    1298:   if (S) nz += q;
    1299:   if (T || b->isTI) nz += (diag && S ? (ncolsaij + ncolsoaij - 1) * q : (ncolsaij + ncolsoaij) * q);

    1301:   if (cols || values) {
    1302:     PetscCall(PetscMalloc2(nz, &idx, nz, &v));
    1303:     for (i = 0; i < q; i++) {
    1304:       /* We need to initialize the v[i] to zero to handle the case in which T is NULL (not the identity matrix). */
    1305:       v[i] = 0.0;
    1306:     }
    1307:     if (b->isTI) {
    1308:       for (i = 0; i < ncolsaij; i++) {
    1309:         for (j = 0; j < q; j++) {
    1310:           idx[i * q + j] = (colsaij[i] + rstart / p) * q + j;
    1311:           v[i * q + j]   = (j == s ? vals[i] : 0.0);
    1312:         }
    1313:       }
    1314:       for (i = 0; i < ncolsoaij; i++) {
    1315:         for (j = 0; j < q; j++) {
    1316:           idx[(i + ncolsaij) * q + j] = garray[colsoaij[i]] * q + j;
    1317:           v[(i + ncolsaij) * q + j]   = (j == s ? ovals[i] : 0.0);
    1318:         }
    1319:       }
    1320:     } else if (T) {
    1321:       for (i = 0; i < ncolsaij; i++) {
    1322:         for (j = 0; j < q; j++) {
    1323:           idx[i * q + j] = (colsaij[i] + rstart / p) * q + j;
    1324:           v[i * q + j]   = vals[i] * T[s + j * p];
    1325:         }
    1326:       }
    1327:       for (i = 0; i < ncolsoaij; i++) {
    1328:         for (j = 0; j < q; j++) {
    1329:           idx[(i + ncolsaij) * q + j] = garray[colsoaij[i]] * q + j;
    1330:           v[(i + ncolsaij) * q + j]   = ovals[i] * T[s + j * p];
    1331:         }
    1332:       }
    1333:     }
    1334:     if (S) {
    1335:       for (j = 0; j < q; j++) {
    1336:         idx[c * q + j] = (r + rstart / p) * q + j;
    1337:         v[c * q + j] += S[s + j * p];
    1338:       }
    1339:     }
    1340:   }

    1342:   if (ncols) *ncols = nz;
    1343:   if (cols) *cols = idx;
    1344:   if (values) *values = v;
    1345:   PetscFunctionReturn(PETSC_SUCCESS);
    1346: }

    1348: static PetscErrorCode MatRestoreRow_MPIKAIJ(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
    1349: {
    1350:   PetscFunctionBegin;
    1351:   PetscCall(PetscFree2(*idx, *v));
    1352:   ((Mat_SeqKAIJ *)A->data)->getrowactive = PETSC_FALSE;
    1353:   PetscFunctionReturn(PETSC_SUCCESS);
    1354: }

    1356: static PetscErrorCode MatCreateSubMatrix_KAIJ(Mat mat, IS isrow, IS iscol, MatReuse cll, Mat *newmat)
    1357: {
    1358:   Mat A;

    1360:   PetscFunctionBegin;
    1361:   PetscCall(MatConvert(mat, MATAIJ, MAT_INITIAL_MATRIX, &A));
    1362:   PetscCall(MatCreateSubMatrix(A, isrow, iscol, cll, newmat));
    1363:   PetscCall(MatDestroy(&A));
    1364:   PetscFunctionReturn(PETSC_SUCCESS);
    1365: }

    1367: /*@C
    1368:   MatCreateKAIJ - Creates a matrix of type `MATKAIJ`.

    1370:   Collective

    1372:   Input Parameters:
    1373: + A - the `MATAIJ` matrix
    1374: . p - number of rows in `S` and `T`
    1375: . q - number of columns in `S` and `T`
    1376: . S - the `S` matrix (can be `NULL`), stored as a `PetscScalar` array (column-major)
    1377: - T - the `T` matrix (can be `NULL`), stored as a `PetscScalar` array (column-major)

    1379:   Output Parameter:
    1380: . kaij - the new `MATKAIJ` matrix

    1382:   Level: advanced

    1384:   Notes:
    1385:   The created matrix is of the following form\:
    1386: .vb
    1387:     [I \otimes S + A \otimes T]
    1388: .ve
    1389:   where
    1390: .vb
    1391:   S is a dense (p \times q) matrix
    1392:   T is a dense (p \times q) matrix
    1393:   A is a `MATAIJ`  (n \times n) matrix
    1394:   I is the identity matrix
    1395: .ve
    1396:   The resulting matrix is (np \times nq)

    1398:   `S` and `T` are always stored independently on all processes as `PetscScalar` arrays in
    1399:   column-major format.

    1401:   This function increases the reference count on the `MATAIJ` matrix, so the user is free to destroy the matrix if it is not needed.

    1403:   Changes to the entries of the `MATAIJ` matrix will immediately affect the `MATKAIJ` matrix.

    1405:   Developer Notes:
    1406:   In the `MATMPIKAIJ` case, the internal 'AIJ' and 'OAIJ' sequential KAIJ matrices are kept up to date by tracking the object state
    1407:   of the AIJ matrix 'A' that describes the blockwise action of the `MATMPIKAIJ` matrix and, if the object state has changed, lazily
    1408:   rebuilding 'AIJ' and 'OAIJ' just before executing operations with the `MATMPIKAIJ` matrix. If new types of operations are added,
    1409:   routines implementing those must also ensure these are rebuilt when needed (by calling the internal MatKAIJ_build_AIJ_OAIJ() routine).

    1411: .seealso: [](ch_matrices), `Mat`, `MatKAIJSetAIJ()`, `MatKAIJSetS()`, `MatKAIJSetT()`, `MatKAIJGetAIJ()`, `MatKAIJGetS()`, `MatKAIJGetT()`, `MATKAIJ`
    1412: @*/
    1413: PetscErrorCode MatCreateKAIJ(Mat A, PetscInt p, PetscInt q, const PetscScalar S[], const PetscScalar T[], Mat *kaij)
    1414: {
    1415:   PetscFunctionBegin;
    1416:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), kaij));
    1417:   PetscCall(MatSetType(*kaij, MATKAIJ));
    1418:   PetscCall(MatKAIJSetAIJ(*kaij, A));
    1419:   PetscCall(MatKAIJSetS(*kaij, p, q, S));
    1420:   PetscCall(MatKAIJSetT(*kaij, p, q, T));
    1421:   PetscCall(MatSetUp(*kaij));
    1422:   PetscFunctionReturn(PETSC_SUCCESS);
    1423: }

    1425: /*MC
    1426:   MATKAIJ - MATKAIJ = "kaij" - A matrix type to be used to evaluate matrices of form
    1427:     [I \otimes S + A \otimes T],
    1428:   where
    1429: .vb
    1430:     S is a dense (p \times q) matrix,
    1431:     T is a dense (p \times q) matrix,
    1432:     A is an AIJ  (n \times n) matrix,
    1433:     and I is the identity matrix.
    1434: .ve
    1435:   The resulting matrix is (np \times nq).

    1437:   S and T are always stored independently on all processes as `PetscScalar` arrays in column-major format.

    1439:   Level: advanced

    1441:   Note:
    1442:   A linear system with multiple right-hand sides, AX = B, can be expressed in the KAIJ-friendly form of (A \otimes I) x = b,
    1443:   where x and b are column vectors containing the row-major representations of X and B.

    1445: .seealso: [](ch_matrices), `Mat`, `MatKAIJSetAIJ()`, `MatKAIJSetS()`, `MatKAIJSetT()`, `MatKAIJGetAIJ()`, `MatKAIJGetS()`, `MatKAIJGetT()`, `MatCreateKAIJ()`
    1446: M*/

    1448: PETSC_EXTERN PetscErrorCode MatCreate_KAIJ(Mat A)
    1449: {
    1450:   Mat_MPIKAIJ *b;
    1451:   PetscMPIInt  size;

    1453:   PetscFunctionBegin;
    1454:   PetscCall(PetscNew(&b));
    1455:   A->data = (void *)b;

    1457:   PetscCall(PetscMemzero(A->ops, sizeof(struct _MatOps)));

    1459:   b->w = NULL;
    1460:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
    1461:   if (size == 1) {
    1462:     PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATSEQKAIJ));
    1463:     A->ops->destroy             = MatDestroy_SeqKAIJ;
    1464:     A->ops->mult                = MatMult_SeqKAIJ;
    1465:     A->ops->multadd             = MatMultAdd_SeqKAIJ;
    1466:     A->ops->invertblockdiagonal = MatInvertBlockDiagonal_SeqKAIJ;
    1467:     A->ops->getrow              = MatGetRow_SeqKAIJ;
    1468:     A->ops->restorerow          = MatRestoreRow_SeqKAIJ;
    1469:     A->ops->sor                 = MatSOR_SeqKAIJ;
    1470:     PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqkaij_seqaij_C", MatConvert_KAIJ_AIJ));
    1471:   } else {
    1472:     PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATMPIKAIJ));
    1473:     A->ops->destroy             = MatDestroy_MPIKAIJ;
    1474:     A->ops->mult                = MatMult_MPIKAIJ;
    1475:     A->ops->multadd             = MatMultAdd_MPIKAIJ;
    1476:     A->ops->invertblockdiagonal = MatInvertBlockDiagonal_MPIKAIJ;
    1477:     A->ops->getrow              = MatGetRow_MPIKAIJ;
    1478:     A->ops->restorerow          = MatRestoreRow_MPIKAIJ;
    1479:     PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatGetDiagonalBlock_C", MatGetDiagonalBlock_MPIKAIJ));
    1480:     PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_mpikaij_mpiaij_C", MatConvert_KAIJ_AIJ));
    1481:   }
    1482:   A->ops->setup           = MatSetUp_KAIJ;
    1483:   A->ops->view            = MatView_KAIJ;
    1484:   A->ops->createsubmatrix = MatCreateSubMatrix_KAIJ;
    1485:   PetscFunctionReturn(PETSC_SUCCESS);
    1486: }