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