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