Actual source code: ex108.c
petsc-3.14.6 2021-03-30
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;
14: PetscMPIInt size;
15: PetscBool flg;
17: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
18: MPI_Comm_size(PETSC_COMM_WORLD,&size);
19: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
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);
50: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"MatMult(A,B) are NOT equal");
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);
55: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"MatMult(A,As) are NOT equal");
57: /* Free spaces */
58: PetscRandomDestroy(&rctx);
59: MatDestroy(&A);
60: MatDestroy(&B);
61: MatDestroy(&As);
62: PetscFinalize();
63: return ierr;
64: }