Actual source code: ex15.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      2: static char help[] = "Tests MatNorm(), MatLUFactor(), MatSolve() and MatSolveAdd().\n\n";

      4:  #include <petscmat.h>

      6: int main(int argc,char **args)
      7: {
      8:   Mat            C;
      9:   PetscInt       i,j,m = 3,n = 3,Ii,J;
     11:   PetscBool      flg;
     12:   PetscScalar    v;
     13:   IS             perm,iperm;
     14:   Vec            x,u,b,y;
     15:   PetscReal      norm,tol=PETSC_SMALL;
     16:   MatFactorInfo  info;
     17:   PetscMPIInt    size;

     19:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
     20:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     21:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor example only!");
     22:   MatCreate(PETSC_COMM_WORLD,&C);
     23:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
     24:   MatSetFromOptions(C);
     25:   MatSetUp(C);
     26:   PetscOptionsHasName(NULL,NULL,"-symmetric",&flg);
     27:   if (flg) {  /* Treat matrix as symmetric only if we set this flag */
     28:     MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE);
     29:     MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
     30:   }

     32:   /* Create the matrix for the five point stencil, YET AGAIN */
     33:   for (i=0; i<m; i++) {
     34:     for (j=0; j<n; j++) {
     35:       v = -1.0;  Ii = j + n*i;
     36:       if (i>0)   {J = Ii - n; MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
     37:       if (i<m-1) {J = Ii + n; MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
     38:       if (j>0)   {J = Ii - 1; MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
     39:       if (j<n-1) {J = Ii + 1; MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
     40:       v = 4.0; MatSetValues(C,1,&Ii,1,&Ii,&v,INSERT_VALUES);
     41:     }
     42:   }
     43:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
     44:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
     45:   MatGetOrdering(C,MATORDERINGRCM,&perm,&iperm);
     46:   MatView(C,PETSC_VIEWER_STDOUT_WORLD);
     47:   ISView(perm,PETSC_VIEWER_STDOUT_SELF);
     48:   VecCreateSeq(PETSC_COMM_SELF,m*n,&u);
     49:   VecSet(u,1.0);
     50:   VecDuplicate(u,&x);
     51:   VecDuplicate(u,&b);
     52:   VecDuplicate(u,&y);
     53:   MatMult(C,u,b);
     54:   VecCopy(b,y);
     55:   VecScale(y,2.0);

     57:   MatNorm(C,NORM_FROBENIUS,&norm);
     58:   PetscPrintf(PETSC_COMM_SELF,"Frobenius norm of matrix %g\n",(double)norm);
     59:   MatNorm(C,NORM_1,&norm);
     60:   PetscPrintf(PETSC_COMM_SELF,"One  norm of matrix %g\n",(double)norm);
     61:   MatNorm(C,NORM_INFINITY,&norm);
     62:   PetscPrintf(PETSC_COMM_SELF,"Infinity norm of matrix %g\n",(double)norm);

     64:   MatFactorInfoInitialize(&info);
     65:   info.fill          = 2.0;
     66:   info.dtcol         = 0.0;
     67:   info.zeropivot     = 1.e-14;
     68:   info.pivotinblocks = 1.0;

     70:   MatLUFactor(C,perm,iperm,&info);

     72:   /* Test MatSolve */
     73:   MatSolve(C,b,x);
     74:   VecView(b,PETSC_VIEWER_STDOUT_SELF);
     75:   VecView(x,PETSC_VIEWER_STDOUT_SELF);
     76:   VecAXPY(x,-1.0,u);
     77:   VecNorm(x,NORM_2,&norm);
     78:   if (norm > tol) {
     79:     PetscPrintf(PETSC_COMM_SELF,"MatSolve: Norm of error %g\n",(double)norm);
     80:   }

     82:   /* Test MatSolveAdd */
     83:   MatSolveAdd(C,b,y,x);
     84:   VecAXPY(x,-1.0,y);
     85:   VecAXPY(x,-1.0,u);
     86:   VecNorm(x,NORM_2,&norm);
     87:   if (norm > tol) {
     88:     PetscPrintf(PETSC_COMM_SELF,"MatSolveAdd(): Norm of error %g\n",(double)norm);
     89:   }

     91:   ISDestroy(&perm);
     92:   ISDestroy(&iperm);
     93:   VecDestroy(&u);
     94:   VecDestroy(&y);
     95:   VecDestroy(&b);
     96:   VecDestroy(&x);
     97:   MatDestroy(&C);
     98:   PetscFinalize();
     99:   return ierr;
    100: }

    103: /*TEST

    105:    test:

    107: TEST*/