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;
18: PetscInitialize(&argc, &argv, (char *)0, help);
19: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_DENSE);
21: MatCreate(PETSC_COMM_WORLD, &C);
22: MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, 6, 18);
23: MatSetFromOptions(C);
24: 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: MatSetValues(C, 2, rows, 2, rows, (PetscScalar *)values, INSERT_VALUES);
33: }
34: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
35: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
36: MatView(C, PETSC_VIEWER_STDOUT_WORLD);
38: MatMatTransposeMult(C, C, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &A);
39: MatView(A, PETSC_VIEWER_STDOUT_WORLD);
41: MatInvertVariableBlockEnvelope(A, MAT_INITIAL_MATRIX, &D);
42: MatView(D, PETSC_VIEWER_STDOUT_WORLD);
44: MatMatMult(A, D, MAT_INITIAL_MATRIX, 1.0, &F);
45: MatView(F, PETSC_VIEWER_STDOUT_WORLD);
46: MatIsDiagonal(F);
48: MatDestroy(&A);
49: MatDestroy(&D);
50: MatDestroy(&C);
51: MatDestroy(&F);
53: PetscRandomCreate(PETSC_COMM_SELF, &rctx);
54: MPI_Comm_size(PETSC_COMM_WORLD, &size);
55: PetscMalloc1(size, &parts);
57: for (j = 0; j < 128; j++) {
58: cnt = 0;
59: for (i = 0; i < size - 1; i++) {
60: 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: PetscRandomGetValueReal(rctx, &rand);
68: nblocks = rand * 10;
69: nblocks = PetscMax(nblocks, 2);
70: cnt = 0;
71: PetscMalloc1(nblocks, &blocksizes);
72: for (i = 0; i < nblocks - 1; i++) {
73: 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: BuildMatrix(parts, nblocks, blocksizes, &A);
85: PetscFree(blocksizes);
87: MatInvertVariableBlockEnvelope(A, MAT_INITIAL_MATRIX, &D);
89: MatMatMult(A, D, MAT_INITIAL_MATRIX, 1.0, &F);
90: MatIsDiagonal(F);
92: MatDestroy(&A);
93: MatDestroy(&D);
94: MatDestroy(&F);
95: }
96: PetscFree(parts);
97: PetscRandomDestroy(&rctx);
99: 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;
111: MatGetOwnershipRange(A, &rstart, &rend);
112: for (i = rstart; i < rend; i++) {
113: founddiag = PETSC_FALSE;
114: MatGetRow(A, i, &ncols, &cols, &vals);
115: for (j = 0; j < ncols; j++) {
116: if (cols[j] == i) {
118: founddiag = PETSC_TRUE;
119: } else {
121: }
122: }
123: MatRestoreRow(A, i, &ncols, &cols, &vals);
125: }
126: return 0;
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;
138: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
139: MatCreateAIJ(PETSC_COMM_WORLD, parts[rank], parts[rank], PETSC_DETERMINE, PETSC_DETERMINE, 0, NULL, 0, NULL, A);
140: MatSetOption(*A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
141: if (rank == 0) {
142: for (i = 0; i < nblocks; i++) {
143: MatSetValue(*A, cnt, cnt + blocksizes[i] - 1, 1.0, INSERT_VALUES);
144: MatSetValue(*A, cnt + blocksizes[i] - 1, cnt, 1.0, INSERT_VALUES);
145: cnt += blocksizes[i];
146: }
147: }
148: MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY);
149: MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY);
150: MatShift(*A, 10);
151: return 0;
152: }
154: /*TEST
156: test:
158: test:
159: suffix: 2
160: nsize: 2
162: test:
163: suffix: 5
164: nsize: 5
166: TEST*/