Actual source code: ex177.c
2: static char help[] = "Tests various routines in MatKAIJ format.\n";
4: #include <petscmat.h>
5: #define IMAX 15
7: int main(int argc,char **args)
8: {
9: Mat A,B,TA;
10: PetscScalar *S,*T;
11: PetscViewer fd;
12: char file[PETSC_MAX_PATH_LEN];
13: PetscInt m,n,M,N,p=1,q=1,i,j;
14: PetscMPIInt rank,size;
15: PetscBool flg;
17: PetscInitialize(&argc,&args,(char*)0,help);
18: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
19: MPI_Comm_size(PETSC_COMM_WORLD,&size);
21: /* Load AIJ matrix A */
22: PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);
24: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
25: MatCreate(PETSC_COMM_WORLD,&A);
26: MatLoad(A,fd);
27: PetscViewerDestroy(&fd);
29: /* Get dof, then create S and T */
30: PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL);
31: PetscOptionsGetInt(NULL,NULL,"-q",&q,NULL);
32: PetscMalloc2(p*q,&S,p*q,&T);
33: for (i=0; i<p*q; i++) S[i] = 0;
35: for (i=0; i<p; i++) {
36: for (j=0; j<q; j++) {
37: /* Set some random non-zero values */
38: S[i+p*j] = ((PetscReal) ((i+1)*(j+1))) / ((PetscReal) (p+q));
39: T[i+p*j] = ((PetscReal) ((p-i)+j)) / ((PetscReal) (p*q));
40: }
41: }
43: /* Test KAIJ when both S & T are not NULL */
45: /* Create KAIJ matrix TA */
46: MatCreateKAIJ(A,p,q,S,T,&TA);
47: MatGetLocalSize(A,&m,&n);
48: MatGetSize(A,&M,&N);
50: MatConvert(TA,MATAIJ,MAT_INITIAL_MATRIX,&B);
52: /* Test MatKAIJGetScaledIdentity() */
53: MatKAIJGetScaledIdentity(TA,&flg);
55: /* Test MatMult() */
56: MatMultEqual(TA,B,10,&flg);
58: /* Test MatMultAdd() */
59: MatMultAddEqual(TA,B,10,&flg);
62: MatDestroy(&TA);
63: MatDestroy(&B);
65: /* Test KAIJ when S is NULL */
67: /* Create KAIJ matrix TA */
68: MatCreateKAIJ(A,p,q,NULL,T,&TA);
69: MatGetLocalSize(A,&m,&n);
70: MatGetSize(A,&M,&N);
72: MatConvert(TA,MATAIJ,MAT_INITIAL_MATRIX,&B);
74: /* Test MatKAIJGetScaledIdentity() */
75: MatKAIJGetScaledIdentity(TA,&flg);
77: /* Test MatMult() */
78: MatMultEqual(TA,B,10,&flg);
80: /* Test MatMultAdd() */
81: MatMultAddEqual(TA,B,10,&flg);
84: MatDestroy(&TA);
85: MatDestroy(&B);
87: /* Test KAIJ when T is NULL */
89: /* Create KAIJ matrix TA */
90: MatCreateKAIJ(A,p,q,S,NULL,&TA);
91: MatGetLocalSize(A,&m,&n);
92: MatGetSize(A,&M,&N);
94: MatConvert(TA,MATAIJ,MAT_INITIAL_MATRIX,&B);
96: /* Test MatKAIJGetScaledIdentity() */
97: MatKAIJGetScaledIdentity(TA,&flg);
99: /* Test MatMult() */
100: MatMultEqual(TA,B,10,&flg);
102: /* Test MatMultAdd() */
103: MatMultAddEqual(TA,B,10,&flg);
106: MatDestroy(&TA);
107: MatDestroy(&B);
109: /* Test KAIJ when T is is an identity matrix */
111: if (p == q) {
112: for (i=0; i<p; i++) {
113: for (j=0; j<q; j++) {
114: if (i==j) T[i+j*p] = 1.0;
115: else T[i+j*p] = 0.0;
116: }
117: }
119: /* Create KAIJ matrix TA */
120: MatCreateKAIJ(A,p,q,S,T,&TA);
121: MatGetLocalSize(A,&m,&n);
122: MatGetSize(A,&M,&N);
124: MatConvert(TA,MATAIJ,MAT_INITIAL_MATRIX,&B);
126: /* Test MatKAIJGetScaledIdentity() */
127: MatKAIJGetScaledIdentity(TA,&flg);
129: /* Test MatMult() */
130: MatMultEqual(TA,B,10,&flg);
132: /* Test MatMultAdd() */
133: MatMultAddEqual(TA,B,10,&flg);
136: MatDestroy(&TA);
137: MatDestroy(&B);
139: MatCreateKAIJ(A,p,q,NULL,T,&TA);
140: /* Test MatKAIJGetScaledIdentity() */
141: MatKAIJGetScaledIdentity(TA,&flg);
143: MatDestroy(&TA);
145: for (i=0; i<p; i++) {
146: for (j=0; j<q; j++) {
147: if (i==j) S[i+j*p] = T[i+j*p] = 2.0;
148: else S[i+j*p] = T[i+j*p] = 0.0;
149: }
150: }
151: MatCreateKAIJ(A,p,q,S,T,&TA);
152: /* Test MatKAIJGetScaledIdentity() */
153: MatKAIJGetScaledIdentity(TA,&flg);
155: MatDestroy(&TA);
156: }
158: /* Done with all tests */
160: PetscFree2(S,T);
161: MatDestroy(&A);
162: PetscFinalize();
163: return 0;
164: }
166: /*TEST
168: build:
169: requires: !complex
171: test:
172: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
173: output_file: output/ex176.out
174: nsize: {{1 2 3 4}}
175: args: -f ${DATAFILESPATH}/matrices/small -p {{2 3 7}} -q {{3 7}} -viewer_binary_skip_info
177: TEST*/