Actual source code: ex19.c
1: /*$Id: ex19.c,v 1.30 2001/08/07 21:30:08 bsmith Exp $*/
3: static char help[] = "Tests reusing MPI parallel matrices and MatGetValues().\n\
4: To test the parallel matrix assembly, this example intentionally lays out\n\
5: the matrix across processors differently from the way it is assembled.\n\
6: This example uses bilinear elements on the unit square. Input arguments are:\n\
7: -m <size> : problem size\n\n";
9: #include petscmat.h
13: int FormElementStiffness(PetscReal H,PetscScalar *Ke)
14: {
16: Ke[0] = H/6.0; Ke[1] = -.125*H; Ke[2] = H/12.0; Ke[3] = -.125*H;
17: Ke[4] = -.125*H; Ke[5] = H/6.0; Ke[6] = -.125*H; Ke[7] = H/12.0;
18: Ke[8] = H/12.0; Ke[9] = -.125*H; Ke[10] = H/6.0; Ke[11] = -.125*H;
19: Ke[12] = -.125*H; Ke[13] = H/12.0; Ke[14] = -.125*H; Ke[15] = H/6.0;
20: return(0);
21: }
25: int main(int argc,char **args)
26: {
27: Mat C;
28: Vec u,b;
29: int i,m = 5,rank,size,N,start,end,M,ierr,idx[4];
30: int j,nrsub,ncsub,*rsub,*csub,mystart,myend;
31: PetscTruth flg;
32: PetscScalar one = 1.0,Ke[16],*vals;
33: PetscReal h,norm;
35: PetscInitialize(&argc,&args,(char *)0,help);
36: PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
38: N = (m+1)*(m+1); /* dimension of matrix */
39: M = m*m; /* number of elements */
40: h = 1.0/m; /* mesh width */
41: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
42: MPI_Comm_size(PETSC_COMM_WORLD,&size);
44: /* Create stiffness matrix */
45: MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,&C);
46: MatSetFromOptions(C);
48: start = rank*(M/size) + ((M%size) < rank ? (M%size) : rank);
49: end = start + M/size + ((M%size) > rank);
51: /* Form the element stiffness for the Laplacian */
52: FormElementStiffness(h*h,Ke);
53: for (i=start; i<end; i++) {
54: /* location of lower left corner of element */
55: /* node numbers for the four corners of element */
56: idx[0] = (m+1)*(i/m) + (i % m);
57: idx[1] = idx[0]+1; idx[2] = idx[1] + m + 1; idx[3] = idx[2] - 1;
58: MatSetValues(C,4,idx,4,idx,Ke,ADD_VALUES);
59: }
60: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
61: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
63: /* Assemble the matrix again */
64: MatZeroEntries(C);
66: for (i=start; i<end; i++) {
67: /* location of lower left corner of element */
68: /* node numbers for the four corners of element */
69: idx[0] = (m+1)*(i/m) + (i % m);
70: idx[1] = idx[0]+1; idx[2] = idx[1] + m + 1; idx[3] = idx[2] - 1;
71: MatSetValues(C,4,idx,4,idx,Ke,ADD_VALUES);
72: }
73: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
74: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
76: /* Create test vectors */
77: VecCreate(PETSC_COMM_WORLD,&u);
78: VecSetSizes(u,PETSC_DECIDE,N);
79: VecSetFromOptions(u);
80: VecDuplicate(u,&b);
81: VecSet(&one,u);
83: /* Check error */
84: MatMult(C,u,b);
85: VecNorm(b,NORM_2,&norm);
86: if (norm > 1.e-10 || norm < -1.e-10) {
87: PetscPrintf(PETSC_COMM_WORLD,"Norm of error b %g should be near 0\n",norm);
88: }
90: /* Now test MatGetValues() */
91: PetscOptionsHasName(PETSC_NULL,"-get_values",&flg);
92: if (flg) {
93: MatGetOwnershipRange(C,&mystart,&myend);
94: nrsub = myend - mystart; ncsub = 4;
95: PetscMalloc(nrsub*ncsub*sizeof(PetscScalar),&vals);
96: PetscMalloc(nrsub*sizeof(int),&rsub);
97: PetscMalloc(ncsub*sizeof(int),&csub);
98: for (i=myend-1; i>=mystart; i--) rsub[myend-i-1] = i;
99: for (i=0; i<ncsub; i++) csub[i] = 2*(ncsub-i) + mystart;
100: MatGetValues(C,nrsub,rsub,ncsub,csub,vals);
101: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
102: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"processor number %d: start=%d, end=%d, mystart=%d, myend=%d\n",
103: rank,start,end,mystart,myend);
104: for (i=0; i<nrsub; i++) {
105: for (j=0; j<ncsub; j++) {
106: #if defined(PETSC_USE_COMPLEX)
107: if (PetscImaginaryPart(vals[i*ncsub+j]) != 0.0) {
108: PetscSynchronizedPrintf(PETSC_COMM_WORLD," C[%d, %d] = %g + %g i\n",rsub[i],csub[j],PetscRealPart(vals[i*ncsub+j]),
109: PetscImaginaryPart(vals[i*ncsub+j]));
110: } else {
111: PetscSynchronizedPrintf(PETSC_COMM_WORLD," C[%d, %d] = %g\n",rsub[i],csub[j],PetscRealPart(vals[i*ncsub+j]));
112: }
113: #else
114: PetscSynchronizedPrintf(PETSC_COMM_WORLD," C[%d, %d] = %g\n",rsub[i],csub[j],vals[i*ncsub+j]);
115: #endif
116: }
117: }
118: PetscSynchronizedFlush(PETSC_COMM_WORLD);
119: PetscFree(rsub);
120: PetscFree(csub);
121: PetscFree(vals);
122: }
124: /* Free data structures */
125: VecDestroy(u);
126: VecDestroy(b);
127: MatDestroy(C);
128: PetscFinalize();
129: return 0;
130: }