Actual source code: ex141.c


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

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

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

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

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

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

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

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

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

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

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

 98: /*TEST

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

103: TEST*/