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