Actual source code: ex2.c
1: /*$Id: ex2.c,v 1.25 2001/08/07 21:30:08 bsmith Exp $*/
3: static char help[] = "Tests MatTranspose(), MatNorm(), MatValid(), and MatAXPY().\n\n";
5: #include petscmat.h
7: #define TestMatNorm
8: #define TestMatTranspose
9: #define TestMatNorm_tmat
10: #define TestMatAXPY
11: #undef TestMatAXPY2
15: int main(int argc,char **argv)
16: {
17: Mat mat,tmat = 0;
18: int m = 7,n,i,j,ierr,size,rank,rstart,rend,rect = 0;
19: PetscTruth flg;
20: PetscScalar v, alpha;
21: PetscReal normf,normi,norm1;
23: PetscInitialize(&argc,&argv,(char*)0,help);
24: PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_COMMON);
25: PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
26: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
27: MPI_Comm_size(PETSC_COMM_WORLD,&size);
28: n = m;
29: PetscOptionsHasName(PETSC_NULL,"-rectA",&flg);
30: if (flg) {n += 2; rect = 1;}
31: PetscOptionsHasName(PETSC_NULL,"-rectB",&flg);
32: if (flg) {n -= 2; rect = 1;}
34: /* ------- Assemble matrix, test MatValid() --------- */
36: MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,m,n,&mat);
37: MatSetFromOptions(mat);
38: MatGetOwnershipRange(mat,&rstart,&rend);
39: for (i=rstart; i<rend; i++) {
40: for (j=0; j<n; j++) {
41: v=10*i+j;
42: MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
43: }
44: }
45: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
46: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
48: /* Test whether matrix has been corrupted (just to demonstrate this
49: routine) not needed in most application codes. */
50: MatValid(mat,(PetscTruth*)&flg);
51: if (!flg) SETERRQ(1,"Corrupted matrix.");
53: /* ----------------- Test MatNorm() ----------------- */
54: #ifdef TestMatNorm
55: MatNorm(mat,NORM_FROBENIUS,&normf);
56: MatNorm(mat,NORM_1,&norm1);
57: MatNorm(mat,NORM_INFINITY,&normi);
58: PetscPrintf(PETSC_COMM_WORLD,"original: Frobenious norm = %g, one norm = %g, infinity norm = %g\n",
59: normf,norm1,normi);
60: MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
61: #endif
62: /* --------------- Test MatTranspose() -------------- */
63: #ifdef TestMatTranspose
64: PetscOptionsHasName(PETSC_NULL,"-in_place",&flg);
65: if (!rect && flg) {
66: MatTranspose(mat,0); /* in-place transpose */
67: tmat = mat; mat = 0;
68: } else { /* out-of-place transpose */
69: MatTranspose(mat,&tmat);
70: }
71: #endif
72: /* ----------------- Test MatNorm() ----------------- */
73: #ifdef TestMatNorm_tmat
74: /* Print info about transpose matrix */
75: MatNorm(tmat,NORM_FROBENIUS,&normf);
76: MatNorm(tmat,NORM_1,&norm1);
77: MatNorm(tmat,NORM_INFINITY,&normi);
78: PetscPrintf(PETSC_COMM_WORLD,"transpose: Frobenious norm = %g, one norm = %g, infinity norm = %g\n",
79: normf,norm1,normi);
80: MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);
81: #endif
82: /* ----------------- Test MatAXPY() ----------------- */
83: #ifdef TestMatAXPY
84: if (mat && !rect) {
85: alpha = 1.0;
86: PetscOptionsGetScalar(PETSC_NULL,"-alpha",&alpha,PETSC_NULL);
87: PetscPrintf(PETSC_COMM_WORLD,"matrix addition: B = B + alpha * A\n");
88: MatAXPY(&alpha,mat,tmat,DIFFERENT_NONZERO_PATTERN);
89: MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);
90: }
91: MatDestroy(tmat);
92: #endif
93: #ifdef TestMatAXPY2
94: Mat matB;
95: alpha = 1.0;
96: PetscPrintf(PETSC_COMM_WORLD,"matrix addition: B = B + alpha * A, SAME_NONZERO_PATTERN\n");
97: MatDuplicate(mat,MAT_COPY_VALUES,&matB);
98: MatAXPY(&alpha,mat,matB,SAME_NONZERO_PATTERN);
99: MatView(matB,PETSC_VIEWER_STDOUT_WORLD);
100: MatDestroy(matB);
102: /* get matB that has nonzeros of mat in all even numbers of row and col */
103: MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,m,n,&matB);
104: MatSetFromOptions(matB);
105: MatGetOwnershipRange(matB,&rstart,&rend);
106: for (i=rstart; i<rend; i += 2) {
107: for (j=0; j<n; j += 2) {
108: v=10*i+j;
109: MatSetValues(matB,1,&i,1,&j,&v,INSERT_VALUES);
110: }
111: }
112: MatAssemblyBegin(matB,MAT_FINAL_ASSEMBLY);
113: MatAssemblyEnd(matB,MAT_FINAL_ASSEMBLY);
114: PetscPrintf(PETSC_COMM_WORLD," B: original matrix:\n");
115: MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
116: PetscPrintf(PETSC_COMM_WORLD," A(a subset of B):\n");
117: MatView(matB,PETSC_VIEWER_STDOUT_WORLD);
118: PetscPrintf(PETSC_COMM_WORLD,"matrix addition: B = B + alpha * A, SUBSET_NONZERO_PATTERN\n");
119: MatAXPY(&alpha,matB,mat,SUBSET_NONZERO_PATTERN);
120: MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
121:
122: PetscPrintf(PETSC_COMM_WORLD,"matrix addition: B = B + alpha * A, SUBSET_NONZERO_PATTERN\n");
123: MatAXPY(&alpha,matB,mat,SUBSET_NONZERO_PATTERN);
124: MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
125:
126: MatDestroy(matB);
127: #endif
128:
130: /* Free data structures */
131: if (mat) {MatDestroy(mat);}
133: PetscFinalize();
134: return 0;
135: }
136: