1 const char help[] = "Test MATDIAGONAL"; 2 3 #include <petsc/private/petscimpl.h> 4 #include <petscmat.h> 5 6 int main(int argc, char **argv) 7 { 8 Vec a, a2, b, b2, c, c2, A_diag, A_inv_diag; 9 Mat A, B; 10 PetscInt n = 10; 11 12 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 13 14 MPI_Comm comm = PETSC_COMM_SELF; 15 PetscCall(VecCreateSeq(comm, n, &a)); 16 PetscCall(VecDuplicate(a, &b)); 17 PetscCall(VecDuplicate(a, &c)); 18 PetscRandom rand; 19 20 PetscCall(PetscRandomCreate(comm, &rand)); 21 PetscCall(VecSetRandom(a, rand)); 22 PetscCall(VecSetRandom(b, rand)); 23 24 PetscCall(VecDuplicate(a, &a2)); 25 PetscCall(VecCopy(a, a2)); 26 PetscCall(VecDuplicate(b, &b2)); 27 PetscCall(VecCopy(b, b2)); 28 PetscCall(VecDuplicate(c, &c2)); 29 30 PetscCall(MatCreateDiagonal(a2, &A)); 31 PetscCall(MatCreateDiagonal(b2, &B)); 32 PetscCall(VecDestroy(&a2)); 33 PetscCall(VecDestroy(&b2)); 34 35 PetscCall(VecDuplicate(a, &a2)); 36 PetscCall(VecDuplicate(b, &b2)); 37 38 PetscCall(MatAXPY(A, 0.5, B, SAME_NONZERO_PATTERN)); 39 PetscCall(VecAXPY(a, 0.5, b)); 40 41 PetscReal mat_norm, vec_norm; 42 PetscCall(VecNorm(a, NORM_2, &vec_norm)); 43 PetscCall(MatNorm(A, NORM_FROBENIUS, &mat_norm)); 44 PetscCheck(vec_norm == mat_norm, comm, PETSC_ERR_PLIB, "Norms don't match"); 45 46 // For diagonal matrix, all operator norms are the max norm of the vector 47 PetscCall(VecNorm(a, NORM_INFINITY, &vec_norm)); 48 PetscCall(MatNorm(A, NORM_INFINITY, &mat_norm)); 49 PetscCheck(vec_norm == mat_norm, comm, PETSC_ERR_PLIB, "Norms don't match"); 50 PetscCall(MatNorm(A, NORM_1, &mat_norm)); 51 PetscCheck(vec_norm == mat_norm, comm, PETSC_ERR_PLIB, "Norms don't match"); 52 53 PetscCall(VecPointwiseMult(c, b, a)); 54 PetscCall(MatMult(A, b, c2)); 55 PetscCall(VecAXPY(c2, -1.0, c)); 56 PetscCall(VecNorm(c2, NORM_INFINITY, &vec_norm)); 57 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatMult is not like VecPointwiseMultiply"); 58 59 PetscCall(VecPointwiseMult(c, b, a)); 60 PetscCall(VecAXPY(c, 1.0, a)); 61 PetscCall(MatMultAdd(A, b, a, c2)); 62 PetscCall(VecAXPY(c2, -1.0, c)); 63 PetscCall(VecNorm(c2, NORM_INFINITY, &vec_norm)); 64 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatMultAdd gave unexpected value"); 65 66 PetscCall(VecSet(c, 1.2)); 67 PetscCall(VecSet(c2, 1.2)); 68 PetscCall(VecPointwiseMult(c, b, a)); 69 PetscCall(VecAXPY(c, 1.0, c2)); 70 PetscCall(MatMultAdd(A, b, c2, c2)); 71 PetscCall(VecAXPY(c2, -1.0, c)); 72 PetscCall(VecNorm(c2, NORM_INFINITY, &vec_norm)); 73 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatMultAdd gave unexpected value"); 74 75 PetscCall(VecPointwiseDivide(c, b, a)); 76 PetscCall(MatSolve(A, b, c2)); 77 PetscCall(VecAXPY(c2, -1.0, c)); 78 PetscCall(VecNorm(c2, NORM_INFINITY, &vec_norm)); 79 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatMult is not like VecPointwiseMultiply"); 80 81 Mat A_dup; 82 PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &A_dup)); 83 PetscCall(MatDestroy(&A_dup)); 84 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A_dup)); 85 PetscCall(MatGetDiagonal(A_dup, a2)); 86 PetscCall(VecAXPY(a2, -1.0, a)); 87 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm)); 88 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDuplicate with MAT_COPY_VALUES did not make a duplicate vector"); 89 PetscCall(MatDestroy(&A_dup)); 90 91 PetscCall(MatShift(A, 1.5)); 92 PetscCall(VecShift(a, 1.5)); 93 PetscCall(MatGetDiagonal(A, a2)); 94 PetscCall(VecAXPY(a2, -1.0, a)); 95 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm)); 96 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatShift gave different result from VecShift"); 97 98 PetscCall(MatScale(A, 0.75)); 99 PetscCall(VecScale(a, 0.75)); 100 PetscCall(MatGetDiagonal(A, a2)); 101 PetscCall(VecAXPY(a2, -1.0, a)); 102 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm)); 103 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatScale gave different result from VecScale"); 104 105 PetscCall(VecPointwiseMult(a, a, b)); 106 PetscCall(MatDiagonalScale(A, b, NULL)); 107 PetscCall(MatGetDiagonal(A, a2)); 108 PetscCall(VecAXPY(a2, -1.0, a)); 109 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm)); 110 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDiagonalScale gave unexpected result"); 111 112 PetscCall(VecPointwiseMult(a, a, b)); 113 PetscCall(MatDiagonalScale(A, NULL, b)); 114 PetscCall(MatGetDiagonal(A, a2)); 115 PetscCall(VecAXPY(a2, -1.0, a)); 116 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm)); 117 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDiagonalScale gave unexpected result"); 118 119 PetscCall(VecCopy(b, a)); 120 PetscCall(MatDiagonalSet(A, b, INSERT_VALUES)); 121 PetscCall(MatGetDiagonal(A, a2)); 122 PetscCall(VecAXPY(a2, -1.0, a)); 123 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm)); 124 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDiagonalSet gave unexpected result"); 125 126 PetscCall(VecSetRandom(a, rand)); 127 PetscCall(VecSetRandom(b, rand)); 128 PetscCall(MatDiagonalSet(A, a, INSERT_VALUES)); 129 PetscCall(VecAXPY(a, 1.0, b)); 130 PetscCall(MatDiagonalSet(A, b, ADD_VALUES)); 131 PetscCall(MatGetDiagonal(A, a2)); 132 PetscCall(VecAXPY(a2, -1.0, a)); 133 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm)); 134 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDiagonalSet gave unexpected result"); 135 136 PetscCall(VecSetRandom(a, rand)); 137 PetscCall(VecSetRandom(b, rand)); 138 PetscCall(MatDiagonalSet(A, a, INSERT_VALUES)); 139 PetscCall(VecPointwiseMax(a, a, b)); 140 PetscCall(MatDiagonalSet(A, b, MAX_VALUES)); 141 PetscCall(MatGetDiagonal(A, a2)); 142 PetscCall(VecAXPY(a2, -1.0, a)); 143 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm)); 144 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDiagonalSet gave unexpected result"); 145 146 PetscCall(VecSetRandom(a, rand)); 147 PetscCall(VecSetRandom(b, rand)); 148 PetscCall(MatDiagonalSet(A, a, INSERT_VALUES)); 149 PetscCall(VecPointwiseMin(a, a, b)); 150 PetscCall(MatDiagonalSet(A, b, MIN_VALUES)); 151 PetscCall(MatGetDiagonal(A, a2)); 152 PetscCall(VecAXPY(a2, -1.0, a)); 153 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm)); 154 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDiagonalSet gave unexpected result"); 155 156 PetscCall(VecSetRandom(a, rand)); 157 PetscCall(VecSetRandom(b, rand)); 158 PetscCall(MatDiagonalSet(A, a, INSERT_VALUES)); 159 PetscCall(MatDiagonalSet(A, b, NOT_SET_VALUES)); 160 PetscCall(MatGetDiagonal(A, a2)); 161 PetscCall(VecAXPY(a2, -1.0, a)); 162 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm)); 163 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDiagonalSet gave unexpected result"); 164 165 PetscCall(VecSet(a2, 0.5)); 166 167 PetscObjectState state_pre, state_post; 168 PetscCall(PetscObjectStateGet((PetscObject)A, &state_pre)); 169 PetscCall(MatDiagonalGetInverseDiagonal(A, &A_inv_diag)); 170 PetscCall(MatDiagonalRestoreInverseDiagonal(A, &A_inv_diag)); 171 PetscCall(PetscObjectStateGet((PetscObject)A, &state_post)); 172 PetscCheck(state_pre == state_post, comm, PETSC_ERR_PLIB, "State changed on noop"); 173 174 PetscCall(PetscObjectStateGet((PetscObject)A, &state_pre)); 175 PetscCall(MatDiagonalGetInverseDiagonal(A, &A_inv_diag)); 176 PetscCall(VecSet(A_inv_diag, 2.0)); 177 PetscCall(MatDiagonalRestoreInverseDiagonal(A, &A_inv_diag)); 178 PetscCall(PetscObjectStateGet((PetscObject)A, &state_post)); 179 PetscCheck(state_pre != state_post, comm, PETSC_ERR_PLIB, "State not changed on mutation"); 180 181 PetscCall(PetscObjectStateGet((PetscObject)A, &state_pre)); 182 PetscCall(MatDiagonalGetDiagonal(A, &A_diag)); 183 PetscCall(MatDiagonalRestoreDiagonal(A, &A_diag)); 184 PetscCall(PetscObjectStateGet((PetscObject)A, &state_post)); 185 PetscCheck(state_pre == state_post, comm, PETSC_ERR_PLIB, "State changed on noop"); 186 187 PetscCall(MatDiagonalGetDiagonal(A, &A_diag)); 188 PetscCall(VecAXPY(a2, -1.0, A_diag)); 189 PetscCall(VecSet(A_diag, 1.0)); 190 PetscCall(MatDiagonalRestoreDiagonal(A, &A_diag)); 191 PetscCall(PetscObjectStateGet((PetscObject)A, &state_post)); 192 PetscCheck(state_pre != state_post, comm, PETSC_ERR_PLIB, "State not changed on mutation"); 193 194 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm)); 195 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDiagonalGetInverse gave unexpected result"); 196 197 PetscCall(MatZeroEntries(A)); 198 PetscCall(MatNorm(A, NORM_INFINITY, &mat_norm)); 199 PetscCheck(mat_norm == 0.0, comm, PETSC_ERR_PLIB, "MatZeroEntries gave unexpected result"); 200 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF)); 201 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF, PETSC_VIEWER_ASCII_INFO)); 202 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF)); 203 PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF)); 204 205 PetscCall(MatDestroy(&A)); 206 PetscCall(MatDestroy(&B)); 207 208 PetscCall(PetscRandomDestroy(&rand)); 209 PetscCall(VecDestroy(&c2)); 210 PetscCall(VecDestroy(&b2)); 211 PetscCall(VecDestroy(&a2)); 212 PetscCall(VecDestroy(&c)); 213 PetscCall(VecDestroy(&b)); 214 PetscCall(VecDestroy(&a)); 215 PetscCall(PetscFinalize()); 216 return 0; 217 } 218 219 /*TEST 220 221 test: 222 suffix: 0 223 224 TEST*/ 225