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 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, (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 { 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 */ 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