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*/