1 2 static char help[] = "Tests Elemental interface.\n\n"; 3 4 #include <petscmat.h> 5 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; 17 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)); 26 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)); 62 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"); 73 74 /* Test MatNorm() */ 75 PetscCall(MatNorm(C, NORM_1, &Cnorm)); 76 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 } 102 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)); 133 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"); 140 141 if (mats_view) { 142 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C MatMatTransposeMult C:\n")); 143 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD)); 144 } 145 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 } 155 156 /*TEST 157 158 test: 159 nsize: 2 160 args: -m 3 -n 2 161 requires: elemental 162 163 test: 164 suffix: 2 165 nsize: 6 166 args: -m 2 -n 3 167 requires: elemental 168 169 test: 170 suffix: 3 171 nsize: 1 172 args: -m 2 -n 3 173 requires: elemental 174 output_file: output/ex39_1.out 175 176 TEST*/ 177