Actual source code: ex55.c
2: static char help[] = "Tests converting a matrix to another format with MatConvert().\n\n";
4: #include petscmat.h
5: /* Usage: mpirun -np <np> ex55 -display <0 or 1> */
9: int main(int argc,char **args)
10: {
11: Mat C,A,B,D;
13: PetscInt i,j,ntypes,bs,mbs,m,block,d_nz=6, o_nz=3,col[3],row,displ=0;
14: PetscMPIInt size,rank;
15: /* const MatType type[9] = {MATMPIAIJ,MATMPIBAIJ,MATMPIROWBS};*/ /* BlockSolve95 is required for MATMPIROWBS */
16: MatType type[9];
17: char file[PETSC_MAX_PATH_LEN];
18: PetscViewer fd;
19: PetscTruth equal,flg_loadmat,flg;
20: PetscScalar value[3];
22: PetscInitialize(&argc,&args,(char *)0,help);
23: PetscOptionsGetInt(PETSC_NULL,"-displ",&displ,PETSC_NULL);
24: PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN-1,&flg_loadmat);
25: MPI_Comm_size(PETSC_COMM_WORLD,&size);
26: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
28: PetscOptionsHasName(PETSC_NULL,"-testseqaij",&flg);
29: if (flg){
30: if (size == 1){
31: type[0] = MATSEQAIJ;
32: } else {
33: type[0] = MATMPIAIJ;
34: }
35: } else {
36: type[0] = MATAIJ;
37: }
38: if (size == 1){
39: ntypes = 3;
40: type[1] = MATSEQBAIJ;
41: type[2] = MATSEQSBAIJ;
42: } else {
43: ntypes = 2;
44: type[1] = MATMPIBAIJ;
45: type[2] = MATMPISBAIJ; /* Matconvert from mpisbaij mat to other formats are not supported */
46: }
48: /* input matrix C */
49: if (flg_loadmat){
50: /* Open binary file. */
51: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
53: /* Load a baij matrix, then destroy the viewer. */
54: if (size == 1){
55: MatLoad(fd,MATSEQBAIJ,&C);
56: } else {
57: MatLoad(fd,MATMPIBAIJ,&C);
58: }
59: PetscViewerDestroy(fd);
60: } else { /* Create a baij mat with bs>1 */
61: bs = 2; mbs=8;
62: PetscOptionsGetInt(PETSC_NULL,"-mbs",&mbs,PETSC_NULL);
63: PetscOptionsGetInt(PETSC_NULL,"-bs",&bs,PETSC_NULL);
64: if (bs <= 1) SETERRQ(PETSC_ERR_ARG_WRONG," bs must be >1 in this case");
65: m = mbs*bs;
66: if (size == 1){
67: MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,m,m,d_nz,PETSC_NULL,&C);
68: } else {
69: MatCreateMPIBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,m,m,d_nz,PETSC_NULL,o_nz,PETSC_NULL,&C);
70: }
71: for (block=0; block<mbs; block++){
72: /* diagonal blocks */
73: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
74: for (i=1+block*bs; i<bs-1+block*bs; i++) {
75: col[0] = i-1; col[1] = i; col[2] = i+1;
76: MatSetValues(C,1,&i,3,col,value,INSERT_VALUES);
77: }
78: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
79: value[0]=-1.0; value[1]=4.0;
80: MatSetValues(C,1,&i,2,col,value,INSERT_VALUES);
82: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
83: value[0]=4.0; value[1] = -1.0;
84: MatSetValues(C,1,&i,2,col,value,INSERT_VALUES);
85: }
86: /* off-diagonal blocks */
87: value[0]=-1.0;
88: for (i=0; i<(mbs-1)*bs; i++){
89: col[0]=i+bs;
90: MatSetValues(C,1,&i,1,col,value,INSERT_VALUES);
91: col[0]=i; row=i+bs;
92: MatSetValues(C,1,&row,1,col,value,INSERT_VALUES);
93: }
94: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
95: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
96: }
97:
98: /* convert C to other formats */
99: for (i=0; i<ntypes; i++) {
100: MatConvert(C,type[i],MAT_INITIAL_MATRIX,&A);
101: MatMultEqual(A,C,10,&equal);
102: if (!equal) SETERRQ1(PETSC_ERR_ARG_NOTSAMETYPE,"Error in conversion from BAIJ to %s",type[i]);
103: for (j=i+1; j<ntypes; j++) {
104: if (displ>0) {
105: PetscPrintf(PETSC_COMM_WORLD," [%d] test conversion between %s and %s\n",rank,type[i],type[j]);
106: }
108: MatConvert(A,type[j],MAT_INITIAL_MATRIX,&B);
109: MatConvert(B,type[i],MAT_INITIAL_MATRIX,&D);
111: if (bs == 1){
112: MatEqual(A,D,&equal);
113: if (!equal){
114: PetscPrintf(PETSC_COMM_SELF," A: %s\n",type[i]);
115: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
116: PetscPrintf(PETSC_COMM_SELF," B: %s\n",type[j]);
117: MatView(B,PETSC_VIEWER_STDOUT_WORLD);
118: PetscPrintf(PETSC_COMM_SELF," D: %s\n",type[i]);
119: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
120: SETERRQ2(1,"Error in conversion from %s to %s",type[i],type[j]);
121: }
122: } else { /* bs > 1 */
123: MatMultEqual(A,B,10,&equal);
124:
125: if (!equal) SETERRQ2(PETSC_ERR_ARG_NOTSAMETYPE,"Error in conversion from %s to %s",type[i],type[j]);
126: }
127: MatDestroy(B);
128: MatDestroy(D);
129: B = PETSC_NULL;
130: D = PETSC_NULL;
131: }
132: /* Test in-place convert */
133: if (size == 1){ /* size > 1 is not working yet! */
134: j = (i+1)%ntypes;
135: /* printf("[%d] i: %d, j: %d\n",rank,i,j); */
136: MatConvert(A,type[j],MAT_REUSE_MATRIX,&A);
137: }
139: MatDestroy(A);
140: }
141: MatDestroy(C);
143: PetscFinalize();
144: return 0;
145: }