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;
20: PetscInitialize(&argc, &args, (char *)0, help);
21: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
22: PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL);
23: N = M;
24: PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);
25: PetscOptionsGetInt(NULL, NULL, "-mb", &mb, NULL);
26: nb = mb;
27: PetscOptionsGetInt(NULL, NULL, "-nb", &nb, NULL);
28: PetscOptionsHasName(NULL, NULL, "-mats_view", &mats_view);
30: MatCreate(PETSC_COMM_WORLD, &C);
31: MatSetType(C, MATSCALAPACK);
32: mloc = PETSC_DECIDE;
33: PetscSplitOwnershipEqual(PETSC_COMM_WORLD, &mloc, &M);
34: nloc = PETSC_DECIDE;
35: PetscSplitOwnershipEqual(PETSC_COMM_WORLD, &nloc, &N);
36: MatSetSizes(C, mloc, nloc, M, N);
37: MatScaLAPACKSetBlockSizes(C, mb, nb);
38: MatSetFromOptions(C);
39: MatSetUp(C);
40: /*MatCreateScaLAPACK(PETSC_COMM_WORLD,mb,nb,M,N,0,0,&C); */
42: MatGetOwnershipIS(C, &isrows, &iscols);
43: ISGetLocalSize(isrows, &nrows);
44: ISGetIndices(isrows, &rows);
45: ISGetLocalSize(iscols, &ncols);
46: ISGetIndices(iscols, &cols);
47: 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: MatSetValues(C, nrows, rows, ncols, cols, v, INSERT_VALUES);
52: PetscFree(v);
53: ISRestoreIndices(isrows, &rows);
54: ISRestoreIndices(iscols, &cols);
55: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
56: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
57: ISDestroy(&isrows);
58: ISDestroy(&iscols);
60: /* Test MatView(), MatDuplicate() and out-of-place MatConvert() */
61: MatDuplicate(C, MAT_COPY_VALUES, &B);
62: if (mats_view) {
63: PetscPrintf(PETSC_COMM_WORLD, "Duplicated C:\n");
64: MatView(B, PETSC_VIEWER_STDOUT_WORLD);
65: }
66: MatDestroy(&B);
67: MatConvert(C, MATDENSE, MAT_INITIAL_MATRIX, &Cdense);
68: MatMultEqual(C, Cdense, 5, &flg);
71: /* Test MatNorm() */
72: MatNorm(C, NORM_1, &Cnorm);
74: /* Test MatTranspose(), MatZeroEntries() and MatGetDiagonal() */
75: MatTranspose(C, MAT_INITIAL_MATRIX, &Ct);
76: MatConjugate(Ct);
77: if (mats_view) {
78: PetscPrintf(PETSC_COMM_WORLD, "C's Transpose Conjugate:\n");
79: MatView(Ct, PETSC_VIEWER_STDOUT_WORLD);
80: }
81: MatZeroEntries(Ct);
82: if (M > N) MatCreateVecs(C, &d, NULL);
83: else MatCreateVecs(C, NULL, &d);
84: MatGetDiagonal(C, d);
85: if (mats_view) {
86: PetscPrintf(PETSC_COMM_WORLD, "Diagonal of C:\n");
87: VecView(d, PETSC_VIEWER_STDOUT_WORLD);
88: }
89: if (M > N) {
90: MatDiagonalScale(C, NULL, d);
91: } else {
92: MatDiagonalScale(C, d, NULL);
93: }
94: if (mats_view) {
95: PetscPrintf(PETSC_COMM_WORLD, "Diagonal Scaled C:\n");
96: MatView(C, PETSC_VIEWER_STDOUT_WORLD);
97: }
99: /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */
100: MatCreate(PETSC_COMM_WORLD, &B);
101: MatSetType(B, MATSCALAPACK);
102: MatSetSizes(B, mloc, nloc, PETSC_DECIDE, PETSC_DECIDE);
103: MatScaLAPACKSetBlockSizes(B, mb, nb);
104: MatSetFromOptions(B);
105: MatSetUp(B);
106: /* MatCreateScaLAPACK(PETSC_COMM_WORLD,mb,nb,M,N,0,0,&B); */
107: MatGetOwnershipIS(B, &isrows, &iscols);
108: ISGetLocalSize(isrows, &nrows);
109: ISGetIndices(isrows, &rows);
110: ISGetLocalSize(iscols, &ncols);
111: ISGetIndices(iscols, &cols);
112: 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: MatSetValues(B, nrows, rows, ncols, cols, v, INSERT_VALUES);
117: PetscFree(v);
118: ISRestoreIndices(isrows, &rows);
119: ISRestoreIndices(iscols, &cols);
120: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
121: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
122: if (mats_view) {
123: PetscPrintf(PETSC_COMM_WORLD, "B:\n");
124: MatView(B, PETSC_VIEWER_STDOUT_WORLD);
125: }
126: MatAXPY(B, 2.5, C, SAME_NONZERO_PATTERN);
127: MatAYPX(B, 3.75, C, SAME_NONZERO_PATTERN);
128: MatConvert(B, MATDENSE, MAT_INPLACE_MATRIX, &B);
129: if (mats_view) {
130: PetscPrintf(PETSC_COMM_WORLD, "B after MatAXPY and MatAYPX:\n");
131: MatView(B, PETSC_VIEWER_STDOUT_WORLD);
132: }
133: ISDestroy(&isrows);
134: ISDestroy(&iscols);
135: MatDestroy(&B);
137: /* Test MatMatTransposeMult(): B = C*C^T */
138: MatMatTransposeMult(C, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &B);
139: MatScale(C, 2.0);
140: MatMatTransposeMult(C, C, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B);
141: MatMatTransposeMultEqual(C, C, B, 10, &flg);
144: if (mats_view) {
145: PetscPrintf(PETSC_COMM_WORLD, "C MatMatTransposeMult C:\n");
146: MatView(B, PETSC_VIEWER_STDOUT_WORLD);
147: }
149: /* Test MatMult() */
150: MatComputeOperator(C, MATAIJ, &Caij);
151: MatMultEqual(C, Caij, 5, &flg);
153: MatMultTransposeEqual(C, Caij, 5, &flg);
156: /* Test MatMultAdd() and MatMultTransposeAddEqual() */
157: MatMultAddEqual(C, Caij, 5, &flg);
159: MatMultTransposeAddEqual(C, Caij, 5, &flg);
162: /* Test MatMatMult() */
163: PetscOptionsHasName(NULL, NULL, "-test_matmatmult", &flg);
164: if (flg) {
165: Mat CC, CCaij;
166: MatMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CC);
167: MatMatMult(Caij, Caij, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CCaij);
168: MatMultEqual(CC, CCaij, 5, &flg);
170: MatDestroy(&CCaij);
171: MatDestroy(&CC);
172: }
174: /* Test MatCreate() on subcomm */
175: color = rank % 2;
176: MPI_Comm_split(PETSC_COMM_WORLD, color, 0, &subcomm);
177: if (color == 0) {
178: MatCreate(subcomm, &Asub);
179: MatSetType(Asub, MATSCALAPACK);
180: mloc = PETSC_DECIDE;
181: PetscSplitOwnershipEqual(subcomm, &mloc, &M);
182: nloc = PETSC_DECIDE;
183: PetscSplitOwnershipEqual(subcomm, &nloc, &N);
184: MatSetSizes(Asub, mloc, nloc, M, N);
185: MatScaLAPACKSetBlockSizes(Asub, mb, nb);
186: MatSetFromOptions(Asub);
187: MatSetUp(Asub);
188: MatDestroy(&Asub);
189: }
191: MatDestroy(&Cdense);
192: MatDestroy(&Caij);
193: MatDestroy(&B);
194: MatDestroy(&C);
195: MatDestroy(&Ct);
196: VecDestroy(&d);
197: MPI_Comm_free(&subcomm);
198: 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*/