Actual source code: ex108.c
1: static char help[] = "Testing MatCreateSeqBAIJWithArrays() and MatCreateSeqSBAIJWithArrays().\n\n";
3: #include <petscmat.h>
5: int main(int argc, char **argv)
6: {
7: Mat A, B, As;
8: const PetscInt *ai, *aj;
9: PetscInt i, j, k, nz, n, asi[] = {0, 2, 3, 4, 6, 7};
10: PetscInt asj[] = {0, 4, 1, 2, 3, 4, 4};
11: PetscScalar asa[7], *aa;
12: PetscRandom rctx;
13: PetscMPIInt size;
14: PetscBool flg;
17: PetscInitialize(&argc, &argv, (char *)0, help);
18: MPI_Comm_size(PETSC_COMM_WORLD, &size);
21: /* Create a aij matrix for checking */
22: MatCreateSeqAIJ(PETSC_COMM_SELF, 5, 5, 2, NULL, &A);
23: PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
24: PetscRandomSetFromOptions(rctx);
26: k = 0;
27: for (i = 0; i < 5; i++) {
28: nz = asi[i + 1] - asi[i]; /* length of i_th row of A */
29: for (j = 0; j < nz; j++) {
30: PetscRandomGetValue(rctx, &asa[k]);
31: MatSetValues(A, 1, &i, 1, &asj[k], &asa[k], INSERT_VALUES);
32: MatSetValues(A, 1, &i, 1, &asj[k], &asa[k], INSERT_VALUES);
33: if (i != asj[k]) { /* insert symmetric entry */
34: MatSetValues(A, 1, &asj[k], 1, &i, &asa[k], INSERT_VALUES);
35: }
36: k++;
37: }
38: }
39: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
40: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
42: /* Create a baij matrix using MatCreateSeqBAIJWithArrays() */
43: MatGetRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &n, &ai, &aj, &flg);
44: MatSeqAIJGetArray(A, &aa);
45: /* WARNING: This sharing is dangerous if either A or B is later assembled */
46: MatCreateSeqBAIJWithArrays(PETSC_COMM_SELF, 1, 5, 5, (PetscInt *)ai, (PetscInt *)aj, aa, &B);
47: MatSeqAIJRestoreArray(A, &aa);
48: MatRestoreRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &n, &ai, &aj, &flg);
49: MatMultEqual(A, B, 10, &flg);
52: /* Create a sbaij matrix using MatCreateSeqSBAIJWithArrays() */
53: MatCreateSeqSBAIJWithArrays(PETSC_COMM_SELF, 1, 5, 5, asi, asj, asa, &As);
54: MatMultEqual(A, As, 10, &flg);
57: /* Free spaces */
58: PetscRandomDestroy(&rctx);
59: MatDestroy(&A);
60: MatDestroy(&B);
61: MatDestroy(&As);
62: PetscFinalize();
63: return 0;
64: }