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: mpiexec -n <np> ex55 -verbose <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,verbose=0;
14: PetscMPIInt size,rank;
15: const MatType type[9];
16: char file[PETSC_MAX_PATH_LEN];
17: PetscViewer fd;
18: PetscBool equal,flg_loadmat,flg;
19: PetscScalar value[3];
21: PetscInitialize(&argc,&args,(char *)0,help);
22: PetscOptionsGetInt(PETSC_NULL,"-verbose",&verbose,PETSC_NULL);
23: PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg_loadmat);
24: MPI_Comm_size(PETSC_COMM_WORLD,&size);
25: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
27: PetscOptionsHasName(PETSC_NULL,"-testseqaij",&flg);
28: if (flg){
29: if (size == 1){
30: type[0] = MATSEQAIJ;
31: } else {
32: type[0] = MATMPIAIJ;
33: }
34: } else {
35: type[0] = MATAIJ;
36: }
37: if (size == 1){
38: ntypes = 3;
39: type[1] = MATSEQBAIJ;
40: type[2] = MATSEQSBAIJ;
41: } else {
42: ntypes = 3;
43: type[1] = MATMPIBAIJ;
44: type[2] = MATMPISBAIJ;
45: }
47: /* input matrix C */
48: if (flg_loadmat){
49: /* Open binary file. */
50: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
52: /* Load a baij matrix, then destroy the viewer. */
53: MatCreate(PETSC_COMM_WORLD,&C);
54: if (size == 1){
55: MatSetType(C,MATSEQBAIJ);
56: MatLoad(C,fd);
57: } else {
58: MatSetType(C,MATMPIBAIJ);
59: MatLoad(C,fd);
60: }
61: PetscViewerDestroy(&fd);
62: } else { /* Create a baij mat with bs>1 */
63: bs = 2; mbs=8;
64: PetscOptionsGetInt(PETSC_NULL,"-mbs",&mbs,PETSC_NULL);
65: PetscOptionsGetInt(PETSC_NULL,"-bs",&bs,PETSC_NULL);
66: if (bs <= 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG," bs must be >1 in this case");
67: m = mbs*bs;
68: if (size == 1){
69: MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,m,m,d_nz,PETSC_NULL,&C);
70: } else {
71: MatCreateMPIBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,m,m,d_nz,PETSC_NULL,o_nz,PETSC_NULL,&C);
72: }
73: for (block=0; block<mbs; block++){
74: /* diagonal blocks */
75: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
76: for (i=1+block*bs; i<bs-1+block*bs; i++) {
77: col[0] = i-1; col[1] = i; col[2] = i+1;
78: MatSetValues(C,1,&i,3,col,value,INSERT_VALUES);
79: }
80: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
81: value[0]=-1.0; value[1]=4.0;
82: MatSetValues(C,1,&i,2,col,value,INSERT_VALUES);
84: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
85: value[0]=4.0; value[1] = -1.0;
86: MatSetValues(C,1,&i,2,col,value,INSERT_VALUES);
87: }
88: /* off-diagonal blocks */
89: value[0]=-1.0;
90: for (i=0; i<(mbs-1)*bs; i++){
91: col[0]=i+bs;
92: MatSetValues(C,1,&i,1,col,value,INSERT_VALUES);
93: col[0]=i; row=i+bs;
94: MatSetValues(C,1,&row,1,col,value,INSERT_VALUES);
95: }
96: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
97: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
98: }
100: {
101: /* Check the symmetry of C because it will be converted to a sbaij matrix */
102: Mat Ctrans;
103: MatTranspose(C, MAT_INITIAL_MATRIX,&Ctrans);
104: MatEqual(C, Ctrans, &flg);
105: if (flg) {
106: MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE);
107: } else {
108: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"C must be symmetric for this example");
109: }
110: MatDestroy(&Ctrans);
111: }
112: //MatView(C,PETSC_VIEWER_STDOUT_WORLD);
113:
114: /* convert C to other formats */
115: for (i=0; i<ntypes; i++) {
116: MatConvert(C,type[i],MAT_INITIAL_MATRIX,&A);
117: MatMultEqual(A,C,10,&equal);
118: if (!equal) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in conversion from BAIJ to %s",type[i]);
119: for (j=i+1; j<ntypes; j++) {
120: if (verbose>0) {
121: PetscPrintf(PETSC_COMM_WORLD," \n[%d] test conversion between %s and %s\n",rank,type[i],type[j]);
122: }
123:
124: if (!rank && verbose) printf("Convert %s A to %s B\n",type[i],type[j]);
125: MatConvert(A,type[j],MAT_INITIAL_MATRIX,&B);
126: /*
127: if (j == 2){
128: PetscPrintf(PETSC_COMM_SELF," A: %s\n",type[i]);
129: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
130: PetscPrintf(PETSC_COMM_SELF," B: %s\n",type[j]);
131: MatView(B,PETSC_VIEWER_STDOUT_WORLD);
132: }
133: */
134: MatMultEqual(A,B,10,&equal);
135: if (!equal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in conversion from %s to %s",type[i],type[j]);
137: if (size == 1 || j != 2){ /* Matconvert from mpisbaij mat to other formats are not supported */
138: if (!rank && verbose) printf("Convert %s B to %s D\n",type[j],type[i]);
139: MatConvert(B,type[i],MAT_INITIAL_MATRIX,&D);
140: MatMultEqual(B,D,10,&equal);
141: if (!equal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in conversion from %s to %s",type[j],type[i]);
143: MatDestroy(&D);
144: }
145: MatDestroy(&B);
146: MatDestroy(&D);
147: }
149: /* Test in-place convert */
150: if (size == 1){ /* size > 1 is not working yet! */
151: j = (i+1)%ntypes;
152: /* printf("[%d] i: %d, j: %d\n",rank,i,j); */
153: MatConvert(A,type[j],MAT_REUSE_MATRIX,&A);
154: }
156: MatDestroy(&A);
157: }
158: MatDestroy(&C);
160: PetscFinalize();
161: return 0;
162: }