Actual source code: ex243.c

  1: static char help[] = "Test conversion of ScaLAPACK matrices.\n\n";

  3: #include <petscmat.h>

  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;

 16:   PetscFunctionBeginUser;
 17:   PetscCall(PetscInitialize(&argc, &argv, NULL, 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));

 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));

 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));

 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));

 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:   }

 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));

 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:   }

 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: }

 90: /*TEST

 92:    build:
 93:       requires: scalapack

 95:    test:
 96:       requires: !single # garbage prints in single precision from sgemr2d
 97:       nsize: 6

 99:    test:
100:       requires: !single # garbage prints in single precision from sgemr2d
101:       suffix: 2
102:       nsize: 6
103:       args: -mat_type aij
104:       output_file: output/ex243_1.out

106:    test:
107:       requires: !single # garbage prints in single precision from sgemr2d
108:       suffix: 3
109:       nsize: 6
110:       args: -mat_type scalapack
111:       output_file: output/ex243_1.out

113: TEST*/