1 static char help[] = "Tests MatInvertVariableBlockEnvelope()\n\n";
2
3 #include <petscmat.h>
4 extern PetscErrorCode MatIsDiagonal(Mat);
5 extern PetscErrorCode BuildMatrix(const PetscInt *, PetscInt, const PetscInt *, Mat *);
6
main(int argc,char ** argv)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;
15
16 PetscFunctionBeginUser;
17 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
18 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_DENSE));
19
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));
36
37 PetscCall(MatMatTransposeMult(C, C, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &A));
38 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
39
40 PetscCall(MatInvertVariableBlockEnvelope(A, MAT_INITIAL_MATRIX, &D));
41 PetscCall(MatView(D, PETSC_VIEWER_STDOUT_WORLD));
42
43 PetscCall(MatMatMult(A, D, MAT_INITIAL_MATRIX, 1.0, &F));
44 PetscCall(MatView(F, PETSC_VIEWER_STDOUT_WORLD));
45 PetscCall(MatIsDiagonal(F));
46
47 PetscCall(MatDestroy(&A));
48 PetscCall(MatDestroy(&D));
49 PetscCall(MatDestroy(&C));
50 PetscCall(MatDestroy(&F));
51
52 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rctx));
53 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
54 PetscCall(PetscMalloc1(size, &parts));
55
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;
65
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;
82
83 PetscCall(BuildMatrix(parts, nblocks, blocksizes, &A));
84 PetscCall(PetscFree(blocksizes));
85
86 PetscCall(MatInvertVariableBlockEnvelope(A, MAT_INITIAL_MATRIX, &D));
87
88 PetscCall(MatMatMult(A, D, MAT_INITIAL_MATRIX, 1.0, &F));
89 PetscCall(MatIsDiagonal(F));
90
91 PetscCall(MatDestroy(&A));
92 PetscCall(MatDestroy(&D));
93 PetscCall(MatDestroy(&F));
94 }
95 PetscCall(PetscFree(parts));
96 PetscCall(PetscRandomDestroy(&rctx));
97
98 PetscCall(PetscFinalize());
99 return 0;
100 }
101
MatIsDiagonal(Mat A)102 PetscErrorCode MatIsDiagonal(Mat A)
103 {
104 PetscInt ncols, i, j, rstart, rend;
105 const PetscInt *cols;
106 const PetscScalar *vals;
107 PetscBool founddiag;
108
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 }
127
128 /*
129 All processes receive all the block information
130 */
BuildMatrix(const PetscInt * parts,PetscInt nblocks,const PetscInt * blocksizes,Mat * A)131 PetscErrorCode BuildMatrix(const PetscInt *parts, PetscInt nblocks, const PetscInt *blocksizes, Mat *A)
132 {
133 PetscInt i, cnt = 0;
134 PetscMPIInt rank;
135
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 }
152
153 /*TEST
154
155 test:
156
157 test:
158 suffix: 2
159 nsize: 2
160
161 test:
162 suffix: 5
163 nsize: 5
164
165 TEST*/
166