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;

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

190: /*TEST

192:    test:
193:       nsize: {{1 3}shared output}
194:       suffix: wo_explicit_schur
195:       requires: !complex
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*/