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*/