1 2 static char help[] = "Tests MatInvertVariableBlockEnvelope()\n\n"; 3 4 #include <petscmat.h> 5 extern PetscErrorCode MatIsDiagonal(Mat); 6 extern PetscErrorCode BuildMatrix(const PetscInt *, PetscInt, const PetscInt *, Mat *); 7 8 int main(int argc, char **argv) { 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, (char *)0, 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 102 PetscErrorCode MatIsDiagonal(Mat A) { 103 PetscInt ncols, i, j, rstart, rend; 104 const PetscInt *cols; 105 const PetscScalar *vals; 106 PetscBool founddiag; 107 108 PetscFunctionBeginUser; 109 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 110 for (i = rstart; i < rend; i++) { 111 founddiag = PETSC_FALSE; 112 PetscCall(MatGetRow(A, i, &ncols, &cols, &vals)); 113 for (j = 0; j < ncols; j++) { 114 if (cols[j] == i) { 115 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])); 116 founddiag = PETSC_TRUE; 117 } else { 118 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]); 119 } 120 } 121 PetscCall(MatRestoreRow(A, i, &ncols, &cols, &vals)); 122 PetscCheck(founddiag, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Row %" PetscInt_FMT " does not have diagonal entrie", i); 123 } 124 PetscFunctionReturn(0); 125 } 126 127 /* 128 All processes receive all the block information 129 */ 130 PetscErrorCode BuildMatrix(const PetscInt *parts, PetscInt nblocks, const PetscInt *blocksizes, Mat *A) { 131 PetscInt i, cnt = 0; 132 PetscMPIInt rank; 133 134 PetscFunctionBeginUser; 135 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 136 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, parts[rank], parts[rank], PETSC_DETERMINE, PETSC_DETERMINE, 0, NULL, 0, NULL, A)); 137 PetscCall(MatSetOption(*A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 138 if (rank == 0) { 139 for (i = 0; i < nblocks; i++) { 140 PetscCall(MatSetValue(*A, cnt, cnt + blocksizes[i] - 1, 1.0, INSERT_VALUES)); 141 PetscCall(MatSetValue(*A, cnt + blocksizes[i] - 1, cnt, 1.0, INSERT_VALUES)); 142 cnt += blocksizes[i]; 143 } 144 } 145 PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY)); 146 PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY)); 147 PetscCall(MatShift(*A, 10)); 148 PetscFunctionReturn(0); 149 } 150 151 /*TEST 152 153 test: 154 155 test: 156 suffix: 2 157 nsize: 2 158 159 test: 160 suffix: 5 161 nsize: 5 162 163 TEST*/ 164