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: