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