Actual source code: ex81.c

      1: #include <petsc.h>

      3: static char help[] = "Solves a linear system with a MatNest and nested fields.\n\n";

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

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

     19:   PetscFunctionBeginUser;
     20:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
     21:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
     22:   PetscCall(PetscMalloc1(found, &outer));
     23:   PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-outer_fieldsplit_sizes", outer, &found, &flg));
     24:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
     25:   if (flg) {
     26:     PetscCheck(found != 1, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must supply more than one field");
     27:     j = 0;
     28:     for (i = 0; i < found; ++i) j += outer[i];
     29:     PetscCheck(j == Q, PETSC_COMM_WORLD, PETSC_ERR_USER, "Sum of outer fieldsplit sizes (%" PetscInt_FMT ") greater than number of blocks in MatNest (%d)", j, Q);
     30:   }
     31:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
     32:   size = PetscMax(3, size);
     33:   for (i = 0; i < Q * Q; ++i) array[i] = NULL;
     34:   for (i = 0; i < Q; ++i) {
     35:     if (i == 0) {
     36:       PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, size, size, 1, NULL, 0, NULL, array + (Q + 1) * i));
     37:     } else if (i == 1 || i == 3) {
     38:       PetscCall(MatCreateSBAIJ(PETSC_COMM_WORLD, 2, PETSC_DECIDE, PETSC_DECIDE, size, size, 1, NULL, 0, NULL, array + (Q + 1) * i));
     39:     } else if (i == 2 || i == 4) {
     40:       PetscCall(MatCreateBAIJ(PETSC_COMM_WORLD, 2, PETSC_DECIDE, PETSC_DECIDE, size, size, 1, NULL, 0, NULL, array + (Q + 1) * i));
     41:     }
     42:     PetscCall(MatAssemblyBegin(array[(Q + 1) * i], MAT_FINAL_ASSEMBLY));
     43:     PetscCall(MatAssemblyEnd(array[(Q + 1) * i], MAT_FINAL_ASSEMBLY));
     44:     PetscCall(MatShift(array[(Q + 1) * i], 100 + i + 1));
     45:     if (i == 3) {
     46:       PetscCall(MatDuplicate(array[(Q + 1) * i], MAT_COPY_VALUES, &a));
     47:       PetscCall(MatDestroy(array + (Q + 1) * i));
     48:       PetscCall(MatCreateHermitianTranspose(a, array + (Q + 1) * i));
     49:       PetscCall(MatDestroy(&a));
     50:     }
     51:     size *= 2;
     52:   }
     53:   PetscCall(MatGetSize(array[0], &M, NULL));
     54:   for (i = 2; i < Q; ++i) {
     55:     PetscCall(MatGetSize(array[(Q + 1) * i], NULL, &N));
     56:     if (i != Q - 1) {
     57:       PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, i == 3 ? N : M, i == 3 ? M : N, 0, NULL, 0, NULL, array + i));
     58:     } else {
     59:       PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, NULL, array + i));
     60:     }
     61:     PetscCall(MatAssemblyBegin(array[i], MAT_FINAL_ASSEMBLY));
     62:     PetscCall(MatAssemblyEnd(array[i], MAT_FINAL_ASSEMBLY));
     63:     PetscCall(MatSetRandom(array[i], rctx));
     64:     if (i == 3) {
     65:       PetscCall(MatDuplicate(array[i], MAT_COPY_VALUES, &a));
     66:       PetscCall(MatDestroy(array + i));
     67:       PetscCall(MatCreateHermitianTranspose(a, array + i));
     68:       PetscCall(MatDestroy(&a));
     69:     }
     70:   }
     71:   PetscCall(MatGetSize(array[0], NULL, &N));
     72:   for (i = 2; i < Q; i += 2) {
     73:     PetscCall(MatGetSize(array[(Q + 1) * i], &M, NULL));
     74:     if (i != Q - 1) {
     75:       PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, 2, NULL, 2, NULL, array + Q * i));
     76:     } else {
     77:       PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, N, M, NULL, array + Q * i));
     78:     }
     79:     PetscCall(MatAssemblyBegin(array[Q * i], MAT_FINAL_ASSEMBLY));
     80:     PetscCall(MatAssemblyEnd(array[Q * i], MAT_FINAL_ASSEMBLY));
     81:     PetscCall(MatSetRandom(array[Q * i], rctx));
     82:     if (i == Q - 1) {
     83:       PetscCall(MatDuplicate(array[Q * i], MAT_COPY_VALUES, &a));
     84:       PetscCall(MatDestroy(array + Q * i));
     85:       PetscCall(MatCreateHermitianTranspose(a, array + Q * i));
     86:       PetscCall(MatDestroy(&a));
     87:     }
     88:   }
     89:   PetscCall(MatGetSize(array[(Q + 1) * 3], &M, NULL));
     90:   for (i = 1; i < 3; ++i) {
     91:     PetscCall(MatGetSize(array[(Q + 1) * i], NULL, &N));
     92:     PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, 2, NULL, 2, NULL, array + Q * 3 + i));
     93:     PetscCall(MatAssemblyBegin(array[Q * 3 + i], MAT_FINAL_ASSEMBLY));
     94:     PetscCall(MatAssemblyEnd(array[Q * 3 + i], MAT_FINAL_ASSEMBLY));
     95:     PetscCall(MatSetRandom(array[Q * 3 + i], rctx));
     96:   }
     97:   PetscCall(MatGetSize(array[(Q + 1) * 1], NULL, &N));
     98:   PetscCall(MatGetSize(array[(Q + 1) * (Q - 1)], &M, NULL));
     99:   PetscCall(MatCreateBAIJ(PETSC_COMM_WORLD, 2, PETSC_DECIDE, PETSC_DECIDE, M, N, 0, NULL, 0, NULL, &a));
    100:   PetscCall(MatAssemblyBegin(a, MAT_FINAL_ASSEMBLY));
    101:   PetscCall(MatAssemblyEnd(a, MAT_FINAL_ASSEMBLY));
    102:   PetscCall(MatCreateHermitianTranspose(a, array + Q + Q - 1));
    103:   PetscCall(MatDestroy(&a));
    104:   PetscCall(MatDestroy(array + Q * Q - 1));
    105:   PetscCall(MatCreateNest(PETSC_COMM_WORLD, Q, NULL, Q, NULL, array, &A));
    106:   for (i = 0; i < Q; ++i) PetscCall(MatDestroy(array + (Q + 1) * i));
    107:   for (i = 2; i < Q; ++i) {
    108:     PetscCall(MatDestroy(array + i));
    109:     PetscCall(MatDestroy(array + Q * i));
    110:   }
    111:   for (i = 1; i < 3; ++i) PetscCall(MatDestroy(array + Q * 3 + i));
    112:   PetscCall(MatDestroy(array + Q + Q - 1));
    113:   PetscCall(KSPSetOperators(ksp, A, A));
    114:   PetscCall(MatNestGetISs(A, rows, NULL));
    115:   PetscCall(KSPGetPC(ksp, &pc));
    116:   PetscCall(PCSetType(pc, PCFIELDSPLIT));
    117:   M = 0;
    118:   for (j = 0; j < found; ++j) {
    119:     IS expand1, expand2;
    120:     PetscCall(ISDuplicate(rows[M], &expand1));
    121:     for (i = 1; i < outer[j]; ++i) {
    122:       PetscCall(ISExpand(expand1, rows[M + i], &expand2));
    123:       PetscCall(ISDestroy(&expand1));
    124:       expand1 = expand2;
    125:     }
    126:     M += outer[j];
    127:     PetscCall(PCFieldSplitSetIS(pc, NULL, expand1));
    128:     PetscCall(ISDestroy(&expand1));
    129:   }
    130:   PetscCall(KSPSetFromOptions(ksp));
    131:   flg = PETSC_FALSE;
    132:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matmult", &flg, NULL));
    133:   if (flg) {
    134:     Mat D, E, F, H;
    135:     PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &D));
    136:     PetscCall(MatMultEqual(A, D, 10, &flg));
    137:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatDense != MatNest");
    138:     PetscCall(MatZeroEntries(D));
    139:     PetscCall(MatConvert(A, MATDENSE, MAT_REUSE_MATRIX, &D));
    140:     PetscCall(MatMultEqual(A, D, 10, &flg));
    141:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatDense != MatNest");
    142:     PetscCall(MatConvert(D, MATAIJ, MAT_INITIAL_MATRIX, &E));
    143:     PetscCall(MatMultEqual(E, D, 10, &flg));
    144:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatDense != MatAIJ");
    145:     PetscCall(MatZeroEntries(E));
    146:     PetscCall(MatConvert(D, MATAIJ, MAT_REUSE_MATRIX, &E));
    147:     PetscCall(MatMultEqual(E, D, 10, &flg));
    148:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatDense != MatAIJ");
    149:     PetscCall(MatConvert(E, MATDENSE, MAT_INPLACE_MATRIX, &E));
    150:     PetscCall(MatMultEqual(A, E, 10, &flg));
    151:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatAIJ != MatNest");
    152:     PetscCall(MatConvert(D, MATAIJ, MAT_INPLACE_MATRIX, &D));
    153:     PetscCall(MatMultEqual(A, D, 10, &flg));
    154:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatDense != MatNest");
    155:     PetscCall(MatDestroy(&E));
    156:     PetscCall(MatCreateHermitianTranspose(D, &E));
    157:     PetscCall(MatConvert(E, MATAIJ, MAT_INITIAL_MATRIX, &F));
    158:     PetscCall(MatMultEqual(E, F, 10, &flg));
    159:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatHermitianTranspose != MatAIJ");
    160:     PetscCall(MatZeroEntries(F));
    161:     PetscCall(MatConvert(E, MATAIJ, MAT_REUSE_MATRIX, &F));
    162:     PetscCall(MatMultEqual(E, F, 10, &flg));
    163:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatHermitianTranspose != MatAIJ");
    164:     PetscCall(MatDestroy(&F));
    165:     PetscCall(MatConvert(E, MATAIJ, MAT_INPLACE_MATRIX, &E));
    166:     PetscCall(MatCreateHermitianTranspose(D, &H));
    167:     PetscCall(MatMultEqual(E, H, 10, &flg));
    168:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatHermitianTranspose != MatAIJ");
    169:     PetscCall(MatDestroy(&H));
    170:     PetscCall(MatDestroy(&E));
    171:     PetscCall(MatDestroy(&D));
    172:   }
    173:   PetscCall(KSPSetUp(ksp));
    174:   PetscCall(MatCreateVecs(A, &b, &x));
    175:   PetscCall(VecSetRandom(b, rctx));
    176:   PetscCall(VecGetSubVector(b, rows[Q - 1], &sub));
    177:   PetscCall(VecSet(sub, 0.0));
    178:   PetscCall(VecRestoreSubVector(b, rows[Q - 1], &sub));
    179:   PetscCall(KSPSolve(ksp, b, x));
    180:   PetscCall(VecDestroy(&b));
    181:   PetscCall(VecDestroy(&x));
    182:   PetscCall(PetscRandomDestroy(&rctx));
    183:   PetscCall(MatDestroy(&A));
    184:   PetscCall(KSPDestroy(&ksp));
    185:   PetscCall(PetscFree(outer));
    186:   PetscCall(PetscFinalize());
    187:   return 0;
    188: }

    190: /*TEST

    192:    test:
    193:       nsize: {{1 3}shared output}
    194:       suffix: wo_explicit_schur
    195:       requires: !complex !single
    196:       filter: sed -e "s/seq/mpi/g" -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" -e "s/ 1 MPI process/ 3 MPI processes/g" -e "s/iterations [4-5]/iterations 4/g" -e "s/hermitiantranspose/transpose/g"
    197:       args: -outer_fieldsplit_sizes {{1,2,2 2,1,2 2,2,1}separate output} -ksp_view_mat -pc_type fieldsplit -ksp_converged_reason -fieldsplit_pc_type jacobi -test_matmult

    199:    test:
    200:       nsize: {{1 3}shared output}
    201:       suffix: w_explicit_schur
    202:       requires: !complex
    203:       filter: sed -e "s/seq/mpi/g" -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" -e "s/ 1 MPI process/ 3 MPI processes/g" -e "s/iterations [1-2]/iterations 1/g" -e "s/hermitiantranspose/transpose/g"
    204:       args: -outer_fieldsplit_sizes {{1,4 2,3 3,2 4,1}separate output} -ksp_view_mat -pc_type fieldsplit -ksp_converged_reason -fieldsplit_pc_type jacobi -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition full

    206: TEST*/