Actual source code: ex81a.c

      1: #include <petsc.h>

      3: static char help[] = "Solves a linear system with a MatNest and PCFIELDSPLIT with fields defined from the command line.\n\n";
      4: /* similar to ex81.c except the PCFIELDSPLIT fields are defined from the command line with fields instead of hardwired IS from the MATNEST */

      6: #define Q 5 /* everything is hardwired for a 5x5 MatNest for now */

      8: int main(int argc, char **args)
      9: {
     10:   KSP         ksp;
     11:   PC          pc;
     12:   Mat         array[Q * Q], A, a;
     13:   Vec         b, x, sub;
     14:   IS          rows[Q];
     15:   PetscInt    i, M, N;
     16:   PetscMPIInt size;
     17:   PetscRandom rctx;

     19:   PetscFunctionBeginUser;
     20:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
     21:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
     22:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
     23:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
     24:   size = PetscMax(3, size);
     25:   for (i = 0; i < Q * Q; ++i) array[i] = NULL;
     26:   for (i = 0; i < Q; ++i) {
     27:     if (i == 0) {
     28:       PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, size, size, 1, NULL, 0, NULL, array + (Q + 1) * i));
     29:     } else if (i == 1 || i == 3) {
     30:       PetscCall(MatCreateSBAIJ(PETSC_COMM_WORLD, 2, PETSC_DECIDE, PETSC_DECIDE, size, size, 1, NULL, 0, NULL, array + (Q + 1) * i));
     31:     } else if (i == 2 || i == 4) {
     32:       PetscCall(MatCreateBAIJ(PETSC_COMM_WORLD, 2, PETSC_DECIDE, PETSC_DECIDE, size, size, 1, NULL, 0, NULL, array + (Q + 1) * i));
     33:     }
     34:     PetscCall(MatAssemblyBegin(array[(Q + 1) * i], MAT_FINAL_ASSEMBLY));
     35:     PetscCall(MatAssemblyEnd(array[(Q + 1) * i], MAT_FINAL_ASSEMBLY));
     36:     PetscCall(MatShift(array[(Q + 1) * i], 100 + i + 1));
     37:     if (i == 3) {
     38:       PetscCall(MatDuplicate(array[(Q + 1) * i], MAT_COPY_VALUES, &a));
     39:       PetscCall(MatDestroy(array + (Q + 1) * i));
     40:       PetscCall(MatCreateHermitianTranspose(a, array + (Q + 1) * i));
     41:       PetscCall(MatDestroy(&a));
     42:     }
     43:     size *= 2;
     44:   }
     45:   PetscCall(MatGetSize(array[0], &M, NULL));
     46:   for (i = 2; i < Q; ++i) {
     47:     PetscCall(MatGetSize(array[(Q + 1) * i], NULL, &N));
     48:     if (i != Q - 1) {
     49:       PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, i == 3 ? N : M, i == 3 ? M : N, 0, NULL, 0, NULL, array + i));
     50:     } else {
     51:       PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, NULL, array + i));
     52:     }
     53:     PetscCall(MatAssemblyBegin(array[i], MAT_FINAL_ASSEMBLY));
     54:     PetscCall(MatAssemblyEnd(array[i], MAT_FINAL_ASSEMBLY));
     55:     PetscCall(MatSetRandom(array[i], rctx));
     56:     if (i == 3) {
     57:       PetscCall(MatDuplicate(array[i], MAT_COPY_VALUES, &a));
     58:       PetscCall(MatDestroy(array + i));
     59:       PetscCall(MatCreateHermitianTranspose(a, array + i));
     60:       PetscCall(MatDestroy(&a));
     61:     }
     62:   }
     63:   PetscCall(MatGetSize(array[0], NULL, &N));
     64:   for (i = 2; i < Q; i += 2) {
     65:     PetscCall(MatGetSize(array[(Q + 1) * i], &M, NULL));
     66:     if (i != Q - 1) {
     67:       PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, 2, NULL, 2, NULL, array + Q * i));
     68:     } else {
     69:       PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, N, M, NULL, array + Q * i));
     70:     }
     71:     PetscCall(MatAssemblyBegin(array[Q * i], MAT_FINAL_ASSEMBLY));
     72:     PetscCall(MatAssemblyEnd(array[Q * i], MAT_FINAL_ASSEMBLY));
     73:     PetscCall(MatSetRandom(array[Q * i], rctx));
     74:     if (i == Q - 1) {
     75:       PetscCall(MatDuplicate(array[Q * i], MAT_COPY_VALUES, &a));
     76:       PetscCall(MatDestroy(array + Q * i));
     77:       PetscCall(MatCreateHermitianTranspose(a, array + Q * i));
     78:       PetscCall(MatDestroy(&a));
     79:     }
     80:   }
     81:   PetscCall(MatGetSize(array[(Q + 1) * 3], &M, NULL));
     82:   for (i = 1; i < 3; ++i) {
     83:     PetscCall(MatGetSize(array[(Q + 1) * i], NULL, &N));
     84:     PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, 2, NULL, 2, NULL, array + Q * 3 + i));
     85:     PetscCall(MatAssemblyBegin(array[Q * 3 + i], MAT_FINAL_ASSEMBLY));
     86:     PetscCall(MatAssemblyEnd(array[Q * 3 + i], MAT_FINAL_ASSEMBLY));
     87:     PetscCall(MatSetRandom(array[Q * 3 + i], rctx));
     88:   }
     89:   PetscCall(MatGetSize(array[(Q + 1) * 1], NULL, &N));
     90:   PetscCall(MatGetSize(array[(Q + 1) * (Q - 1)], &M, NULL));
     91:   PetscCall(MatCreateBAIJ(PETSC_COMM_WORLD, 2, PETSC_DECIDE, PETSC_DECIDE, M, N, 0, NULL, 0, NULL, &a));
     92:   PetscCall(MatAssemblyBegin(a, MAT_FINAL_ASSEMBLY));
     93:   PetscCall(MatAssemblyEnd(a, MAT_FINAL_ASSEMBLY));
     94:   PetscCall(MatCreateHermitianTranspose(a, array + Q + Q - 1));
     95:   PetscCall(MatDestroy(&a));
     96:   PetscCall(MatDestroy(array + Q * Q - 1));
     97:   PetscCall(MatCreateNest(PETSC_COMM_WORLD, Q, NULL, Q, NULL, array, &A));
     98:   for (i = 0; i < Q; ++i) PetscCall(MatDestroy(array + (Q + 1) * i));
     99:   for (i = 2; i < Q; ++i) {
    100:     PetscCall(MatDestroy(array + i));
    101:     PetscCall(MatDestroy(array + Q * i));
    102:   }
    103:   for (i = 1; i < 3; ++i) PetscCall(MatDestroy(array + Q * 3 + i));
    104:   PetscCall(MatDestroy(array + Q + Q - 1));
    105:   PetscCall(KSPSetOperators(ksp, A, A));
    106:   PetscCall(MatNestGetISs(A, rows, NULL));
    107:   PetscCall(KSPGetPC(ksp, &pc));
    108:   PetscCall(PCSetType(pc, PCFIELDSPLIT));
    109:   PetscCall(KSPSetFromOptions(ksp));
    110:   PetscCall(MatCreateVecs(A, &b, &x));
    111:   PetscCall(VecSetRandom(b, rctx));
    112:   PetscCall(VecGetSubVector(b, rows[Q - 1], &sub));
    113:   PetscCall(VecSet(sub, 0.0));
    114:   PetscCall(VecRestoreSubVector(b, rows[Q - 1], &sub));
    115:   PetscCall(KSPSolve(ksp, b, x));
    116:   PetscCall(VecDestroy(&b));
    117:   PetscCall(VecDestroy(&x));
    118:   PetscCall(PetscRandomDestroy(&rctx));
    119:   PetscCall(MatDestroy(&A));
    120:   PetscCall(KSPDestroy(&ksp));
    121:   PetscCall(PetscFinalize());
    122:   return 0;
    123: }

    125: /*TEST

    127:    test:
    128:       nsize: 3
    129:       suffix: 1
    130:       requires: !complex !single
    131:       filter: sed  -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" -e "s/iterations [4-5]/iterations 4/g" -e "s/hermitiantranspose/transpose/g"
    132:       args: -pc_fieldsplit_0_fields 0,1 -pc_fieldsplit_1_fields 2,3 -pc_fieldsplit_2_fields 4 -pc_type fieldsplit -ksp_converged_reason -fieldsplit_pc_type jacobi -ksp_view

    134:    test:
    135:       suffix: 2
    136:       nsize: 3
    137:       requires: !complex !single
    138:       filter: sed  -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" -e "s/iterations [4-5]/iterations 4/g" -e "s/hermitiantranspose/transpose/g"
    139:       args: -pc_type fieldsplit -ksp_converged_reason -fieldsplit_pc_type jacobi -ksp_view

    141: TEST*/