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;

 17:   PetscInitialize(&argc, &argv, (char *)0, help);
 18:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 19:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 20:   PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL);
 21:   PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);

 23:   /* Create a matrix */
 24:   MatCreate(PETSC_COMM_WORLD, &A);
 25:   mloc = PETSC_DECIDE;
 26:   PetscSplitOwnershipEqual(PETSC_COMM_WORLD, &mloc, &M);
 27:   nloc = PETSC_DECIDE;
 28:   PetscSplitOwnershipEqual(PETSC_COMM_WORLD, &nloc, &N);
 29:   MatSetSizes(A, mloc, nloc, M, N);
 30:   MatSetType(A, MATDENSE);
 31:   MatSetFromOptions(A);
 32:   MatSetUp(A);

 34:   /* Set local matrix entries */
 35:   MatGetOwnershipIS(A, &isrows, &iscols);
 36:   ISGetLocalSize(isrows, &nrows);
 37:   ISGetIndices(isrows, &rows);
 38:   ISGetLocalSize(iscols, &ncols);
 39:   ISGetIndices(iscols, &cols);
 40:   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:   MatSetValues(A, nrows, rows, ncols, cols, v, INSERT_VALUES);
 52:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 53:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

 55:   /* Test MatSetValues() by converting A to A_scalapack */
 56:   MatGetType(A, &type);
 57:   if (size == 1) {
 58:     PetscObjectTypeCompare((PetscObject)A, MATSEQDENSE, &isDense);
 59:     PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &isAIJ);
 60:   } else {
 61:     PetscObjectTypeCompare((PetscObject)A, MATMPIDENSE, &isDense);
 62:     PetscObjectTypeCompare((PetscObject)A, MATMPIAIJ, &isAIJ);
 63:   }

 65:   if (isDense || isAIJ) {
 66:     Mat Aexplicit;
 67:     MatConvert(A, MATSCALAPACK, MAT_INITIAL_MATRIX, &A_scalapack);
 68:     MatComputeOperator(A_scalapack, isAIJ ? MATAIJ : MATDENSE, &Aexplicit);
 69:     MatMultEqual(Aexplicit, A_scalapack, 5, &flg);
 71:     MatDestroy(&Aexplicit);

 73:     /* Test MAT_REUSE_MATRIX which is only supported for inplace conversion */
 74:     MatConvert(A, MATSCALAPACK, MAT_INPLACE_MATRIX, &A);
 75:     MatMultEqual(A_scalapack, A, 5, &flg);
 77:     MatDestroy(&A_scalapack);
 78:   }

 80:   ISRestoreIndices(isrows, &rows);
 81:   ISRestoreIndices(iscols, &cols);
 82:   ISDestroy(&isrows);
 83:   ISDestroy(&iscols);
 84:   PetscFree(v);
 85:   MatDestroy(&A);
 86:   PetscFinalize();
 87:   return 0;
 88: }

 90: /*TEST

 92:    build:
 93:       requires: scalapack

 95:    test:
 96:       nsize: 6

 98:    test:
 99:       suffix: 2
100:       nsize: 6
101:       args: -mat_type aij
102:       output_file: output/ex243_1.out

104:    test:
105:       suffix: 3
106:       nsize: 6
107:       args: -mat_type scalapack
108:       output_file: output/ex243_1.out

110: TEST*/