Actual source code: ex235.c
1: static char help[] = "Test combinations of scalings, shifts and get diagonal of MATSHELL\n\n";
3: #include <petscmat.h>
5: static PetscErrorCode myMult(Mat S, Vec x, Vec y)
6: {
7: Mat A;
9: MatShellGetContext(S, &A);
10: MatMult(A, x, y);
11: return 0;
12: }
14: static PetscErrorCode myGetDiagonal(Mat S, Vec d)
15: {
16: Mat A;
18: MatShellGetContext(S, &A);
19: MatGetDiagonal(A, d);
20: return 0;
21: }
23: static PetscErrorCode shiftandscale(Mat A, Vec *D)
24: {
25: Vec ll, d, rr;
27: MatCreateVecs(A, &ll, &rr);
28: MatCreateVecs(A, &d, NULL);
29: VecSetRandom(ll, NULL);
30: VecSetRandom(rr, NULL);
31: VecSetRandom(d, NULL);
32: MatScale(A, 3.0);
33: MatShift(A, -4.0);
34: MatScale(A, 8.0);
35: MatDiagonalSet(A, d, ADD_VALUES);
36: MatShift(A, 9.0);
37: MatScale(A, 8.0);
38: VecSetRandom(ll, NULL);
39: VecSetRandom(rr, NULL);
40: MatDiagonalScale(A, ll, rr);
41: MatShift(A, 2.0);
42: MatScale(A, 11.0);
43: VecSetRandom(d, NULL);
44: MatDiagonalSet(A, d, ADD_VALUES);
45: VecSetRandom(ll, NULL);
46: VecSetRandom(rr, NULL);
47: MatDiagonalScale(A, ll, rr);
48: MatShift(A, 5.0);
49: MatScale(A, 7.0);
50: MatGetDiagonal(A, d);
51: *D = d;
52: VecDestroy(&ll);
53: VecDestroy(&rr);
54: return 0;
55: }
57: int main(int argc, char **args)
58: {
59: Mat A, Aij, B;
60: Vec Adiag, Aijdiag;
61: PetscInt m = 3;
62: PetscReal Aijnorm, Aijdiagnorm, Bnorm, dnorm;
65: PetscInitialize(&argc, &args, (char *)0, help);
66: PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
68: MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, m, m, 7, NULL, 6, NULL, &Aij);
69: MatSetRandom(Aij, NULL);
70: MatSetOption(Aij, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE);
72: MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, m, m, Aij, &A);
73: MatShellSetOperation(A, MATOP_MULT, (void (*)(void))myMult);
74: MatShellSetOperation(A, MATOP_GET_DIAGONAL, (void (*)(void))myGetDiagonal);
76: shiftandscale(A, &Adiag);
77: MatComputeOperator(A, NULL, &B);
78: shiftandscale(Aij, &Aijdiag);
79: MatAXPY(Aij, -1.0, B, DIFFERENT_NONZERO_PATTERN);
80: MatNorm(Aij, NORM_FROBENIUS, &Aijnorm);
81: MatNorm(B, NORM_FROBENIUS, &Bnorm);
83: VecAXPY(Aijdiag, -1.0, Adiag);
84: VecNorm(Adiag, NORM_2, &dnorm);
85: VecNorm(Aijdiag, NORM_2, &Aijdiagnorm);
87: MatDestroy(&A);
88: MatDestroy(&Aij);
89: VecDestroy(&Adiag);
90: VecDestroy(&Aijdiag);
91: MatDestroy(&B);
92: PetscFinalize();
93: return 0;
94: }
96: /*TEST
98: test:
99: nsize: {{1 2 3 4}}
101: TEST*/