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