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