1 static char help[] = "Test MatSetValues() by converting MATDENSE to MATELEMENTAL. \n\ 2 Modified from the code contributed by Yaning Liu @lbl.gov \n\n"; 3 /* 4 Example: 5 mpiexec -n <np> ./ex103 6 mpiexec -n <np> ./ex103 -mat_type elemental -mat_view 7 mpiexec -n <np> ./ex103 -mat_type aij 8 */ 9 10 #include <petscmat.h> 11 12 int main(int argc, char **argv) 13 { 14 Mat A, A_elemental; 15 PetscInt i, j, M = 10, N = 5, nrows, ncols; 16 PetscMPIInt rank, size; 17 IS isrows, iscols; 18 const PetscInt *rows, *cols; 19 PetscScalar *v; 20 MatType type; 21 PetscBool isDense, isAIJ, flg; 22 23 PetscFunctionBeginUser; 24 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 25 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 26 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 27 28 /* Creat a matrix */ 29 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL)); 30 PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL)); 31 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 32 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N)); 33 PetscCall(MatSetType(A, MATDENSE)); 34 PetscCall(MatSetFromOptions(A)); 35 PetscCall(MatSetUp(A)); 36 37 /* Set local matrix entries */ 38 PetscCall(MatGetOwnershipIS(A, &isrows, &iscols)); 39 PetscCall(ISGetLocalSize(isrows, &nrows)); 40 PetscCall(ISGetIndices(isrows, &rows)); 41 PetscCall(ISGetLocalSize(iscols, &ncols)); 42 PetscCall(ISGetIndices(iscols, &cols)); 43 PetscCall(PetscMalloc1(nrows * ncols, &v)); 44 45 for (i = 0; i < nrows; i++) { 46 for (j = 0; j < ncols; j++) { 47 if (size == 1) { 48 v[i * ncols + j] = (PetscScalar)(i + j); 49 } else { 50 v[i * ncols + j] = (PetscScalar)rank + j * 0.1; 51 } 52 } 53 } 54 PetscCall(MatSetValues(A, nrows, rows, ncols, cols, v, INSERT_VALUES)); 55 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 56 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 57 //PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%" PetscInt_FMT "] local nrows %" PetscInt_FMT ", ncols %" PetscInt_FMT "\n",rank,nrows,ncols)); 58 //PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT)); 59 60 /* Test MatSetValues() by converting A to A_elemental */ 61 PetscCall(MatGetType(A, &type)); 62 if (size == 1) { 63 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQDENSE, &isDense)); 64 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &isAIJ)); 65 } else { 66 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIDENSE, &isDense)); 67 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIAIJ, &isAIJ)); 68 } 69 70 if (isDense || isAIJ) { 71 Mat Aexplicit; 72 PetscCall(MatConvert(A, MATELEMENTAL, MAT_INITIAL_MATRIX, &A_elemental)); 73 PetscCall(MatComputeOperator(A_elemental, isAIJ ? MATAIJ : MATDENSE, &Aexplicit)); 74 PetscCall(MatMultEqual(Aexplicit, A_elemental, 5, &flg)); 75 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Aexplicit != A_elemental."); 76 PetscCall(MatDestroy(&Aexplicit)); 77 78 /* Test MAT_REUSE_MATRIX which is only supported for inplace conversion */ 79 PetscCall(MatConvert(A, MATELEMENTAL, MAT_INPLACE_MATRIX, &A)); 80 PetscCall(MatMultEqual(A_elemental, A, 5, &flg)); 81 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "A_elemental != A."); 82 PetscCall(MatDestroy(&A_elemental)); 83 } 84 85 PetscCall(ISRestoreIndices(isrows, &rows)); 86 PetscCall(ISRestoreIndices(iscols, &cols)); 87 PetscCall(ISDestroy(&isrows)); 88 PetscCall(ISDestroy(&iscols)); 89 PetscCall(PetscFree(v)); 90 PetscCall(MatDestroy(&A)); 91 PetscCall(PetscFinalize()); 92 return 0; 93 } 94 95 /*TEST 96 97 build: 98 requires: elemental 99 100 test: 101 nsize: 6 102 103 test: 104 suffix: 2 105 nsize: 6 106 args: -mat_type aij 107 output_file: output/ex103_1.out 108 109 test: 110 suffix: 3 111 nsize: 6 112 args: -mat_type elemental 113 output_file: output/ex103_1.out 114 115 TEST*/ 116