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