Actual source code: ex185.c

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

  4: #include <petscmat.h>

  6: int main(int argc, char **args)
  7: {
  8:   Vec       X, Y;
  9:   Mat       A, Adup, B, Af;
 10:   PetscBool flg;
 11:   PetscReal xnorm, ynorm, anorm;

 13:   PetscFunctionBeginUser;
 14:   PetscCall(PetscInitialize(&argc, &args, NULL, help));

 16:   PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, 20, 20, 3.0, &A));
 17:   PetscCall(MatCreateVecs(A, &X, &Y));
 18:   PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));

 20:   PetscCall(VecSetRandom(X, NULL));
 21:   PetscCall(VecNorm(X, NORM_2, &xnorm));
 22:   PetscCall(MatMult(A, X, Y));
 23:   PetscCall(VecNorm(Y, NORM_2, &ynorm));
 24:   PetscCheck(PetscAbsReal(ynorm - 3 * xnorm) <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Expected norm %g actual norm %g", (double)(3 * xnorm), (double)ynorm);
 25:   PetscCall(MatShift(A, 5.0));
 26:   PetscCall(MatScale(A, .5));
 27:   PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
 28:   PetscCall(MatNorm(A, NORM_FROBENIUS, &anorm));
 29:   PetscCheck(PetscAbsReal(anorm - 4) <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Expected norm 4.0 actual norm %g", (double)anorm);

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

 46:   PetscCall(MatGetDiagonal(A, Y));
 47:   PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &Af));
 48:   PetscCall(MatLUFactorSymbolic(Af, A, NULL, NULL, NULL));
 49:   PetscCall(MatLUFactorNumeric(Af, A, NULL));
 50:   PetscCall(MatSolve(Af, X, Y));
 51:   PetscCall(VecNorm(Y, NORM_2, &ynorm));
 52:   PetscCheck(PetscAbsReal(ynorm - xnorm / 4) <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Expected norm %g actual norm %g", (double)(.25 * xnorm), (double)ynorm);

 54:   // Solve can be called without factorization
 55:   PetscCall(MatSolve(A, X, Y));
 56:   PetscCall(VecNorm(Y, NORM_2, &ynorm));
 57:   PetscCheck(PetscAbsReal(ynorm - xnorm / 4) <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Expected norm %g actual norm %g", (double)(.25 * xnorm), (double)ynorm);

 59:   // For a scalar multiple of the identity  smoothing is equivalent to solving
 60:   PetscCall(MatSOR(A, X, 1.5, SOR_FORWARD_SWEEP, 0.0, 1, 1, Y));
 61:   PetscCall(VecNorm(Y, NORM_2, &ynorm));
 62:   PetscCheck(PetscAbsReal(ynorm - xnorm / 4) <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Expected norm %g actual norm %g", (double)(.25 * xnorm), (double)ynorm);

 64:   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &Adup));
 65:   PetscCall(MatEqual(A, Adup, &flg));
 66:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "MatEqual after copy failure");

 68:   PetscCall(MatDestroy(&Adup));
 69:   PetscCall(MatDestroy(&A));
 70:   PetscCall(MatDestroy(&B));
 71:   PetscCall(MatDestroy(&Af));
 72:   PetscCall(VecDestroy(&X));
 73:   PetscCall(VecDestroy(&Y));

 75:   PetscCall(PetscFinalize());
 76:   return 0;
 77: }

 79: /*TEST

 81:   test:
 82:     nsize: 2

 84: TEST*/