1 static char help[] = "Test MatConvert() with a MATNEST with scaled and shifted MATTRANSPOSEVIRTUAL blocks.\n\n"; 2 3 #include <petscmat.h> 4 5 /* 6 This example builds the matrix 7 8 H = [ R C 9 alpha C^H + beta I gamma R^T + delta I ], 10 11 where R is Hermitian and C is complex symmetric. In particular, R and C have the 12 following Toeplitz structure: 13 14 R = pentadiag{a,b,c,conj(b),conj(a)} 15 C = tridiag{b,d,b} 16 17 where a,b,d are complex scalars, and c is real. 18 */ 19 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; 26 27 PetscFunctionBeginUser; 28 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 29 30 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 31 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); 36 37 PetscCall(MatCreate(PETSC_COMM_WORLD, &R)); 38 PetscCall(MatSetSizes(R, PETSC_DECIDE, PETSC_DECIDE, n, n)); 39 PetscCall(MatSetFromOptions(R)); 40 41 PetscCall(MatCreate(PETSC_COMM_WORLD, &C)); 42 PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, n, n)); 43 PetscCall(MatSetFromOptions(C)); 44 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 } 53 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 } 60 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)); 65 66 block[0] = R; 67 block[1] = C; 68 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])); 78 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"); 82 83 PetscCall(MatDestroy(&R)); 84 PetscCall(MatDestroy(&C)); 85 PetscCall(MatDestroy(&H)); 86 PetscCall(MatDestroy(&M)); 87 PetscCall(PetscFinalize()); 88 return 0; 89 } 90 91 /*TEST 92 93 build: 94 requires: complex 95 96 test: 97 output_file: output/ex109.out 98 nsize: {{1 4}} 99 100 TEST*/ 101