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