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