Actual source code: ex264.c

  1: static char help[] = "Test MatConvert() with a MATNEST with scaled and shifted MATTRANSPOSEVIRTUAL blocks.\n\n";

  3: #include <petscmat.h>

  5: /*
  6:    This example builds the matrix

  8:         H = [  R                    C
  9:               alpha C^H + beta I    gamma R^T + delta I ],

 11:    where R is Hermitian and C is complex symmetric. In particular, R and C have the
 12:    following Toeplitz structure:

 14:         R = pentadiag{a,b,c,conj(b),conj(a)}
 15:         C = tridiag{b,d,b}

 17:    where a,b,d are complex scalars, and c is real.
 18: */

 20: int main(int argc, char **argv)
 21: {
 22:   Mat         block[4], H, R, C, M;
 23:   PetscScalar a, b, c, d;
 24:   PetscInt    n = 13, Istart, Iend, i;
 25:   PetscBool   flg;

 27:   PetscFunctionBeginUser;
 28:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

 30:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));

 32:   a = PetscCMPLX(-0.1, 0.2);
 33:   b = PetscCMPLX(1.0, 0.5);
 34:   c = 4.5;
 35:   d = PetscCMPLX(2.0, 0.2);

 37:   PetscCall(MatCreate(PETSC_COMM_WORLD, &R));
 38:   PetscCall(MatSetSizes(R, PETSC_DECIDE, PETSC_DECIDE, n, n));
 39:   PetscCall(MatSetFromOptions(R));

 41:   PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
 42:   PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, n, n));
 43:   PetscCall(MatSetFromOptions(C));

 45:   PetscCall(MatGetOwnershipRange(R, &Istart, &Iend));
 46:   for (i = Istart; i < Iend; i++) {
 47:     if (i > 1) PetscCall(MatSetValue(R, i, i - 2, a, INSERT_VALUES));
 48:     if (i > 0) PetscCall(MatSetValue(R, i, i - 1, b, INSERT_VALUES));
 49:     PetscCall(MatSetValue(R, i, i, c, INSERT_VALUES));
 50:     if (i < n - 1) PetscCall(MatSetValue(R, i, i + 1, PetscConj(b), INSERT_VALUES));
 51:     if (i < n - 2) PetscCall(MatSetValue(R, i, i + 2, PetscConj(a), INSERT_VALUES));
 52:   }

 54:   PetscCall(MatGetOwnershipRange(C, &Istart, &Iend));
 55:   for (i = Istart; i < Iend; i++) {
 56:     if (i > 0) PetscCall(MatSetValue(C, i, i - 1, b, INSERT_VALUES));
 57:     PetscCall(MatSetValue(C, i, i, d, INSERT_VALUES));
 58:     if (i < n - 1) PetscCall(MatSetValue(C, i, i + 1, b, INSERT_VALUES));
 59:   }

 61:   PetscCall(MatAssemblyBegin(R, MAT_FINAL_ASSEMBLY));
 62:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
 63:   PetscCall(MatAssemblyEnd(R, MAT_FINAL_ASSEMBLY));
 64:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));

 66:   block[0] = R;
 67:   block[1] = C;

 69:   PetscCall(MatCreateHermitianTranspose(C, &block[2]));
 70:   PetscCall(MatScale(block[2], PetscConj(b)));
 71:   PetscCall(MatShift(block[2], d));
 72:   PetscCall(MatCreateTranspose(R, &block[3]));
 73:   PetscCall(MatScale(block[3], PetscConj(d)));
 74:   PetscCall(MatShift(block[3], b));
 75:   PetscCall(MatCreateNest(PetscObjectComm((PetscObject)R), 2, NULL, 2, NULL, block, &H));
 76:   PetscCall(MatDestroy(&block[2]));
 77:   PetscCall(MatDestroy(&block[3]));

 79:   PetscCall(MatConvert(H, MATAIJ, MAT_INITIAL_MATRIX, &M));
 80:   PetscCall(MatMultEqual(H, M, 20, &flg));
 81:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "MatNest != MatAIJ");

 83:   PetscCall(MatDestroy(&R));
 84:   PetscCall(MatDestroy(&C));
 85:   PetscCall(MatDestroy(&H));
 86:   PetscCall(MatDestroy(&M));
 87:   PetscCall(PetscFinalize());
 88:   return 0;
 89: }

 91: /*TEST

 93:    build:
 94:       requires: complex

 96:    test:
 97:       output_file: output/ex109.out
 98:       nsize: {{1 4}}

100: TEST*/