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*/