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, NULL, 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) PetscCall(MatDiagonalScale(C, NULL, d));
89: else PetscCall(MatDiagonalScale(C, d, NULL));
90: if (mats_view) {
91: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Diagonal Scaled C:\n"));
92: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
93: }
95: /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */
96: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
97: PetscCall(MatSetType(B, MATSCALAPACK));
98: PetscCall(MatSetSizes(B, mloc, nloc, PETSC_DECIDE, PETSC_DECIDE));
99: PetscCall(MatScaLAPACKSetBlockSizes(B, mb, nb));
100: PetscCall(MatSetFromOptions(B));
101: PetscCall(MatSetUp(B));
102: /* PetscCall(MatCreateScaLAPACK(PETSC_COMM_WORLD,mb,nb,M,N,0,0,&B)); */
103: PetscCall(MatGetOwnershipIS(B, &isrows, &iscols));
104: PetscCall(ISGetLocalSize(isrows, &nrows));
105: PetscCall(ISGetIndices(isrows, &rows));
106: PetscCall(ISGetLocalSize(iscols, &ncols));
107: PetscCall(ISGetIndices(iscols, &cols));
108: PetscCall(PetscMalloc1(nrows * ncols, &v));
109: for (i = 0; i < nrows; i++) {
110: for (j = 0; j < ncols; j++) v[i * ncols + j] = (PetscReal)(1000 * rows[i] + cols[j]);
111: }
112: PetscCall(MatSetValues(B, nrows, rows, ncols, cols, v, INSERT_VALUES));
113: PetscCall(PetscFree(v));
114: PetscCall(ISRestoreIndices(isrows, &rows));
115: PetscCall(ISRestoreIndices(iscols, &cols));
116: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
117: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
118: if (mats_view) {
119: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B:\n"));
120: PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
121: }
122: PetscCall(MatAXPY(B, 2.5, C, SAME_NONZERO_PATTERN));
123: PetscCall(MatAYPX(B, 3.75, C, SAME_NONZERO_PATTERN));
124: PetscCall(MatConvert(B, MATDENSE, MAT_INPLACE_MATRIX, &B));
125: if (mats_view) {
126: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B after MatAXPY and MatAYPX:\n"));
127: PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
128: }
129: PetscCall(ISDestroy(&isrows));
130: PetscCall(ISDestroy(&iscols));
131: PetscCall(MatDestroy(&B));
133: /* Test MatMatTransposeMult(): B = C*C^T */
134: PetscCall(MatMatTransposeMult(C, C, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &B));
135: PetscCall(MatScale(C, 2.0));
136: PetscCall(MatMatTransposeMult(C, C, MAT_REUSE_MATRIX, PETSC_DETERMINE, &B));
137: PetscCall(MatMatTransposeMultEqual(C, C, B, 10, &flg));
138: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Check fails: B != C*C^T");
140: if (mats_view) {
141: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C MatMatTransposeMult C:\n"));
142: PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
143: }
144: PetscCall(MatDestroy(&B));
146: /* Test MatTransposeMatMult(): B = C^T*C */
147: PetscCall(MatTransposeMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &B));
148: PetscCall(MatScale(C, 2.0));
149: PetscCall(MatTransposeMatMult(C, C, MAT_REUSE_MATRIX, PETSC_DETERMINE, &B));
150: PetscCall(MatTransposeMatMultEqual(C, C, B, 10, &flg));
151: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Check fails: B != C^T*C");
153: if (mats_view) {
154: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C MatTransposeMatMult C:\n"));
155: PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
156: }
158: /* Test MatMult() */
159: PetscCall(MatComputeOperator(C, MATAIJ, &Caij));
160: PetscCall(MatMultEqual(C, Caij, 5, &flg));
161: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultEqual() fails");
162: PetscCall(MatMultTransposeEqual(C, Caij, 5, &flg));
163: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultTransposeEqual() fails");
165: /* Test MatMultAdd() and MatMultTransposeAddEqual() */
166: PetscCall(MatMultAddEqual(C, Caij, 5, &flg));
167: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultAddEqual() fails");
168: PetscCall(MatMultTransposeAddEqual(C, Caij, 5, &flg));
169: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultTransposeAddEqual() fails");
171: /* Test MatMatMult() */
172: PetscCall(PetscOptionsHasName(NULL, NULL, "-test_matmatmult", &flg));
173: if (flg) {
174: Mat CC, CCaij;
175: PetscCall(MatMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &CC));
176: PetscCall(MatMatMult(Caij, Caij, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &CCaij));
177: PetscCall(MatMultEqual(CC, CCaij, 5, &flg));
178: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "CC != CCaij. MatMatMult() fails");
179: PetscCall(MatDestroy(&CCaij));
180: PetscCall(MatDestroy(&CC));
181: }
183: /* Test MatCreate() on subcomm */
184: color = rank % 2;
185: PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, color, 0, &subcomm));
186: if (color == 0) {
187: PetscCall(MatCreate(subcomm, &Asub));
188: PetscCall(MatSetType(Asub, MATSCALAPACK));
189: mloc = PETSC_DECIDE;
190: PetscCall(PetscSplitOwnershipEqual(subcomm, &mloc, &M));
191: nloc = PETSC_DECIDE;
192: PetscCall(PetscSplitOwnershipEqual(subcomm, &nloc, &N));
193: PetscCall(MatSetSizes(Asub, mloc, nloc, M, N));
194: PetscCall(MatScaLAPACKSetBlockSizes(Asub, mb, nb));
195: PetscCall(MatSetFromOptions(Asub));
196: PetscCall(MatSetUp(Asub));
197: PetscCall(MatDestroy(&Asub));
198: }
200: PetscCall(MatDestroy(&Cdense));
201: PetscCall(MatDestroy(&Caij));
202: PetscCall(MatDestroy(&B));
203: PetscCall(MatDestroy(&C));
204: PetscCall(MatDestroy(&Ct));
205: PetscCall(VecDestroy(&d));
206: PetscCallMPI(MPI_Comm_free(&subcomm));
207: PetscCall(PetscFinalize());
208: return 0;
209: }
211: /*TEST
213: build:
214: requires: scalapack double
216: testset:
217: output_file: output/empty.out
219: test:
220: nsize: 2
221: args: -mb 5 -nb 5 -M 12 -N 10
223: test:
224: suffix: 2
225: nsize: 6
226: args: -mb 8 -nb 6 -M 20 -N 50
228: test:
229: suffix: 3
230: nsize: 3
231: args: -mb 2 -nb 2 -M 20 -N 20 -test_matmatmult
233: TEST*/