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;

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

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

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

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

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

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

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

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

146:   MatDestroy(&Cdense);
147:   PetscFree(v);
148:   MatDestroy(&B);
149:   MatDestroy(&C);
150:   MatDestroy(&Ct);
151:   VecDestroy(&d);
152:   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*/