Actual source code: ex39.c


  2: static char help[] = "Tests Elemental interface.\n\n";

  4: #include <petscmat.h>

  6: int main(int argc, char **args)
  7: {
  8:   Mat             Cdense, B, C, Ct;
  9:   Vec             d;
 10:   PetscInt        i, j, m = 5, n, nrows, ncols;
 11:   const PetscInt *rows, *cols;
 12:   IS              isrows, iscols;
 13:   PetscScalar    *v;
 14:   PetscMPIInt     rank, size;
 15:   PetscReal       Cnorm;
 16:   PetscBool       flg, mats_view = PETSC_FALSE;

 18:   PetscFunctionBeginUser;
 19:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 20:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 21:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 22:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 23:   n = m;
 24:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 25:   PetscCall(PetscOptionsHasName(NULL, NULL, "-mats_view", &mats_view));

 27:   PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
 28:   PetscCall(MatSetSizes(C, m, n, PETSC_DECIDE, PETSC_DECIDE));
 29:   PetscCall(MatSetType(C, MATELEMENTAL));
 30:   PetscCall(MatSetFromOptions(C));
 31:   PetscCall(MatSetUp(C));
 32:   PetscCall(MatGetOwnershipIS(C, &isrows, &iscols));
 33:   PetscCall(ISGetLocalSize(isrows, &nrows));
 34:   PetscCall(ISGetIndices(isrows, &rows));
 35:   PetscCall(ISGetLocalSize(iscols, &ncols));
 36:   PetscCall(ISGetIndices(iscols, &cols));
 37:   PetscCall(PetscMalloc1(nrows * ncols, &v));
 38: #if defined(PETSC_USE_COMPLEX)
 39:   PetscRandom rand;
 40:   PetscScalar rval;
 41:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
 42:   PetscCall(PetscRandomSetFromOptions(rand));
 43:   for (i = 0; i < nrows; i++) {
 44:     for (j = 0; j < ncols; j++) {
 45:       PetscCall(PetscRandomGetValue(rand, &rval));
 46:       v[i * ncols + j] = rval;
 47:     }
 48:   }
 49:   PetscCall(PetscRandomDestroy(&rand));
 50: #else
 51:   for (i = 0; i < nrows; i++) {
 52:     for (j = 0; j < ncols; j++) v[i * ncols + j] = (PetscReal)(10000 * rank + 100 * rows[i] + cols[j]);
 53:   }
 54: #endif
 55:   PetscCall(MatSetValues(C, nrows, rows, ncols, cols, v, INSERT_VALUES));
 56:   PetscCall(ISRestoreIndices(isrows, &rows));
 57:   PetscCall(ISRestoreIndices(iscols, &cols));
 58:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
 59:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
 60:   PetscCall(ISDestroy(&isrows));
 61:   PetscCall(ISDestroy(&iscols));

 63:   /* Test MatView(), MatDuplicate() and out-of-place MatConvert() */
 64:   PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &B));
 65:   if (mats_view) {
 66:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Duplicated C:\n"));
 67:     PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
 68:   }
 69:   PetscCall(MatDestroy(&B));
 70:   PetscCall(MatConvert(C, MATMPIDENSE, MAT_INITIAL_MATRIX, &Cdense));
 71:   PetscCall(MatMultEqual(C, Cdense, 5, &flg));
 72:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "Cdense != C. MatConvert() fails");

 74:   /* Test MatNorm() */
 75:   PetscCall(MatNorm(C, NORM_1, &Cnorm));

 77:   /* Test MatTranspose(), MatZeroEntries() and MatGetDiagonal() */
 78:   PetscCall(MatTranspose(C, MAT_INITIAL_MATRIX, &Ct));
 79:   PetscCall(MatConjugate(Ct));
 80:   if (mats_view) {
 81:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C's Transpose Conjugate:\n"));
 82:     PetscCall(MatView(Ct, PETSC_VIEWER_STDOUT_WORLD));
 83:   }
 84:   PetscCall(MatZeroEntries(Ct));
 85:   PetscCall(VecCreate(PETSC_COMM_WORLD, &d));
 86:   PetscCall(VecSetSizes(d, m > n ? n : m, PETSC_DECIDE));
 87:   PetscCall(VecSetFromOptions(d));
 88:   PetscCall(MatGetDiagonal(C, d));
 89:   if (mats_view) {
 90:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Diagonal of C:\n"));
 91:     PetscCall(VecView(d, PETSC_VIEWER_STDOUT_WORLD));
 92:   }
 93:   if (m > n) {
 94:     PetscCall(MatDiagonalScale(C, NULL, d));
 95:   } else {
 96:     PetscCall(MatDiagonalScale(C, d, NULL));
 97:   }
 98:   if (mats_view) {
 99:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Diagonal Scaled C:\n"));
100:     PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
101:   }

103:   /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */
104:   PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
105:   PetscCall(MatSetSizes(B, m, n, PETSC_DECIDE, PETSC_DECIDE));
106:   PetscCall(MatSetType(B, MATELEMENTAL));
107:   PetscCall(MatSetFromOptions(B));
108:   PetscCall(MatSetUp(B));
109:   PetscCall(MatGetOwnershipIS(B, &isrows, &iscols));
110:   PetscCall(ISGetLocalSize(isrows, &nrows));
111:   PetscCall(ISGetIndices(isrows, &rows));
112:   PetscCall(ISGetLocalSize(iscols, &ncols));
113:   PetscCall(ISGetIndices(iscols, &cols));
114:   for (i = 0; i < nrows; i++) {
115:     for (j = 0; j < ncols; j++) v[i * ncols + j] = (PetscReal)(1000 * rows[i] + cols[j]);
116:   }
117:   PetscCall(MatSetValues(B, nrows, rows, ncols, cols, v, INSERT_VALUES));
118:   PetscCall(PetscFree(v));
119:   PetscCall(ISRestoreIndices(isrows, &rows));
120:   PetscCall(ISRestoreIndices(iscols, &cols));
121:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
122:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
123:   PetscCall(MatAXPY(B, 2.5, C, SAME_NONZERO_PATTERN));
124:   PetscCall(MatAYPX(B, 3.75, C, SAME_NONZERO_PATTERN));
125:   PetscCall(MatConvert(B, MATDENSE, MAT_INPLACE_MATRIX, &B));
126:   if (mats_view) {
127:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B after MatAXPY and MatAYPX:\n"));
128:     PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
129:   }
130:   PetscCall(ISDestroy(&isrows));
131:   PetscCall(ISDestroy(&iscols));
132:   PetscCall(MatDestroy(&B));

134:   /* Test MatMatTransposeMult(): B = C*C^T */
135:   PetscCall(MatMatTransposeMult(C, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &B));
136:   PetscCall(MatScale(C, 2.0));
137:   PetscCall(MatMatTransposeMult(C, C, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B));
138:   PetscCall(MatMatTransposeMultEqual(C, C, B, 10, &flg));
139:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MatMatTransposeMult: B != C*B^T");

141:   if (mats_view) {
142:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C MatMatTransposeMult C:\n"));
143:     PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
144:   }

146:   PetscCall(MatDestroy(&Cdense));
147:   PetscCall(PetscFree(v));
148:   PetscCall(MatDestroy(&B));
149:   PetscCall(MatDestroy(&C));
150:   PetscCall(MatDestroy(&Ct));
151:   PetscCall(VecDestroy(&d));
152:   PetscCall(PetscFinalize());
153:   return 0;
154: }

156: /*TEST

158:    test:
159:       nsize: 2
160:       args: -m 3 -n 2
161:       requires: elemental

163:    test:
164:       suffix: 2
165:       nsize: 6
166:       args: -m 2 -n 3
167:       requires: elemental

169:    test:
170:       suffix: 3
171:       nsize: 1
172:       args: -m 2 -n 3
173:       requires: elemental
174:       output_file: output/ex39_1.out

176: TEST*/