1 2 static char help[] = "Tests ScaLAPACK interface.\n\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc, char **args) 7 { 8 Mat Cdense, Caij, B, C, Ct, Asub; 9 Vec d; 10 PetscInt i, j, M = 5, N, mb = 2, nb, nrows, ncols, mloc, nloc; 11 const PetscInt *rows, *cols; 12 IS isrows, iscols; 13 PetscScalar *v; 14 PetscMPIInt rank, color; 15 PetscReal Cnorm; 16 PetscBool flg, mats_view = PETSC_FALSE; 17 MPI_Comm subcomm; 18 19 PetscFunctionBeginUser; 20 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 21 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 22 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL)); 23 N = M; 24 PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL)); 25 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mb", &mb, NULL)); 26 nb = mb; 27 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nb", &nb, NULL)); 28 PetscCall(PetscOptionsHasName(NULL, NULL, "-mats_view", &mats_view)); 29 30 PetscCall(MatCreate(PETSC_COMM_WORLD, &C)); 31 PetscCall(MatSetType(C, MATSCALAPACK)); 32 mloc = PETSC_DECIDE; 33 PetscCall(PetscSplitOwnershipEqual(PETSC_COMM_WORLD, &mloc, &M)); 34 nloc = PETSC_DECIDE; 35 PetscCall(PetscSplitOwnershipEqual(PETSC_COMM_WORLD, &nloc, &N)); 36 PetscCall(MatSetSizes(C, mloc, nloc, M, N)); 37 PetscCall(MatScaLAPACKSetBlockSizes(C, mb, nb)); 38 PetscCall(MatSetFromOptions(C)); 39 PetscCall(MatSetUp(C)); 40 /*PetscCall(MatCreateScaLAPACK(PETSC_COMM_WORLD,mb,nb,M,N,0,0,&C)); */ 41 42 PetscCall(MatGetOwnershipIS(C, &isrows, &iscols)); 43 PetscCall(ISGetLocalSize(isrows, &nrows)); 44 PetscCall(ISGetIndices(isrows, &rows)); 45 PetscCall(ISGetLocalSize(iscols, &ncols)); 46 PetscCall(ISGetIndices(iscols, &cols)); 47 PetscCall(PetscMalloc1(nrows * ncols, &v)); 48 for (i = 0; i < nrows; i++) { 49 for (j = 0; j < ncols; j++) v[i * ncols + j] = (PetscReal)(rows[i] + 1 + (cols[j] + 1) * 0.01); 50 } 51 PetscCall(MatSetValues(C, nrows, rows, ncols, cols, v, INSERT_VALUES)); 52 PetscCall(PetscFree(v)); 53 PetscCall(ISRestoreIndices(isrows, &rows)); 54 PetscCall(ISRestoreIndices(iscols, &cols)); 55 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 56 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 57 PetscCall(ISDestroy(&isrows)); 58 PetscCall(ISDestroy(&iscols)); 59 60 /* Test MatView(), MatDuplicate() and out-of-place MatConvert() */ 61 PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &B)); 62 if (mats_view) { 63 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Duplicated C:\n")); 64 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD)); 65 } 66 PetscCall(MatDestroy(&B)); 67 PetscCall(MatConvert(C, MATDENSE, MAT_INITIAL_MATRIX, &Cdense)); 68 PetscCall(MatMultEqual(C, Cdense, 5, &flg)); 69 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Check fails: Cdense != C"); 70 71 /* Test MatNorm() */ 72 PetscCall(MatNorm(C, NORM_1, &Cnorm)); 73 74 /* Test MatTranspose(), MatZeroEntries() and MatGetDiagonal() */ 75 PetscCall(MatTranspose(C, MAT_INITIAL_MATRIX, &Ct)); 76 PetscCall(MatConjugate(Ct)); 77 if (mats_view) { 78 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C's Transpose Conjugate:\n")); 79 PetscCall(MatView(Ct, PETSC_VIEWER_STDOUT_WORLD)); 80 } 81 PetscCall(MatZeroEntries(Ct)); 82 if (M > N) PetscCall(MatCreateVecs(C, &d, NULL)); 83 else PetscCall(MatCreateVecs(C, NULL, &d)); 84 PetscCall(MatGetDiagonal(C, d)); 85 if (mats_view) { 86 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Diagonal of C:\n")); 87 PetscCall(VecView(d, PETSC_VIEWER_STDOUT_WORLD)); 88 } 89 if (M > N) { 90 PetscCall(MatDiagonalScale(C, NULL, d)); 91 } else { 92 PetscCall(MatDiagonalScale(C, d, NULL)); 93 } 94 if (mats_view) { 95 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Diagonal Scaled C:\n")); 96 PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD)); 97 } 98 99 /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */ 100 PetscCall(MatCreate(PETSC_COMM_WORLD, &B)); 101 PetscCall(MatSetType(B, MATSCALAPACK)); 102 PetscCall(MatSetSizes(B, mloc, nloc, PETSC_DECIDE, PETSC_DECIDE)); 103 PetscCall(MatScaLAPACKSetBlockSizes(B, mb, nb)); 104 PetscCall(MatSetFromOptions(B)); 105 PetscCall(MatSetUp(B)); 106 /* PetscCall(MatCreateScaLAPACK(PETSC_COMM_WORLD,mb,nb,M,N,0,0,&B)); */ 107 PetscCall(MatGetOwnershipIS(B, &isrows, &iscols)); 108 PetscCall(ISGetLocalSize(isrows, &nrows)); 109 PetscCall(ISGetIndices(isrows, &rows)); 110 PetscCall(ISGetLocalSize(iscols, &ncols)); 111 PetscCall(ISGetIndices(iscols, &cols)); 112 PetscCall(PetscMalloc1(nrows * ncols, &v)); 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 if (mats_view) { 123 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B:\n")); 124 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD)); 125 } 126 PetscCall(MatAXPY(B, 2.5, C, SAME_NONZERO_PATTERN)); 127 PetscCall(MatAYPX(B, 3.75, C, SAME_NONZERO_PATTERN)); 128 PetscCall(MatConvert(B, MATDENSE, MAT_INPLACE_MATRIX, &B)); 129 if (mats_view) { 130 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B after MatAXPY and MatAYPX:\n")); 131 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD)); 132 } 133 PetscCall(ISDestroy(&isrows)); 134 PetscCall(ISDestroy(&iscols)); 135 PetscCall(MatDestroy(&B)); 136 137 /* Test MatMatTransposeMult(): B = C*C^T */ 138 PetscCall(MatMatTransposeMult(C, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &B)); 139 PetscCall(MatScale(C, 2.0)); 140 PetscCall(MatMatTransposeMult(C, C, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B)); 141 PetscCall(MatMatTransposeMultEqual(C, C, B, 10, &flg)); 142 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Check fails: B != C*C^T"); 143 144 if (mats_view) { 145 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C MatMatTransposeMult C:\n")); 146 PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD)); 147 } 148 149 /* Test MatMult() */ 150 PetscCall(MatComputeOperator(C, MATAIJ, &Caij)); 151 PetscCall(MatMultEqual(C, Caij, 5, &flg)); 152 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultEqual() fails"); 153 PetscCall(MatMultTransposeEqual(C, Caij, 5, &flg)); 154 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultTransposeEqual() fails"); 155 156 /* Test MatMultAdd() and MatMultTransposeAddEqual() */ 157 PetscCall(MatMultAddEqual(C, Caij, 5, &flg)); 158 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultAddEqual() fails"); 159 PetscCall(MatMultTransposeAddEqual(C, Caij, 5, &flg)); 160 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultTransposeAddEqual() fails"); 161 162 /* Test MatMatMult() */ 163 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_matmatmult", &flg)); 164 if (flg) { 165 Mat CC, CCaij; 166 PetscCall(MatMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CC)); 167 PetscCall(MatMatMult(Caij, Caij, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CCaij)); 168 PetscCall(MatMultEqual(CC, CCaij, 5, &flg)); 169 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_NOTSAMETYPE, "CC != CCaij. MatMatMult() fails"); 170 PetscCall(MatDestroy(&CCaij)); 171 PetscCall(MatDestroy(&CC)); 172 } 173 174 /* Test MatCreate() on subcomm */ 175 color = rank % 2; 176 PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, color, 0, &subcomm)); 177 if (color == 0) { 178 PetscCall(MatCreate(subcomm, &Asub)); 179 PetscCall(MatSetType(Asub, MATSCALAPACK)); 180 mloc = PETSC_DECIDE; 181 PetscCall(PetscSplitOwnershipEqual(subcomm, &mloc, &M)); 182 nloc = PETSC_DECIDE; 183 PetscCall(PetscSplitOwnershipEqual(subcomm, &nloc, &N)); 184 PetscCall(MatSetSizes(Asub, mloc, nloc, M, N)); 185 PetscCall(MatScaLAPACKSetBlockSizes(Asub, mb, nb)); 186 PetscCall(MatSetFromOptions(Asub)); 187 PetscCall(MatSetUp(Asub)); 188 PetscCall(MatDestroy(&Asub)); 189 } 190 191 PetscCall(MatDestroy(&Cdense)); 192 PetscCall(MatDestroy(&Caij)); 193 PetscCall(MatDestroy(&B)); 194 PetscCall(MatDestroy(&C)); 195 PetscCall(MatDestroy(&Ct)); 196 PetscCall(VecDestroy(&d)); 197 PetscCallMPI(MPI_Comm_free(&subcomm)); 198 PetscCall(PetscFinalize()); 199 return 0; 200 } 201 202 /*TEST 203 204 build: 205 requires: scalapack 206 207 test: 208 nsize: 2 209 args: -mb 5 -nb 5 -M 12 -N 10 210 requires: scalapack 211 212 test: 213 suffix: 2 214 nsize: 6 215 args: -mb 8 -nb 6 -M 20 -N 50 216 requires: scalapack 217 output_file: output/ex242_1.out 218 219 test: 220 suffix: 3 221 nsize: 3 222 args: -mb 2 -nb 2 -M 20 -N 20 -test_matmatmult 223 requires: scalapack 224 output_file: output/ex242_1.out 225 226 TEST*/ 227