Actual source code: ex2.c
petsc-3.4.5 2014-06-29
2: static char help[] = "Tests MatTranspose(), MatNorm(), MatAXPY() and MatAYPX().\n\n";
4: #include <petscmat.h>
8: int main(int argc,char **argv)
9: {
10: Mat mat,tmat = 0;
11: PetscInt m = 7,n,i,j,rstart,rend,rect = 0;
13: PetscMPIInt size,rank;
14: PetscBool flg;
15: PetscScalar v, alpha;
16: PetscReal normf,normi,norm1;
18: PetscInitialize(&argc,&argv,(char*)0,help);
19: PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_COMMON);
20: PetscOptionsGetInt(NULL,"-m",&m,NULL);
21: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
22: MPI_Comm_size(PETSC_COMM_WORLD,&size);
23: n = m;
24: PetscOptionsHasName(NULL,"-rectA",&flg);
25: if (flg) {n += 2; rect = 1;}
26: PetscOptionsHasName(NULL,"-rectB",&flg);
27: if (flg) {n -= 2; rect = 1;}
29: /* ------- Assemble matrix, test MatValid() --------- */
30: MatCreate(PETSC_COMM_WORLD,&mat);
31: MatSetSizes(mat,PETSC_DECIDE,PETSC_DECIDE,m,n);
32: MatSetFromOptions(mat);
33: MatSetUp(mat);
34: MatGetOwnershipRange(mat,&rstart,&rend);
35: for (i=rstart; i<rend; i++) {
36: for (j=0; j<n; j++) {
37: v = 10.0*i+j;
38: MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
39: }
40: }
41: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
42: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
44: /* ----------------- Test MatNorm() ----------------- */
45: MatNorm(mat,NORM_FROBENIUS,&normf);
46: MatNorm(mat,NORM_1,&norm1);
47: MatNorm(mat,NORM_INFINITY,&normi);
48: PetscPrintf(PETSC_COMM_WORLD,"original A: Frobenious norm = %G, one norm = %G, infinity norm = %G\n",
49: normf,norm1,normi);
50: MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
52: /* --------------- Test MatTranspose() -------------- */
53: PetscOptionsHasName(NULL,"-in_place",&flg);
54: if (!rect && flg) {
55: MatTranspose(mat,MAT_REUSE_MATRIX,&mat); /* in-place transpose */
56: tmat = mat; mat = 0;
57: } else { /* out-of-place transpose */
58: MatTranspose(mat,MAT_INITIAL_MATRIX,&tmat);
59: }
61: /* ----------------- Test MatNorm() ----------------- */
62: /* Print info about transpose matrix */
63: MatNorm(tmat,NORM_FROBENIUS,&normf);
64: MatNorm(tmat,NORM_1,&norm1);
65: MatNorm(tmat,NORM_INFINITY,&normi);
66: PetscPrintf(PETSC_COMM_WORLD,"B = A^T: Frobenious norm = %G, one norm = %G, infinity norm = %G\n",
67: normf,norm1,normi);
68: MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);
70: /* ----------------- Test MatAXPY(), MatAYPX() ----------------- */
71: if (mat && !rect) {
72: alpha = 1.0;
73: PetscOptionsGetScalar(NULL,"-alpha",&alpha,NULL);
74: PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: B = B + alpha * A\n");
75: MatAXPY(tmat,alpha,mat,DIFFERENT_NONZERO_PATTERN);
76: MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);
78: PetscPrintf(PETSC_COMM_WORLD,"MatAYPX: B = alpha*B + A\n");
79: MatAYPX(tmat,alpha,mat,DIFFERENT_NONZERO_PATTERN);
80: MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);
81: }
83: {
84: Mat C;
85: alpha = 1.0;
86: PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: C = C + alpha * A, C=A, SAME_NONZERO_PATTERN\n");
87: MatDuplicate(mat,MAT_COPY_VALUES,&C);
88: MatAXPY(C,alpha,mat,SAME_NONZERO_PATTERN);
89: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
90: MatDestroy(&C);
91: }
93: {
94: Mat matB;
95: /* get matB that has nonzeros of mat in all even numbers of row and col */
96: MatCreate(PETSC_COMM_WORLD,&matB);
97: MatSetSizes(matB,PETSC_DECIDE,PETSC_DECIDE,m,n);
98: MatSetFromOptions(matB);
99: MatSetUp(matB);
100: MatGetOwnershipRange(matB,&rstart,&rend);
101: if (rstart % 2 != 0) rstart++;
102: for (i=rstart; i<rend; i += 2) {
103: for (j=0; j<n; j += 2) {
104: v = 10.0*i+j;
105: MatSetValues(matB,1,&i,1,&j,&v,INSERT_VALUES);
106: }
107: }
108: MatAssemblyBegin(matB,MAT_FINAL_ASSEMBLY);
109: MatAssemblyEnd(matB,MAT_FINAL_ASSEMBLY);
110: PetscPrintf(PETSC_COMM_WORLD," A: original matrix:\n");
111: MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
112: PetscPrintf(PETSC_COMM_WORLD," B(a subset of A):\n");
113: MatView(matB,PETSC_VIEWER_STDOUT_WORLD);
114: PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: B = B + alpha * A, SUBSET_NONZERO_PATTERN\n");
115: MatAXPY(mat,alpha,matB,SUBSET_NONZERO_PATTERN);
116: MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
117: MatDestroy(&matB);
118: }
120: /* Free data structures */
121: if (mat) {MatDestroy(&mat);}
122: if (tmat) {MatDestroy(&tmat);}
124: PetscFinalize();
125: return 0;
126: }