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 { 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; 16 17 PetscFunctionBeginUser; 18 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 19 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_DENSE)); 20 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)); 37 38 PetscCall(MatMatTransposeMult(C, C, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &A)); 39 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 40 41 PetscCall(MatInvertVariableBlockEnvelope(A, MAT_INITIAL_MATRIX, &D)); 42 PetscCall(MatView(D, PETSC_VIEWER_STDOUT_WORLD)); 43 44 PetscCall(MatMatMult(A, D, MAT_INITIAL_MATRIX, 1.0, &F)); 45 PetscCall(MatView(F, PETSC_VIEWER_STDOUT_WORLD)); 46 PetscCall(MatIsDiagonal(F)); 47 48 PetscCall(MatDestroy(&A)); 49 PetscCall(MatDestroy(&D)); 50 PetscCall(MatDestroy(&C)); 51 PetscCall(MatDestroy(&F)); 52 53 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rctx)); 54 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 55 PetscCall(PetscMalloc1(size, &parts)); 56 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; 66 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; 83 84 PetscCall(BuildMatrix(parts, nblocks, blocksizes, &A)); 85 PetscCall(PetscFree(blocksizes)); 86 87 PetscCall(MatInvertVariableBlockEnvelope(A, MAT_INITIAL_MATRIX, &D)); 88 89 PetscCall(MatMatMult(A, D, MAT_INITIAL_MATRIX, 1.0, &F)); 90 PetscCall(MatIsDiagonal(F)); 91 92 PetscCall(MatDestroy(&A)); 93 PetscCall(MatDestroy(&D)); 94 PetscCall(MatDestroy(&F)); 95 } 96 PetscCall(PetscFree(parts)); 97 PetscCall(PetscRandomDestroy(&rctx)); 98 99 PetscCall(PetscFinalize()); 100 return 0; 101 } 102 103 PetscErrorCode MatIsDiagonal(Mat A) 104 { 105 PetscInt ncols, i, j, rstart, rend; 106 const PetscInt *cols; 107 const PetscScalar *vals; 108 PetscBool founddiag; 109 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 entrie", i); 125 } 126 PetscFunctionReturn(0); 127 } 128 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; 136 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(0); 152 } 153 154 /*TEST 155 156 test: 157 158 test: 159 suffix: 2 160 nsize: 2 161 162 test: 163 suffix: 5 164 nsize: 5 165 166 TEST*/ 167