Actual source code: ex141.c

  1: static char help[] = "Tests converting a SBAIJ matrix to BAIJ format with MatConvert. Modified from ex55.c\n\n";

  3: #include <petscmat.h>
  4: /* Usage: ./ex141 -mat_view */

  6: int main(int argc, char **args)
  7: {
  8:   Mat         C, B;
  9:   PetscInt    i, bs = 2, mbs, m, block, d_nz = 6, col[3];
 10:   PetscMPIInt size;
 11:   char        file[PETSC_MAX_PATH_LEN];
 12:   PetscViewer fd;
 13:   PetscBool   equal, loadmat;
 14:   PetscScalar value[4];
 15:   PetscInt    ridx[2], cidx[2];

 17:   PetscFunctionBeginUser;
 18:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 19:   PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &loadmat));
 20:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 21:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

 23:   /* input matrix C */
 24:   if (loadmat) {
 25:     /* Open binary file. Load a sbaij matrix, then destroy the viewer. */
 26:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
 27:     PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
 28:     PetscCall(MatSetType(C, MATSEQSBAIJ));
 29:     PetscCall(MatLoad(C, fd));
 30:     PetscCall(PetscViewerDestroy(&fd));
 31:   } else { /* Create a sbaij mat with bs>1  */
 32:     mbs = 8;
 33:     PetscCall(PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL));
 34:     m = mbs * bs;
 35:     PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
 36:     PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, m, m));
 37:     PetscCall(MatSetType(C, MATSBAIJ));
 38:     PetscCall(MatSetFromOptions(C));
 39:     PetscCall(MatSeqSBAIJSetPreallocation(C, bs, d_nz, NULL));
 40:     PetscCall(MatSetUp(C));
 41:     PetscCall(MatSetOption(C, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE));

 43:     for (block = 0; block < mbs; block++) {
 44:       /* diagonal blocks */
 45:       value[0] = -1.0;
 46:       value[1] = 4.0;
 47:       value[2] = -1.0;
 48:       for (i = 1 + block * bs; i < bs - 1 + block * bs; i++) {
 49:         col[0] = i - 1;
 50:         col[1] = i;
 51:         col[2] = i + 1;
 52:         PetscCall(MatSetValues(C, 1, &i, 3, col, value, INSERT_VALUES));
 53:       }
 54:       i      = bs - 1 + block * bs;
 55:       col[0] = bs - 2 + block * bs;
 56:       col[1] = bs - 1 + block * bs;

 58:       value[0] = -1.0;
 59:       value[1] = 4.0;
 60:       PetscCall(MatSetValues(C, 1, &i, 2, col, value, INSERT_VALUES));

 62:       i      = 0 + block * bs;
 63:       col[0] = 0 + block * bs;
 64:       col[1] = 1 + block * bs;

 66:       value[0] = 4.0;
 67:       value[1] = -1.0;
 68:       PetscCall(MatSetValues(C, 1, &i, 2, col, value, INSERT_VALUES));
 69:     }
 70:     /* off-diagonal blocks */
 71:     value[0] = -1.0;
 72:     value[1] = -0.1;
 73:     value[2] = 0.0;
 74:     value[3] = -1.0; /* row-oriented */
 75:     for (block = 0; block < mbs - 1; block++) {
 76:       for (i = 0; i < bs; i++) {
 77:         ridx[i] = block * bs + i;
 78:         cidx[i] = (block + 1) * bs + i;
 79:       }
 80:       PetscCall(MatSetValues(C, bs, ridx, bs, cidx, value, INSERT_VALUES));
 81:     }
 82:     PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
 83:     PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
 84:   }

 86:   /* convert C to BAIJ format */
 87:   PetscCall(MatConvert(C, MATSEQBAIJ, MAT_INITIAL_MATRIX, &B));
 88:   PetscCall(MatMultEqual(B, C, 10, &equal));
 89:   PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatConvert fails!");

 91:   PetscCall(MatDestroy(&B));
 92:   PetscCall(MatDestroy(&C));
 93:   PetscCall(PetscFinalize());
 94:   return 0;
 95: }

 97: /*TEST

 99:    test:
100:      output_file: output/ex141.out

102: TEST*/