Actual source code: ex39.c

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

  3: #include <petscmat.h>

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

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

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

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

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

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

102:   /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */
103:   PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
104:   PetscCall(MatSetSizes(B, m, n, PETSC_DECIDE, PETSC_DECIDE));
105:   PetscCall(MatSetType(B, MATELEMENTAL));
106:   PetscCall(MatSetFromOptions(B));
107:   PetscCall(MatSetUp(B));
108:   PetscCall(MatGetOwnershipIS(B, &isrows, &iscols));
109:   PetscCall(ISGetLocalSize(isrows, &nrows));
110:   PetscCall(ISGetIndices(isrows, &rows));
111:   PetscCall(ISGetLocalSize(iscols, &ncols));
112:   PetscCall(ISGetIndices(iscols, &cols));
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:   PetscCall(MatAXPY(B, 2.5, C, SAME_NONZERO_PATTERN));
123:   PetscCall(MatAYPX(B, 3.75, C, SAME_NONZERO_PATTERN));
124:   PetscCall(MatConvert(B, MATDENSE, MAT_INPLACE_MATRIX, &B));
125:   if (mats_view) {
126:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B after MatAXPY and MatAYPX:\n"));
127:     PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
128:   }
129:   PetscCall(ISDestroy(&isrows));
130:   PetscCall(ISDestroy(&iscols));
131:   PetscCall(MatDestroy(&B));

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

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

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

155: /*TEST

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

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

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

175: TEST*/