Actual source code: ex178.c
1: static char help[] = "Tests MatInvertVariableBlockEnvelope()\n\n";
3: #include <petscmat.h>
4: extern PetscErrorCode MatIsDiagonal(Mat);
5: extern PetscErrorCode BuildMatrix(const PetscInt *, PetscInt, const PetscInt *, Mat *);
7: int main(int argc, char **argv)
8: {
9: Mat A, C, D, F;
10: PetscInt i, j, rows[2], *parts, cnt, N = 21, nblocks, *blocksizes;
11: PetscScalar values[2][2];
12: PetscReal rand;
13: PetscRandom rctx;
14: PetscMPIInt size;
16: PetscFunctionBeginUser;
17: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
18: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_DENSE));
20: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
21: PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, 6, 18));
22: PetscCall(MatSetFromOptions(C));
23: PetscCall(MatSetUp(C));
24: values[0][0] = 2;
25: values[0][1] = 1;
26: values[1][0] = 1;
27: values[1][1] = 2;
28: for (i = 0; i < 3; i++) {
29: rows[0] = 2 * i;
30: rows[1] = 2 * i + 1;
31: PetscCall(MatSetValues(C, 2, rows, 2, rows, (PetscScalar *)values, INSERT_VALUES));
32: }
33: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
34: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
35: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
37: PetscCall(MatMatTransposeMult(C, C, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &A));
38: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
40: PetscCall(MatInvertVariableBlockEnvelope(A, MAT_INITIAL_MATRIX, &D));
41: PetscCall(MatView(D, PETSC_VIEWER_STDOUT_WORLD));
43: PetscCall(MatMatMult(A, D, MAT_INITIAL_MATRIX, 1.0, &F));
44: PetscCall(MatView(F, PETSC_VIEWER_STDOUT_WORLD));
45: PetscCall(MatIsDiagonal(F));
47: PetscCall(MatDestroy(&A));
48: PetscCall(MatDestroy(&D));
49: PetscCall(MatDestroy(&C));
50: PetscCall(MatDestroy(&F));
52: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rctx));
53: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
54: PetscCall(PetscMalloc1(size, &parts));
56: for (j = 0; j < 128; j++) {
57: cnt = 0;
58: for (i = 0; i < size - 1; i++) {
59: PetscCall(PetscRandomGetValueReal(rctx, &rand));
60: parts[i] = (PetscInt)N * rand;
61: parts[i] = PetscMin(parts[i], N - cnt);
62: cnt += parts[i];
63: }
64: parts[size - 1] = N - cnt;
66: PetscCall(PetscRandomGetValueReal(rctx, &rand));
67: nblocks = rand * 10;
68: nblocks = PetscMax(nblocks, 2);
69: cnt = 0;
70: PetscCall(PetscMalloc1(nblocks, &blocksizes));
71: for (i = 0; i < nblocks - 1; i++) {
72: PetscCall(PetscRandomGetValueReal(rctx, &rand));
73: blocksizes[i] = PetscMax(1, (PetscInt)N * rand);
74: blocksizes[i] = PetscMin(blocksizes[i], N - cnt);
75: cnt += blocksizes[i];
76: if (cnt == N) {
77: nblocks = i + 1;
78: break;
79: }
80: }
81: if (cnt < N) blocksizes[nblocks - 1] = N - cnt;
83: PetscCall(BuildMatrix(parts, nblocks, blocksizes, &A));
84: PetscCall(PetscFree(blocksizes));
86: PetscCall(MatInvertVariableBlockEnvelope(A, MAT_INITIAL_MATRIX, &D));
88: PetscCall(MatMatMult(A, D, MAT_INITIAL_MATRIX, 1.0, &F));
89: PetscCall(MatIsDiagonal(F));
91: PetscCall(MatDestroy(&A));
92: PetscCall(MatDestroy(&D));
93: PetscCall(MatDestroy(&F));
94: }
95: PetscCall(PetscFree(parts));
96: PetscCall(PetscRandomDestroy(&rctx));
98: PetscCall(PetscFinalize());
99: return 0;
100: }
102: PetscErrorCode MatIsDiagonal(Mat A)
103: {
104: PetscInt ncols, i, j, rstart, rend;
105: const PetscInt *cols;
106: const PetscScalar *vals;
107: PetscBool founddiag;
109: PetscFunctionBeginUser;
110: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
111: for (i = rstart; i < rend; i++) {
112: founddiag = PETSC_FALSE;
113: PetscCall(MatGetRow(A, i, &ncols, &cols, &vals));
114: for (j = 0; j < ncols; j++) {
115: if (cols[j] == i) {
116: PetscCheck(PetscAbsScalar(vals[j] - 1) < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Row %" PetscInt_FMT " does not have 1 on the diagonal, it has %g", i, (double)PetscAbsScalar(vals[j]));
117: founddiag = PETSC_TRUE;
118: } else {
119: PetscCheck(PetscAbsScalar(vals[j]) < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Row %" PetscInt_FMT " has off-diagonal value %g at %" PetscInt_FMT, i, (double)PetscAbsScalar(vals[j]), cols[j]);
120: }
121: }
122: PetscCall(MatRestoreRow(A, i, &ncols, &cols, &vals));
123: PetscCheck(founddiag, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Row %" PetscInt_FMT " does not have diagonal entry", i);
124: }
125: PetscFunctionReturn(PETSC_SUCCESS);
126: }
128: /*
129: All processes receive all the block information
130: */
131: PetscErrorCode BuildMatrix(const PetscInt *parts, PetscInt nblocks, const PetscInt *blocksizes, Mat *A)
132: {
133: PetscInt i, cnt = 0;
134: PetscMPIInt rank;
136: PetscFunctionBeginUser;
137: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
138: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, parts[rank], parts[rank], PETSC_DETERMINE, PETSC_DETERMINE, 0, NULL, 0, NULL, A));
139: PetscCall(MatSetOption(*A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
140: if (rank == 0) {
141: for (i = 0; i < nblocks; i++) {
142: PetscCall(MatSetValue(*A, cnt, cnt + blocksizes[i] - 1, 1.0, INSERT_VALUES));
143: PetscCall(MatSetValue(*A, cnt + blocksizes[i] - 1, cnt, 1.0, INSERT_VALUES));
144: cnt += blocksizes[i];
145: }
146: }
147: PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
148: PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
149: PetscCall(MatShift(*A, 10));
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: /*TEST
155: test:
157: test:
158: suffix: 2
159: nsize: 2
161: test:
162: suffix: 5
163: nsize: 5
165: TEST*/