1 static char help[] = "Test conversion of ScaLAPACK matrices.\n\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** argv)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, 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));
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 double
94
95 testset:
96 output_file: output/empty.out
97 nsize: 6
98
99 test:
100
101 test:
102 suffix: 2
103 args: -mat_type aij
104
105 test:
106 suffix: 3
107 args: -mat_type scalapack
108 TEST*/
109