1: static char help[] = "Tests MatCreateConstantDiagonal().\n"
  2: "\n";

  4: #include <petscmat.h>

  6: /*T
  7:     Concepts: Mat
  8: T*/

 10: int main(int argc, char **args)
 11: {
 13:   Vec            X, Y;
 14:   Mat            A,B,Af;
 15:   PetscBool      flg;
 16:   PetscReal      xnorm,ynorm;

 18:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;

 20:   MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,20,20,3.0,&A);
 21:   MatCreateVecs(A,&X,&Y);
 22:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);

 24:   VecSetRandom(X,NULL);
 25:   VecNorm(X,NORM_2,&xnorm);
 26:   MatMult(A,X,Y);
 27:   VecNorm(Y,NORM_2,&ynorm);
 28:   if (PetscAbsReal(ynorm - 3*xnorm) > PETSC_SMALL) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Expected norm %g actual norm %g\n",(double)3*xnorm,(double)ynorm);
 29:   MatShift(A,5.0);
 30:   MatScale(A,.5);
 31:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);

 33:   /* Convert to AIJ (exercises MatGetRow/MatRestoreRow) */
 34:   MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&B);
 35:   MatMultEqual(A,B,10,&flg);
 36:   if (!flg) { PetscPrintf(PETSC_COMM_WORLD,"Error MatMult\n"); }
 37:   MatMultAddEqual(A,B,10,&flg);
 38:   if (!flg) { PetscPrintf(PETSC_COMM_WORLD,"Error MatMultAdd\n"); }
 39:   MatMultTransposeEqual(A,B,10,&flg);
 40:   if (!flg) { PetscPrintf(PETSC_COMM_WORLD,"Error MatMultTranspose\n"); }
 41:   MatMultTransposeAddEqual(A,B,10,&flg);
 42:   if (!flg) { PetscPrintf(PETSC_COMM_WORLD,"Error MatMultTransposeAdd\n"); }

 44:   MatGetDiagonal(A,Y);
 45:   MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&Af);
 46:   MatLUFactorSymbolic(Af,A,NULL,NULL,NULL);
 47:   MatLUFactorNumeric(Af,A,NULL);
 48:   MatSolve(Af,X,Y);
 49:   VecNorm(Y,NORM_2,&ynorm);
 50:   if (PetscAbsReal(ynorm - xnorm/4) > PETSC_SMALL) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Expected norm %g actual norm %g\n",(double).25*xnorm,(double)ynorm);

 52:   MatDestroy(&A);
 53:   MatDestroy(&B);
 54:   MatDestroy(&Af);
 55:   VecDestroy(&X);
 56:   VecDestroy(&Y);

 58:   PetscFinalize();
 59:   return ierr;
 60: }

 62: /*TEST

 64:   test:
 65:     nsize: 2

 67: TEST*/