Actual source code: ex242.c
1: static char help[] = "Tests ScaLAPACK interface.\n\n";
3: #include <petscmat.h>
5: int main(int argc, char **args)
6: {
7: Mat Cdense, Caij, B, C, Ct, Asub;
8: Vec d;
9: PetscInt i, j, M = 5, N, mb = 2, nb, nrows, ncols, mloc, nloc;
10: const PetscInt *rows, *cols;
11: IS isrows, iscols;
12: PetscScalar *v;
13: PetscMPIInt rank, color;
14: PetscReal Cnorm;
15: PetscBool flg, mats_view = PETSC_FALSE;
16: MPI_Comm subcomm;
18: PetscFunctionBeginUser;
19: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
20: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
21: PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
22: N = M;
23: PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
24: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mb", &mb, NULL));
25: nb = mb;
26: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nb", &nb, NULL));
27: PetscCall(PetscOptionsHasName(NULL, NULL, "-mats_view", &mats_view));
29: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
30: PetscCall(MatSetType(C, MATSCALAPACK));
31: mloc = PETSC_DECIDE;
32: PetscCall(PetscSplitOwnershipEqual(PETSC_COMM_WORLD, &mloc, &M));
33: nloc = PETSC_DECIDE;
34: PetscCall(PetscSplitOwnershipEqual(PETSC_COMM_WORLD, &nloc, &N));
35: PetscCall(MatSetSizes(C, mloc, nloc, M, N));
36: PetscCall(MatScaLAPACKSetBlockSizes(C, mb, nb));
37: PetscCall(MatSetFromOptions(C));
38: PetscCall(MatSetUp(C));
39: /*PetscCall(MatCreateScaLAPACK(PETSC_COMM_WORLD,mb,nb,M,N,0,0,&C)); */
41: PetscCall(MatGetOwnershipIS(C, &isrows, &iscols));
42: PetscCall(ISGetLocalSize(isrows, &nrows));
43: PetscCall(ISGetIndices(isrows, &rows));
44: PetscCall(ISGetLocalSize(iscols, &ncols));
45: PetscCall(ISGetIndices(iscols, &cols));
46: PetscCall(PetscMalloc1(nrows * ncols, &v));
47: for (i = 0; i < nrows; i++) {
48: for (j = 0; j < ncols; j++) v[i * ncols + j] = (PetscReal)(rows[i] + 1 + (cols[j] + 1) * 0.01);
49: }
50: PetscCall(MatSetValues(C, nrows, rows, ncols, cols, v, INSERT_VALUES));
51: PetscCall(PetscFree(v));
52: PetscCall(ISRestoreIndices(isrows, &rows));
53: PetscCall(ISRestoreIndices(iscols, &cols));
54: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
55: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
56: PetscCall(ISDestroy(&isrows));
57: PetscCall(ISDestroy(&iscols));
59: /* Test MatView(), MatDuplicate() and out-of-place MatConvert() */
60: PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &B));
61: if (mats_view) {
62: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Duplicated C:\n"));
63: PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
64: }
65: PetscCall(MatDestroy(&B));
66: PetscCall(MatConvert(C, MATDENSE, MAT_INITIAL_MATRIX, &Cdense));
67: PetscCall(MatMultEqual(C, Cdense, 5, &flg));
68: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Check fails: Cdense != C");
70: /* Test MatNorm() */
71: PetscCall(MatNorm(C, NORM_1, &Cnorm));
73: /* Test MatTranspose(), MatZeroEntries() and MatGetDiagonal() */
74: PetscCall(MatTranspose(C, MAT_INITIAL_MATRIX, &Ct));
75: PetscCall(MatConjugate(Ct));
76: if (mats_view) {
77: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C's Transpose Conjugate:\n"));
78: PetscCall(MatView(Ct, PETSC_VIEWER_STDOUT_WORLD));
79: }
80: PetscCall(MatZeroEntries(Ct));
81: if (M > N) PetscCall(MatCreateVecs(C, &d, NULL));
82: else PetscCall(MatCreateVecs(C, NULL, &d));
83: PetscCall(MatGetDiagonal(C, d));
84: if (mats_view) {
85: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Diagonal of C:\n"));
86: PetscCall(VecView(d, PETSC_VIEWER_STDOUT_WORLD));
87: }
88: if (M > N) {
89: PetscCall(MatDiagonalScale(C, NULL, d));
90: } else {
91: PetscCall(MatDiagonalScale(C, d, NULL));
92: }
93: if (mats_view) {
94: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Diagonal Scaled C:\n"));
95: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
96: }
98: /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */
99: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
100: PetscCall(MatSetType(B, MATSCALAPACK));
101: PetscCall(MatSetSizes(B, mloc, nloc, PETSC_DECIDE, PETSC_DECIDE));
102: PetscCall(MatScaLAPACKSetBlockSizes(B, mb, nb));
103: PetscCall(MatSetFromOptions(B));
104: PetscCall(MatSetUp(B));
105: /* PetscCall(MatCreateScaLAPACK(PETSC_COMM_WORLD,mb,nb,M,N,0,0,&B)); */
106: PetscCall(MatGetOwnershipIS(B, &isrows, &iscols));
107: PetscCall(ISGetLocalSize(isrows, &nrows));
108: PetscCall(ISGetIndices(isrows, &rows));
109: PetscCall(ISGetLocalSize(iscols, &ncols));
110: PetscCall(ISGetIndices(iscols, &cols));
111: PetscCall(PetscMalloc1(nrows * ncols, &v));
112: for (i = 0; i < nrows; i++) {
113: for (j = 0; j < ncols; j++) v[i * ncols + j] = (PetscReal)(1000 * rows[i] + cols[j]);
114: }
115: PetscCall(MatSetValues(B, nrows, rows, ncols, cols, v, INSERT_VALUES));
116: PetscCall(PetscFree(v));
117: PetscCall(ISRestoreIndices(isrows, &rows));
118: PetscCall(ISRestoreIndices(iscols, &cols));
119: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
120: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
121: if (mats_view) {
122: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B:\n"));
123: PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
124: }
125: PetscCall(MatAXPY(B, 2.5, C, SAME_NONZERO_PATTERN));
126: PetscCall(MatAYPX(B, 3.75, C, SAME_NONZERO_PATTERN));
127: PetscCall(MatConvert(B, MATDENSE, MAT_INPLACE_MATRIX, &B));
128: if (mats_view) {
129: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B after MatAXPY and MatAYPX:\n"));
130: PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
131: }
132: PetscCall(ISDestroy(&isrows));
133: PetscCall(ISDestroy(&iscols));
134: PetscCall(MatDestroy(&B));
136: /* Test MatMatTransposeMult(): B = C*C^T */
137: PetscCall(MatMatTransposeMult(C, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &B));
138: PetscCall(MatScale(C, 2.0));
139: PetscCall(MatMatTransposeMult(C, C, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B));
140: PetscCall(MatMatTransposeMultEqual(C, C, B, 10, &flg));
141: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Check fails: B != C*C^T");
143: if (mats_view) {
144: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C MatMatTransposeMult C:\n"));
145: PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
146: }
148: /* Test MatMult() */
149: PetscCall(MatComputeOperator(C, MATAIJ, &Caij));
150: PetscCall(MatMultEqual(C, Caij, 5, &flg));
151: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultEqual() fails");
152: PetscCall(MatMultTransposeEqual(C, Caij, 5, &flg));
153: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultTransposeEqual() fails");
155: /* Test MatMultAdd() and MatMultTransposeAddEqual() */
156: PetscCall(MatMultAddEqual(C, Caij, 5, &flg));
157: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultAddEqual() fails");
158: PetscCall(MatMultTransposeAddEqual(C, Caij, 5, &flg));
159: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultTransposeAddEqual() fails");
161: /* Test MatMatMult() */
162: PetscCall(PetscOptionsHasName(NULL, NULL, "-test_matmatmult", &flg));
163: if (flg) {
164: Mat CC, CCaij;
165: PetscCall(MatMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CC));
166: PetscCall(MatMatMult(Caij, Caij, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CCaij));
167: PetscCall(MatMultEqual(CC, CCaij, 5, &flg));
168: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "CC != CCaij. MatMatMult() fails");
169: PetscCall(MatDestroy(&CCaij));
170: PetscCall(MatDestroy(&CC));
171: }
173: /* Test MatCreate() on subcomm */
174: color = rank % 2;
175: PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, color, 0, &subcomm));
176: if (color == 0) {
177: PetscCall(MatCreate(subcomm, &Asub));
178: PetscCall(MatSetType(Asub, MATSCALAPACK));
179: mloc = PETSC_DECIDE;
180: PetscCall(PetscSplitOwnershipEqual(subcomm, &mloc, &M));
181: nloc = PETSC_DECIDE;
182: PetscCall(PetscSplitOwnershipEqual(subcomm, &nloc, &N));
183: PetscCall(MatSetSizes(Asub, mloc, nloc, M, N));
184: PetscCall(MatScaLAPACKSetBlockSizes(Asub, mb, nb));
185: PetscCall(MatSetFromOptions(Asub));
186: PetscCall(MatSetUp(Asub));
187: PetscCall(MatDestroy(&Asub));
188: }
190: PetscCall(MatDestroy(&Cdense));
191: PetscCall(MatDestroy(&Caij));
192: PetscCall(MatDestroy(&B));
193: PetscCall(MatDestroy(&C));
194: PetscCall(MatDestroy(&Ct));
195: PetscCall(VecDestroy(&d));
196: PetscCallMPI(MPI_Comm_free(&subcomm));
197: PetscCall(PetscFinalize());
198: return 0;
199: }
201: /*TEST
203: build:
204: requires: scalapack
206: test:
207: nsize: 2
208: args: -mb 5 -nb 5 -M 12 -N 10
209: requires: scalapack
211: test:
212: suffix: 2
213: nsize: 6
214: args: -mb 8 -nb 6 -M 20 -N 50
215: requires: scalapack
216: output_file: output/ex242_1.out
218: test:
219: suffix: 3
220: nsize: 3
221: args: -mb 2 -nb 2 -M 20 -N 20 -test_matmatmult
222: requires: scalapack
223: output_file: output/ex242_1.out
225: TEST*/