Actual source code: ex2.c
2: static char help[] = "Tests MatTranspose(), MatNorm(), MatValid(), 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;
12: PetscErrorCode ierr;
13: PetscMPIInt size,rank;
14: PetscTruth 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(PETSC_NULL,"-m",&m,PETSC_NULL);
21: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
22: MPI_Comm_size(PETSC_COMM_WORLD,&size);
23: n = m;
24: PetscOptionsHasName(PETSC_NULL,"-rectA",&flg);
25: if (flg) {n += 2; rect = 1;}
26: PetscOptionsHasName(PETSC_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: MatGetOwnershipRange(mat,&rstart,&rend);
34: for (i=rstart; i<rend; i++) {
35: for (j=0; j<n; j++) {
36: v=10.0*i+j;
37: MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
38: }
39: }
40: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
41: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
43: /* Test whether matrix has been corrupted (just to demonstrate this
44: routine) not needed in most application codes. */
45: MatValid(mat,(PetscTruth*)&flg);
46: if (!flg) SETERRQ(1,"Corrupted matrix.");
48: /* ----------------- Test MatNorm() ----------------- */
49: MatNorm(mat,NORM_FROBENIUS,&normf);
50: MatNorm(mat,NORM_1,&norm1);
51: MatNorm(mat,NORM_INFINITY,&normi);
52: PetscPrintf(PETSC_COMM_WORLD,"original A: Frobenious norm = %G, one norm = %G, infinity norm = %G\n",
53: normf,norm1,normi);
54: MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
56: /* --------------- Test MatTranspose() -------------- */
57: PetscOptionsHasName(PETSC_NULL,"-in_place",&flg);
58: if (!rect && flg) {
59: MatTranspose(mat,0); /* in-place transpose */
60: tmat = mat; mat = 0;
61: } else { /* out-of-place transpose */
62: MatTranspose(mat,&tmat);
63: }
65: /* ----------------- Test MatNorm() ----------------- */
66: /* Print info about transpose matrix */
67: MatNorm(tmat,NORM_FROBENIUS,&normf);
68: MatNorm(tmat,NORM_1,&norm1);
69: MatNorm(tmat,NORM_INFINITY,&normi);
70: PetscPrintf(PETSC_COMM_WORLD,"B = A^T: Frobenious norm = %G, one norm = %G, infinity norm = %G\n",
71: normf,norm1,normi);
72: MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);
74: /* ----------------- Test MatAXPY(), MatAYPX() ----------------- */
75: if (mat && !rect) {
76: alpha = 1.0;
77: PetscOptionsGetScalar(PETSC_NULL,"-alpha",&alpha,PETSC_NULL);
78: PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: B = B + alpha * A\n");
79: MatAXPY(tmat,alpha,mat,DIFFERENT_NONZERO_PATTERN);
80: MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);
82: PetscPrintf(PETSC_COMM_WORLD,"MatAYPX: B = alpha*B + A\n");
83: MatAYPX(tmat,alpha,mat,DIFFERENT_NONZERO_PATTERN);
84: MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);
85: }
87: {
88: Mat C;
89: alpha = 1.0;
90: PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: C = C + alpha * A, C=A, SAME_NONZERO_PATTERN\n");
91: MatDuplicate(mat,MAT_COPY_VALUES,&C);
92: MatAXPY(C,alpha,mat,SAME_NONZERO_PATTERN);
93: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
94: MatDestroy(C);
95: }
97: {
98: Mat matB;
99: /* get matB that has nonzeros of mat in all even numbers of row and col */
100: MatCreate(PETSC_COMM_WORLD,&matB);
101: MatSetSizes(matB,PETSC_DECIDE,PETSC_DECIDE,m,n);
102: MatSetFromOptions(matB);
103: MatGetOwnershipRange(matB,&rstart,&rend);
104: if (rstart % 2 != 0) rstart++;
105: for (i=rstart; i<rend; i += 2) {
106: for (j=0; j<n; j += 2) {
107: v=10.0*i+j;
108: MatSetValues(matB,1,&i,1,&j,&v,INSERT_VALUES);
109: }
110: }
111: MatAssemblyBegin(matB,MAT_FINAL_ASSEMBLY);
112: MatAssemblyEnd(matB,MAT_FINAL_ASSEMBLY);
113: PetscPrintf(PETSC_COMM_WORLD," A: original matrix:\n");
114: MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
115: PetscPrintf(PETSC_COMM_WORLD," B(a subset of A):\n");
116: MatView(matB,PETSC_VIEWER_STDOUT_WORLD);
117: PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: B = B + alpha * A, SUBSET_NONZERO_PATTERN\n");
118: MatAXPY(mat,alpha,matB,SUBSET_NONZERO_PATTERN);
119: MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
120: MatDestroy(matB);
121: }
123: /* Free data structures */
124: if (mat) {MatDestroy(mat);}
125: if (tmat) {MatDestroy(tmat);}
127: PetscFinalize();
128: return 0;
129: }
130: