1 2 static char help[] = "Reads a matrix from PETSc binary file. Use for view or investigating matrix data structure. \n\n"; 3 /* 4 Example: 5 ./ex16 -f <matrix file> -a_mat_view draw -draw_pause -1 6 ./ex16 -f <matrix file> -a_mat_view ascii::ascii_info 7 */ 8 9 #include <petscmat.h> 10 int main(int argc, char **args) 11 { 12 Mat A, Asp; 13 PetscViewer fd; /* viewer */ 14 char file[PETSC_MAX_PATH_LEN]; /* input file name */ 15 PetscInt m, n, rstart, rend; 16 PetscBool flg; 17 PetscInt row, ncols, j, nrows, nnzA = 0, nnzAsp = 0; 18 const PetscInt *cols; 19 const PetscScalar *vals; 20 PetscReal norm, percent, val, dtol = 1.e-16; 21 PetscMPIInt rank; 22 MatInfo matinfo; 23 PetscInt Dnnz, Onnz; 24 25 PetscFunctionBeginUser; 26 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 27 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 28 29 /* Determine files from which we read the linear systems. */ 30 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg)); 31 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option"); 32 33 /* Open binary file. Note that we use FILE_MODE_READ to indicate 34 reading from this file. */ 35 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd)); 36 37 /* Load the matrix; then destroy the viewer. */ 38 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 39 PetscCall(MatSetOptionsPrefix(A, "a_")); 40 PetscCall(MatSetFromOptions(A)); 41 PetscCall(MatLoad(A, fd)); 42 PetscCall(PetscViewerDestroy(&fd)); 43 PetscCall(MatGetSize(A, &m, &n)); 44 PetscCall(MatGetInfo(A, MAT_LOCAL, &matinfo)); 45 /*printf("matinfo.nz_used %g\n",matinfo.nz_used);*/ 46 47 /* Get a sparse matrix Asp by dumping zero entries of A */ 48 PetscCall(MatCreate(PETSC_COMM_WORLD, &Asp)); 49 PetscCall(MatSetSizes(Asp, m, n, PETSC_DECIDE, PETSC_DECIDE)); 50 PetscCall(MatSetOptionsPrefix(Asp, "asp_")); 51 PetscCall(MatSetFromOptions(Asp)); 52 Dnnz = (PetscInt)matinfo.nz_used / m + 1; 53 Onnz = Dnnz / 2; 54 printf("Dnnz %d %d\n", Dnnz, Onnz); 55 PetscCall(MatSeqAIJSetPreallocation(Asp, Dnnz, NULL)); 56 PetscCall(MatMPIAIJSetPreallocation(Asp, Dnnz, NULL, Onnz, NULL)); 57 /* The allocation above is approximate so we must set this option to be permissive. 58 * Real code should preallocate exactly. */ 59 PetscCall(MatSetOption(Asp, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE)); 60 61 /* Check zero rows */ 62 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 63 nrows = 0; 64 for (row = rstart; row < rend; row++) { 65 PetscCall(MatGetRow(A, row, &ncols, &cols, &vals)); 66 nnzA += ncols; 67 norm = 0.0; 68 for (j = 0; j < ncols; j++) { 69 val = PetscAbsScalar(vals[j]); 70 if (norm < val) norm = norm; 71 if (val > dtol) { 72 PetscCall(MatSetValues(Asp, 1, &row, 1, &cols[j], &vals[j], INSERT_VALUES)); 73 nnzAsp++; 74 } 75 } 76 if (!norm) nrows++; 77 PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals)); 78 } 79 PetscCall(MatAssemblyBegin(Asp, MAT_FINAL_ASSEMBLY)); 80 PetscCall(MatAssemblyEnd(Asp, MAT_FINAL_ASSEMBLY)); 81 82 percent = (PetscReal)nnzA * 100 / (m * n); 83 PetscCall(PetscPrintf(PETSC_COMM_SELF, " [%d] Matrix A local size %d,%d; nnzA %d, %g percent; No. of zero rows: %d\n", rank, m, n, nnzA, percent, nrows)); 84 percent = (PetscReal)nnzAsp * 100 / (m * n); 85 PetscCall(PetscPrintf(PETSC_COMM_SELF, " [%d] Matrix Asp nnzAsp %d, %g percent\n", rank, nnzAsp, percent)); 86 87 /* investigate matcoloring for Asp */ 88 PetscBool Asp_coloring = PETSC_FALSE; 89 PetscCall(PetscOptionsHasName(NULL, NULL, "-Asp_color", &Asp_coloring)); 90 if (Asp_coloring) { 91 MatColoring mc; 92 ISColoring iscoloring; 93 MatFDColoring matfdcoloring; 94 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Create coloring of Asp...\n")); 95 PetscCall(MatColoringCreate(Asp, &mc)); 96 PetscCall(MatColoringSetType(mc, MATCOLORINGSL)); 97 PetscCall(MatColoringSetFromOptions(mc)); 98 PetscCall(MatColoringApply(mc, &iscoloring)); 99 PetscCall(MatColoringDestroy(&mc)); 100 PetscCall(MatFDColoringCreate(Asp, iscoloring, &matfdcoloring)); 101 PetscCall(MatFDColoringSetFromOptions(matfdcoloring)); 102 PetscCall(MatFDColoringSetUp(Asp, iscoloring, matfdcoloring)); 103 /*PetscCall(MatFDColoringView(matfdcoloring,PETSC_VIEWER_STDOUT_WORLD));*/ 104 PetscCall(ISColoringDestroy(&iscoloring)); 105 PetscCall(MatFDColoringDestroy(&matfdcoloring)); 106 } 107 108 /* Write Asp in binary for study - see ~petsc/src/mat/tests/ex124.c */ 109 PetscBool Asp_write = PETSC_FALSE; 110 PetscCall(PetscOptionsHasName(NULL, NULL, "-Asp_write", &Asp_write)); 111 if (Asp_write) { 112 PetscViewer viewer; 113 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Write Asp into file Asp.dat ...\n")); 114 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "Asp.dat", FILE_MODE_WRITE, &viewer)); 115 PetscCall(MatView(Asp, viewer)); 116 PetscCall(PetscViewerDestroy(&viewer)); 117 } 118 119 PetscCall(MatDestroy(&A)); 120 PetscCall(MatDestroy(&Asp)); 121 PetscCall(PetscFinalize()); 122 return 0; 123 } 124