1 2 static char help[] = "Saves 4by4 block matrix.\n\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc, char **args) 7 { 8 Mat A; 9 PetscInt i, j; 10 PetscMPIInt size; 11 PetscViewer fd; 12 PetscScalar values[16], one = 1.0; 13 Vec x; 14 15 PetscFunctionBeginUser; 16 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 17 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 18 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, m "Can only run on one processor"); 19 20 /* 21 Open binary file. Note that we use FILE_MODE_WRITE to indicate 22 writing to this file. 23 */ 24 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "4by4", FILE_MODE_WRITE, &fd)); 25 26 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_WORLD, 4, 12, 12, 0, 0, &A)); 27 28 for (i = 0; i < 16; i++) values[i] = i; 29 for (i = 0; i < 4; i++) values[4 * i + i] += 5; 30 i = 0; 31 j = 0; 32 PetscCall(MatSetValuesBlocked(A, 1, &i, 1, &j, values, INSERT_VALUES)); 33 34 for (i = 0; i < 16; i++) values[i] = i; 35 i = 0; 36 j = 2; 37 PetscCall(MatSetValuesBlocked(A, 1, &i, 1, &j, values, INSERT_VALUES)); 38 39 for (i = 0; i < 16; i++) values[i] = i; 40 i = 1; 41 j = 0; 42 PetscCall(MatSetValuesBlocked(A, 1, &i, 1, &j, values, INSERT_VALUES)); 43 44 for (i = 0; i < 16; i++) values[i] = i; 45 for (i = 0; i < 4; i++) values[4 * i + i] += 6; 46 i = 1; 47 j = 1; 48 PetscCall(MatSetValuesBlocked(A, 1, &i, 1, &j, values, INSERT_VALUES)); 49 50 for (i = 0; i < 16; i++) values[i] = i; 51 i = 2; 52 j = 0; 53 PetscCall(MatSetValuesBlocked(A, 1, &i, 1, &j, values, INSERT_VALUES)); 54 55 for (i = 0; i < 16; i++) values[i] = i; 56 for (i = 0; i < 4; i++) values[4 * i + i] += 7; 57 i = 2; 58 j = 2; 59 PetscCall(MatSetValuesBlocked(A, 1, &i, 1, &j, values, INSERT_VALUES)); 60 61 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 62 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 63 PetscCall(MatView(A, fd)); 64 PetscCall(MatDestroy(&A)); 65 66 PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 12, &x)); 67 PetscCall(VecSet(x, one)); 68 PetscCall(VecView(x, fd)); 69 PetscCall(VecDestroy(&x)); 70 71 PetscCall(PetscViewerDestroy(&fd)); 72 PetscCall(PetscFinalize()); 73 return 0; 74 } 75