2: static char help[] = "Tests reusing MPI parallel matrices and MatGetValues().\n\
  3: To test the parallel matrix assembly, this example intentionally lays out\n\
  4: the matrix across processors differently from the way it is assembled.\n\
  5: This example uses bilinear elements on the unit square.  Input arguments are:\n\
  6:   -m <size> : problem size\n\n";

 8:  #include petscmat.h

 12: int FormElementStiffness(PetscReal H,PetscScalar *Ke)
 13: {
 15:   Ke[0]  = H/6.0;    Ke[1]  = -.125*H; Ke[2]  = H/12.0;   Ke[3]  = -.125*H;
 16:   Ke[4]  = -.125*H;  Ke[5]  = H/6.0;   Ke[6]  = -.125*H;  Ke[7]  = H/12.0;
 17:   Ke[8]  = H/12.0;   Ke[9]  = -.125*H; Ke[10] = H/6.0;    Ke[11] = -.125*H;
 18:   Ke[12] = -.125*H;  Ke[13] = H/12.0;  Ke[14] = -.125*H;  Ke[15] = H/6.0;
 19:   return(0);
 20: }

 24: int main(int argc,char **args)
 25: {
 26:   Mat            C;
 27:   Vec            u,b;
 29:   PetscMPIInt    size,rank;
 30:   PetscInt       i,m = 5,N,start,end,M,idx[4];
 31:   PetscInt       j,nrsub,ncsub,*rsub,*csub,mystart,myend;
 32:   PetscTruth     flg;
 33:   PetscScalar    one = 1.0,Ke[16],*vals;
 34:   PetscReal      h,norm;

 36:   PetscInitialize(&argc,&args,(char *)0,help);
 37:   PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);

 39:   N = (m+1)*(m+1); /* dimension of matrix */
 40:   M = m*m;         /* number of elements */
 41:   h = 1.0/m;       /* mesh width */
 42:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 43:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 45:   /* Create stiffness matrix */
 46:   MatCreate(PETSC_COMM_WORLD,&C);
 47:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,N,N);
 48:   MatSetFromOptions(C);

 50:   start = rank*(M/size) + ((M%size) < rank ? (M%size) : rank);
 51:   end   = start + M/size + ((M%size) > rank);

 53:   /* Form the element stiffness for the Laplacian */
 54:   FormElementStiffness(h*h,Ke);
 55:   for (i=start; i<end; i++) {
 56:      /* location of lower left corner of element */
 57:      /* node numbers for the four corners of element */
 58:      idx[0] = (m+1)*(i/m) + (i % m);
 59:      idx[1] = idx[0]+1; idx[2] = idx[1] + m + 1; idx[3] = idx[2] - 1;
 60:      MatSetValues(C,4,idx,4,idx,Ke,ADD_VALUES);
 61:   }
 62:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 63:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

 65:   /* Assemble the matrix again */
 66:   MatZeroEntries(C);

 68:   for (i=start; i<end; i++) {
 69:      /* location of lower left corner of element */
 70:      /* node numbers for the four corners of element */
 71:      idx[0] = (m+1)*(i/m) + (i % m);
 72:      idx[1] = idx[0]+1; idx[2] = idx[1] + m + 1; idx[3] = idx[2] - 1;
 73:      MatSetValues(C,4,idx,4,idx,Ke,ADD_VALUES);
 74:   }
 75:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 76:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

 78:   /* Create test vectors */
 79:   VecCreate(PETSC_COMM_WORLD,&u);
 80:   VecSetSizes(u,PETSC_DECIDE,N);
 81:   VecSetFromOptions(u);
 82:   VecDuplicate(u,&b);
 83:   VecSet(u,one);

 85:   /* Check error */
 86:   MatMult(C,u,b);
 87:   VecNorm(b,NORM_2,&norm);
 88:   if (norm > 1.e-10 || norm < -1.e-10) {
 89:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error b %G should be near 0\n",norm);
 90:   }

 92:   /* Now test MatGetValues() */
 93:   PetscOptionsHasName(PETSC_NULL,"-get_values",&flg);
 94:   if (flg) {
 95:     MatGetOwnershipRange(C,&mystart,&myend);
 96:     nrsub = myend - mystart; ncsub = 4;
 97:     PetscMalloc(nrsub*ncsub*sizeof(PetscScalar),&vals);
 98:     PetscMalloc(nrsub*sizeof(PetscInt),&rsub);
 99:     PetscMalloc(ncsub*sizeof(PetscInt),&csub);
100:     for (i=myend-1; i>=mystart; i--) rsub[myend-i-1] = i;
101:     for (i=0; i<ncsub; i++) csub[i] = 2*(ncsub-i) + mystart;
102:     MatGetValues(C,nrsub,rsub,ncsub,csub,vals);
103:     MatView(C,PETSC_VIEWER_STDOUT_WORLD);
104:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"processor number %d: start=%D, end=%D, mystart=%D, myend=%D\n",
105:             rank,start,end,mystart,myend);
106:     for (i=0; i<nrsub; i++) {
107:       for (j=0; j<ncsub; j++) {
108: #if defined(PETSC_USE_COMPLEX)
109:         if (PetscImaginaryPart(vals[i*ncsub+j]) != 0.0) {
110:            PetscSynchronizedPrintf(PETSC_COMM_WORLD,"  C[%D, %D] = %G + %G i\n",rsub[i],csub[j],PetscRealPart(vals[i*ncsub+j]),
111:                                        PetscImaginaryPart(vals[i*ncsub+j]));
112:         } else {
113:            PetscSynchronizedPrintf(PETSC_COMM_WORLD,"  C[%D, %D] = %G\n",rsub[i],csub[j],PetscRealPart(vals[i*ncsub+j]));
114:         }
115: #else
116:          PetscSynchronizedPrintf(PETSC_COMM_WORLD,"  C[%D, %D] = %G\n",rsub[i],csub[j],vals[i*ncsub+j]);
117: #endif
118:       }
119:     }
120:     PetscSynchronizedFlush(PETSC_COMM_WORLD);
121:     PetscFree(rsub);
122:     PetscFree(csub);
123:     PetscFree(vals);
124:   }

126:   /* Free data structures */
127:   VecDestroy(u);
128:   VecDestroy(b);
129:   MatDestroy(C);
130:   PetscFinalize();
131:   return 0;
132: }