1 #include <petsc/private/petscscalapack.h> /*I "petscmat.h" I*/ 2 3 const char ScaLAPACKCitation[] = "@BOOK{scalapack-user-guide,\n" 4 " AUTHOR = {L. S. Blackford and J. Choi and A. Cleary and E. D'Azevedo and\n" 5 " J. Demmel and I. Dhillon and J. Dongarra and S. Hammarling and\n" 6 " G. Henry and A. Petitet and K. Stanley and D. Walker and R. C. Whaley},\n" 7 " TITLE = {Sca{LAPACK} Users' Guide},\n" 8 " PUBLISHER = {SIAM},\n" 9 " ADDRESS = {Philadelphia, PA},\n" 10 " YEAR = 1997\n" 11 "}\n"; 12 static PetscBool ScaLAPACKCite = PETSC_FALSE; 13 14 #define DEFAULT_BLOCKSIZE 64 15 16 /* 17 The variable Petsc_ScaLAPACK_keyval is used to indicate an MPI attribute that 18 is attached to a communicator, in this case the attribute is a Mat_ScaLAPACK_Grid 19 */ 20 static PetscMPIInt Petsc_ScaLAPACK_keyval = MPI_KEYVAL_INVALID; 21 22 static PetscErrorCode Petsc_ScaLAPACK_keyval_free(void) 23 { 24 PetscFunctionBegin; 25 PetscCall(PetscInfo(NULL, "Freeing Petsc_ScaLAPACK_keyval\n")); 26 PetscCallMPI(MPI_Comm_free_keyval(&Petsc_ScaLAPACK_keyval)); 27 PetscFunctionReturn(PETSC_SUCCESS); 28 } 29 30 static PetscErrorCode MatView_ScaLAPACK(Mat A, PetscViewer viewer) 31 { 32 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 33 PetscBool isascii; 34 PetscViewerFormat format; 35 Mat Adense; 36 37 PetscFunctionBegin; 38 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 39 if (isascii) { 40 PetscCall(PetscViewerGetFormat(viewer, &format)); 41 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 42 PetscCall(PetscViewerASCIIPrintf(viewer, "block sizes: %d,%d\n", (int)a->mb, (int)a->nb)); 43 PetscCall(PetscViewerASCIIPrintf(viewer, "grid height=%d, grid width=%d\n", (int)a->grid->nprow, (int)a->grid->npcol)); 44 PetscCall(PetscViewerASCIIPrintf(viewer, "coordinates of process owning first row and column: (%d,%d)\n", (int)a->rsrc, (int)a->csrc)); 45 PetscCall(PetscViewerASCIIPrintf(viewer, "dimension of largest local matrix: %d x %d\n", (int)a->locr, (int)a->locc)); 46 PetscFunctionReturn(PETSC_SUCCESS); 47 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 48 PetscFunctionReturn(PETSC_SUCCESS); 49 } 50 } 51 /* convert to dense format and call MatView() */ 52 PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &Adense)); 53 PetscCall(MatView(Adense, viewer)); 54 PetscCall(MatDestroy(&Adense)); 55 PetscFunctionReturn(PETSC_SUCCESS); 56 } 57 58 static PetscErrorCode MatGetInfo_ScaLAPACK(Mat A, MatInfoType flag, MatInfo *info) 59 { 60 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 61 PetscLogDouble isend[2], irecv[2]; 62 63 PetscFunctionBegin; 64 info->block_size = 1.0; 65 66 isend[0] = a->lld * a->locc; /* locally allocated */ 67 isend[1] = a->locr * a->locc; /* used submatrix */ 68 if (flag == MAT_LOCAL || flag == MAT_GLOBAL_MAX) { 69 info->nz_allocated = isend[0]; 70 info->nz_used = isend[1]; 71 } else if (flag == MAT_GLOBAL_MAX) { 72 PetscCallMPI(MPIU_Allreduce(isend, irecv, 2, MPIU_PETSCLOGDOUBLE, MPI_MAX, PetscObjectComm((PetscObject)A))); 73 info->nz_allocated = irecv[0]; 74 info->nz_used = irecv[1]; 75 } else if (flag == MAT_GLOBAL_SUM) { 76 PetscCallMPI(MPIU_Allreduce(isend, irecv, 2, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)A))); 77 info->nz_allocated = irecv[0]; 78 info->nz_used = irecv[1]; 79 } 80 81 info->nz_unneeded = 0; 82 info->assemblies = A->num_ass; 83 info->mallocs = 0; 84 info->memory = 0; /* REVIEW ME */ 85 info->fill_ratio_given = 0; 86 info->fill_ratio_needed = 0; 87 info->factor_mallocs = 0; 88 PetscFunctionReturn(PETSC_SUCCESS); 89 } 90 91 static PetscErrorCode MatSetOption_ScaLAPACK(Mat A, MatOption op, PetscBool flg) 92 { 93 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 94 95 PetscFunctionBegin; 96 switch (op) { 97 case MAT_NEW_NONZERO_LOCATIONS: 98 case MAT_NEW_NONZERO_LOCATION_ERR: 99 case MAT_NEW_NONZERO_ALLOCATION_ERR: 100 case MAT_SYMMETRIC: 101 case MAT_SORTED_FULL: 102 case MAT_HERMITIAN: 103 break; 104 case MAT_ROW_ORIENTED: 105 a->roworiented = flg; 106 break; 107 default: 108 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported option %s", MatOptions[op]); 109 } 110 PetscFunctionReturn(PETSC_SUCCESS); 111 } 112 113 static PetscErrorCode MatSetValues_ScaLAPACK(Mat A, PetscInt nr, const PetscInt *rows, PetscInt nc, const PetscInt *cols, const PetscScalar *vals, InsertMode imode) 114 { 115 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 116 PetscInt i, j; 117 PetscBLASInt gridx, gcidx, lridx, lcidx, rsrc, csrc; 118 PetscBool roworiented = a->roworiented; 119 120 PetscFunctionBegin; 121 PetscCheck(imode == INSERT_VALUES || imode == ADD_VALUES, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for InsertMode %d", (int)imode); 122 for (i = 0; i < nr; i++) { 123 if (rows[i] < 0) continue; 124 PetscCall(PetscBLASIntCast(rows[i] + 1, &gridx)); 125 for (j = 0; j < nc; j++) { 126 if (cols[j] < 0) continue; 127 PetscCall(PetscBLASIntCast(cols[j] + 1, &gcidx)); 128 PetscCallBLAS("SCALAPACKinfog2l", SCALAPACKinfog2l_(&gridx, &gcidx, a->desc, &a->grid->nprow, &a->grid->npcol, &a->grid->myrow, &a->grid->mycol, &lridx, &lcidx, &rsrc, &csrc)); 129 if (rsrc == a->grid->myrow && csrc == a->grid->mycol) { 130 if (roworiented) { 131 switch (imode) { 132 case INSERT_VALUES: 133 a->loc[lridx - 1 + (lcidx - 1) * a->lld] = vals[i * nc + j]; 134 break; 135 default: 136 a->loc[lridx - 1 + (lcidx - 1) * a->lld] += vals[i * nc + j]; 137 break; 138 } 139 } else { 140 switch (imode) { 141 case INSERT_VALUES: 142 a->loc[lridx - 1 + (lcidx - 1) * a->lld] = vals[i + j * nr]; 143 break; 144 default: 145 a->loc[lridx - 1 + (lcidx - 1) * a->lld] += vals[i + j * nr]; 146 break; 147 } 148 } 149 } else { 150 PetscCheck(!A->nooffprocentries, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Setting off process entry even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set"); 151 A->assembled = PETSC_FALSE; 152 PetscCall(MatStashValuesRow_Private(&A->stash, rows[i], 1, cols + j, roworiented ? vals + i * nc + j : vals + i + j * nr, (PetscBool)(imode == ADD_VALUES))); 153 } 154 } 155 } 156 PetscFunctionReturn(PETSC_SUCCESS); 157 } 158 159 static PetscErrorCode MatMultXXXYYY_ScaLAPACK(Mat A, PetscBool transpose, PetscBool hermitian, PetscScalar beta, const PetscScalar *x, PetscScalar *y) 160 { 161 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 162 PetscScalar *x2d, *y2d, alpha = 1.0; 163 const PetscInt *ranges; 164 PetscBLASInt xdesc[9], ydesc[9], x2desc[9], y2desc[9], mb, nb, lszx, lszy, zero = 0, one = 1, xlld, ylld, info; 165 166 PetscFunctionBegin; 167 if (transpose) { 168 /* create ScaLAPACK descriptors for vectors (1d block distribution) */ 169 PetscCall(PetscLayoutGetRanges(A->rmap, &ranges)); 170 PetscCall(PetscBLASIntCast(ranges[1], &mb)); /* x block size */ 171 PetscCall(PetscBLASIntCast(PetscMax(1, A->rmap->n), &xlld)); 172 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(xdesc, &a->M, &one, &mb, &one, &zero, &zero, &a->grid->ictxcol, &xlld, &info)); 173 PetscCheckScaLapackInfo("descinit", info); 174 PetscCall(PetscLayoutGetRanges(A->cmap, &ranges)); 175 PetscCall(PetscBLASIntCast(ranges[1], &nb)); /* y block size */ 176 ylld = 1; 177 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(ydesc, &one, &a->N, &one, &nb, &zero, &zero, &a->grid->ictxrow, &ylld, &info)); 178 PetscCheckScaLapackInfo("descinit", info); 179 180 /* allocate 2d vectors */ 181 lszx = SCALAPACKnumroc_(&a->M, &a->mb, &a->grid->myrow, &a->rsrc, &a->grid->nprow); 182 lszy = SCALAPACKnumroc_(&a->N, &a->nb, &a->grid->mycol, &a->csrc, &a->grid->npcol); 183 PetscCall(PetscMalloc2(lszx, &x2d, lszy, &y2d)); 184 PetscCall(PetscBLASIntCast(PetscMax(1, lszx), &xlld)); 185 186 /* create ScaLAPACK descriptors for vectors (2d block distribution) */ 187 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(x2desc, &a->M, &one, &a->mb, &one, &zero, &zero, &a->grid->ictxt, &xlld, &info)); 188 PetscCheckScaLapackInfo("descinit", info); 189 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(y2desc, &one, &a->N, &one, &a->nb, &zero, &zero, &a->grid->ictxt, &ylld, &info)); 190 PetscCheckScaLapackInfo("descinit", info); 191 192 /* redistribute x as a column of a 2d matrix */ 193 PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &one, x, &one, &one, xdesc, x2d, &one, &one, x2desc, &a->grid->ictxcol)); 194 195 /* redistribute y as a row of a 2d matrix */ 196 if (beta != 0.0) PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&one, &a->N, y, &one, &one, ydesc, y2d, &one, &one, y2desc, &a->grid->ictxrow)); 197 198 /* call PBLAS subroutine */ 199 if (hermitian) PetscCallBLAS("PBLASgemv", PBLASgemv_("C", &a->M, &a->N, &alpha, a->loc, &one, &one, a->desc, x2d, &one, &one, x2desc, &one, &beta, y2d, &one, &one, y2desc, &one)); 200 else PetscCallBLAS("PBLASgemv", PBLASgemv_("T", &a->M, &a->N, &alpha, a->loc, &one, &one, a->desc, x2d, &one, &one, x2desc, &one, &beta, y2d, &one, &one, y2desc, &one)); 201 202 /* redistribute y from a row of a 2d matrix */ 203 PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&one, &a->N, y2d, &one, &one, y2desc, y, &one, &one, ydesc, &a->grid->ictxrow)); 204 205 } else { /* non-transpose */ 206 207 /* create ScaLAPACK descriptors for vectors (1d block distribution) */ 208 PetscCall(PetscLayoutGetRanges(A->cmap, &ranges)); 209 PetscCall(PetscBLASIntCast(ranges[1], &nb)); /* x block size */ 210 xlld = 1; 211 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(xdesc, &one, &a->N, &one, &nb, &zero, &zero, &a->grid->ictxrow, &xlld, &info)); 212 PetscCheckScaLapackInfo("descinit", info); 213 PetscCall(PetscLayoutGetRanges(A->rmap, &ranges)); 214 PetscCall(PetscBLASIntCast(ranges[1], &mb)); /* y block size */ 215 PetscCall(PetscBLASIntCast(PetscMax(1, A->rmap->n), &ylld)); 216 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(ydesc, &a->M, &one, &mb, &one, &zero, &zero, &a->grid->ictxcol, &ylld, &info)); 217 PetscCheckScaLapackInfo("descinit", info); 218 219 /* allocate 2d vectors */ 220 lszy = SCALAPACKnumroc_(&a->M, &a->mb, &a->grid->myrow, &a->rsrc, &a->grid->nprow); 221 lszx = SCALAPACKnumroc_(&a->N, &a->nb, &a->grid->mycol, &a->csrc, &a->grid->npcol); 222 PetscCall(PetscMalloc2(lszx, &x2d, lszy, &y2d)); 223 PetscCall(PetscBLASIntCast(PetscMax(1, lszy), &ylld)); 224 225 /* create ScaLAPACK descriptors for vectors (2d block distribution) */ 226 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(x2desc, &one, &a->N, &one, &a->nb, &zero, &zero, &a->grid->ictxt, &xlld, &info)); 227 PetscCheckScaLapackInfo("descinit", info); 228 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(y2desc, &a->M, &one, &a->mb, &one, &zero, &zero, &a->grid->ictxt, &ylld, &info)); 229 PetscCheckScaLapackInfo("descinit", info); 230 231 /* redistribute x as a row of a 2d matrix */ 232 PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&one, &a->N, x, &one, &one, xdesc, x2d, &one, &one, x2desc, &a->grid->ictxrow)); 233 234 /* redistribute y as a column of a 2d matrix */ 235 if (beta != 0.0) PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &one, y, &one, &one, ydesc, y2d, &one, &one, y2desc, &a->grid->ictxcol)); 236 237 /* call PBLAS subroutine */ 238 PetscCallBLAS("PBLASgemv", PBLASgemv_("N", &a->M, &a->N, &alpha, a->loc, &one, &one, a->desc, x2d, &one, &one, x2desc, &one, &beta, y2d, &one, &one, y2desc, &one)); 239 240 /* redistribute y from a column of a 2d matrix */ 241 PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &one, y2d, &one, &one, y2desc, y, &one, &one, ydesc, &a->grid->ictxcol)); 242 } 243 PetscCall(PetscFree2(x2d, y2d)); 244 PetscFunctionReturn(PETSC_SUCCESS); 245 } 246 247 static PetscErrorCode MatMult_ScaLAPACK(Mat A, Vec x, Vec y) 248 { 249 const PetscScalar *xarray; 250 PetscScalar *yarray; 251 252 PetscFunctionBegin; 253 PetscCall(VecGetArrayRead(x, &xarray)); 254 PetscCall(VecGetArray(y, &yarray)); 255 PetscCall(MatMultXXXYYY_ScaLAPACK(A, PETSC_FALSE, PETSC_FALSE, 0.0, xarray, yarray)); 256 PetscCall(VecRestoreArrayRead(x, &xarray)); 257 PetscCall(VecRestoreArray(y, &yarray)); 258 PetscFunctionReturn(PETSC_SUCCESS); 259 } 260 261 static PetscErrorCode MatMultTranspose_ScaLAPACK(Mat A, Vec x, Vec y) 262 { 263 const PetscScalar *xarray; 264 PetscScalar *yarray; 265 266 PetscFunctionBegin; 267 PetscCall(VecGetArrayRead(x, &xarray)); 268 PetscCall(VecGetArray(y, &yarray)); 269 PetscCall(MatMultXXXYYY_ScaLAPACK(A, PETSC_TRUE, PETSC_FALSE, 0.0, xarray, yarray)); 270 PetscCall(VecRestoreArrayRead(x, &xarray)); 271 PetscCall(VecRestoreArray(y, &yarray)); 272 PetscFunctionReturn(PETSC_SUCCESS); 273 } 274 275 static PetscErrorCode MatMultHermitianTranspose_ScaLAPACK(Mat A, Vec x, Vec y) 276 { 277 const PetscScalar *xarray; 278 PetscScalar *yarray; 279 280 PetscFunctionBegin; 281 PetscCall(VecGetArrayRead(x, &xarray)); 282 PetscCall(VecGetArrayWrite(y, &yarray)); 283 PetscCall(MatMultXXXYYY_ScaLAPACK(A, PETSC_TRUE, PETSC_TRUE, 0.0, xarray, yarray)); 284 PetscCall(VecRestoreArrayRead(x, &xarray)); 285 PetscCall(VecRestoreArrayWrite(y, &yarray)); 286 PetscFunctionReturn(PETSC_SUCCESS); 287 } 288 289 static PetscErrorCode MatMultAdd_ScaLAPACK(Mat A, Vec x, Vec y, Vec z) 290 { 291 const PetscScalar *xarray; 292 PetscScalar *zarray; 293 294 PetscFunctionBegin; 295 if (y != z) PetscCall(VecCopy(y, z)); 296 PetscCall(VecGetArrayRead(x, &xarray)); 297 PetscCall(VecGetArray(z, &zarray)); 298 PetscCall(MatMultXXXYYY_ScaLAPACK(A, PETSC_FALSE, PETSC_FALSE, 1.0, xarray, zarray)); 299 PetscCall(VecRestoreArrayRead(x, &xarray)); 300 PetscCall(VecRestoreArray(z, &zarray)); 301 PetscFunctionReturn(PETSC_SUCCESS); 302 } 303 304 static PetscErrorCode MatMultTransposeAdd_ScaLAPACK(Mat A, Vec x, Vec y, Vec z) 305 { 306 const PetscScalar *xarray; 307 PetscScalar *zarray; 308 309 PetscFunctionBegin; 310 if (y != z) PetscCall(VecCopy(y, z)); 311 PetscCall(VecGetArrayRead(x, &xarray)); 312 PetscCall(VecGetArray(z, &zarray)); 313 PetscCall(MatMultXXXYYY_ScaLAPACK(A, PETSC_TRUE, PETSC_FALSE, 1.0, xarray, zarray)); 314 PetscCall(VecRestoreArrayRead(x, &xarray)); 315 PetscCall(VecRestoreArray(z, &zarray)); 316 PetscFunctionReturn(PETSC_SUCCESS); 317 } 318 319 static PetscErrorCode MatMultHermitianTransposeAdd_ScaLAPACK(Mat A, Vec x, Vec y, Vec z) 320 { 321 const PetscScalar *xarray; 322 PetscScalar *zarray; 323 324 PetscFunctionBegin; 325 if (y != z) PetscCall(VecCopy(y, z)); 326 PetscCall(VecGetArrayRead(x, &xarray)); 327 PetscCall(VecGetArray(z, &zarray)); 328 PetscCall(MatMultXXXYYY_ScaLAPACK(A, PETSC_TRUE, PETSC_TRUE, 1.0, xarray, zarray)); 329 PetscCall(VecRestoreArrayRead(x, &xarray)); 330 PetscCall(VecRestoreArray(z, &zarray)); 331 PetscFunctionReturn(PETSC_SUCCESS); 332 } 333 334 PetscErrorCode MatMatMultNumeric_ScaLAPACK(Mat A, Mat B, Mat C) 335 { 336 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 337 Mat_ScaLAPACK *b = (Mat_ScaLAPACK *)B->data; 338 Mat_ScaLAPACK *c = (Mat_ScaLAPACK *)C->data; 339 PetscScalar sone = 1.0, zero = 0.0; 340 PetscBLASInt one = 1; 341 342 PetscFunctionBegin; 343 PetscCallBLAS("PBLASgemm", PBLASgemm_("N", "N", &a->M, &b->N, &a->N, &sone, a->loc, &one, &one, a->desc, b->loc, &one, &one, b->desc, &zero, c->loc, &one, &one, c->desc)); 344 C->assembled = PETSC_TRUE; 345 PetscFunctionReturn(PETSC_SUCCESS); 346 } 347 348 PetscErrorCode MatMatMultSymbolic_ScaLAPACK(Mat A, Mat B, PetscReal fill, Mat C) 349 { 350 PetscFunctionBegin; 351 PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, PETSC_DECIDE, PETSC_DECIDE)); 352 PetscCall(MatSetType(C, MATSCALAPACK)); 353 PetscCall(MatSetUp(C)); 354 C->ops->matmultnumeric = MatMatMultNumeric_ScaLAPACK; 355 PetscFunctionReturn(PETSC_SUCCESS); 356 } 357 358 static PetscErrorCode MatTransposeMatMultNumeric_ScaLAPACK(Mat A, Mat B, Mat C) 359 { 360 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 361 Mat_ScaLAPACK *b = (Mat_ScaLAPACK *)B->data; 362 Mat_ScaLAPACK *c = (Mat_ScaLAPACK *)C->data; 363 PetscScalar sone = 1.0, zero = 0.0; 364 PetscBLASInt one = 1; 365 366 PetscFunctionBegin; 367 PetscCallBLAS("PBLASgemm", PBLASgemm_("T", "N", &a->N, &b->N, &a->M, &sone, a->loc, &one, &one, a->desc, b->loc, &one, &one, b->desc, &zero, c->loc, &one, &one, c->desc)); 368 C->assembled = PETSC_TRUE; 369 PetscFunctionReturn(PETSC_SUCCESS); 370 } 371 372 static PetscErrorCode MatTransposeMatMultSymbolic_ScaLAPACK(Mat A, Mat B, PetscReal fill, Mat C) 373 { 374 PetscFunctionBegin; 375 PetscCall(MatSetSizes(C, A->cmap->n, B->cmap->n, PETSC_DECIDE, PETSC_DECIDE)); 376 PetscCall(MatSetType(C, MATSCALAPACK)); 377 PetscCall(MatSetUp(C)); 378 C->ops->transposematmultnumeric = MatTransposeMatMultNumeric_ScaLAPACK; 379 PetscFunctionReturn(PETSC_SUCCESS); 380 } 381 382 static PetscErrorCode MatMatTransposeMultNumeric_ScaLAPACK(Mat A, Mat B, Mat C) 383 { 384 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 385 Mat_ScaLAPACK *b = (Mat_ScaLAPACK *)B->data; 386 Mat_ScaLAPACK *c = (Mat_ScaLAPACK *)C->data; 387 PetscScalar sone = 1.0, zero = 0.0; 388 PetscBLASInt one = 1; 389 390 PetscFunctionBegin; 391 PetscCallBLAS("PBLASgemm", PBLASgemm_("N", "T", &a->M, &b->M, &a->N, &sone, a->loc, &one, &one, a->desc, b->loc, &one, &one, b->desc, &zero, c->loc, &one, &one, c->desc)); 392 C->assembled = PETSC_TRUE; 393 PetscFunctionReturn(PETSC_SUCCESS); 394 } 395 396 static PetscErrorCode MatMatTransposeMultSymbolic_ScaLAPACK(Mat A, Mat B, PetscReal fill, Mat C) 397 { 398 PetscFunctionBegin; 399 PetscCall(MatSetSizes(C, A->rmap->n, B->rmap->n, PETSC_DECIDE, PETSC_DECIDE)); 400 PetscCall(MatSetType(C, MATSCALAPACK)); 401 PetscCall(MatSetUp(C)); 402 C->ops->mattransposemultnumeric = MatMatTransposeMultNumeric_ScaLAPACK; 403 PetscFunctionReturn(PETSC_SUCCESS); 404 } 405 406 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_ScaLAPACK(Mat C) 407 { 408 Mat_Product *product = C->product; 409 410 PetscFunctionBegin; 411 switch (product->type) { 412 case MATPRODUCT_AB: 413 C->ops->matmultsymbolic = MatMatMultSymbolic_ScaLAPACK; 414 C->ops->productsymbolic = MatProductSymbolic_AB; 415 break; 416 case MATPRODUCT_AtB: 417 C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_ScaLAPACK; 418 C->ops->productsymbolic = MatProductSymbolic_AtB; 419 break; 420 case MATPRODUCT_ABt: 421 C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_ScaLAPACK; 422 C->ops->productsymbolic = MatProductSymbolic_ABt; 423 break; 424 default: 425 SETERRQ(PetscObjectComm((PetscObject)C), PETSC_ERR_SUP, "MatProduct type %s is not supported for ScaLAPACK and ScaLAPACK matrices", MatProductTypes[product->type]); 426 } 427 PetscFunctionReturn(PETSC_SUCCESS); 428 } 429 430 static PetscErrorCode MatGetDiagonal_ScaLAPACK(Mat A, Vec D) 431 { 432 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 433 PetscScalar *darray, *d2d, v; 434 const PetscInt *ranges; 435 PetscBLASInt j, ddesc[9], d2desc[9], mb, nb, lszd, zero = 0, one = 1, dlld, info; 436 437 PetscFunctionBegin; 438 PetscCall(VecGetArray(D, &darray)); 439 440 if (A->rmap->N <= A->cmap->N) { /* row version */ 441 442 /* create ScaLAPACK descriptor for vector (1d block distribution) */ 443 PetscCall(PetscLayoutGetRanges(A->rmap, &ranges)); 444 PetscCall(PetscBLASIntCast(ranges[1], &mb)); /* D block size */ 445 PetscCall(PetscBLASIntCast(PetscMax(1, A->rmap->n), &dlld)); 446 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(ddesc, &a->M, &one, &mb, &one, &zero, &zero, &a->grid->ictxcol, &dlld, &info)); 447 PetscCheckScaLapackInfo("descinit", info); 448 449 /* allocate 2d vector */ 450 lszd = SCALAPACKnumroc_(&a->M, &a->mb, &a->grid->myrow, &a->rsrc, &a->grid->nprow); 451 PetscCall(PetscCalloc1(lszd, &d2d)); 452 PetscCall(PetscBLASIntCast(PetscMax(1, lszd), &dlld)); 453 454 /* create ScaLAPACK descriptor for vector (2d block distribution) */ 455 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(d2desc, &a->M, &one, &a->mb, &one, &zero, &zero, &a->grid->ictxt, &dlld, &info)); 456 PetscCheckScaLapackInfo("descinit", info); 457 458 /* collect diagonal */ 459 for (j = 1; j <= a->M; j++) { 460 PetscCallBLAS("SCALAPACKelget", SCALAPACKelget_("R", " ", &v, a->loc, &j, &j, a->desc)); 461 PetscCallBLAS("SCALAPACKelset", SCALAPACKelset_(d2d, &j, &one, d2desc, &v)); 462 } 463 464 /* redistribute d from a column of a 2d matrix */ 465 PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &one, d2d, &one, &one, d2desc, darray, &one, &one, ddesc, &a->grid->ictxcol)); 466 PetscCall(PetscFree(d2d)); 467 468 } else { /* column version */ 469 470 /* create ScaLAPACK descriptor for vector (1d block distribution) */ 471 PetscCall(PetscLayoutGetRanges(A->cmap, &ranges)); 472 PetscCall(PetscBLASIntCast(ranges[1], &nb)); /* D block size */ 473 dlld = 1; 474 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(ddesc, &one, &a->N, &one, &nb, &zero, &zero, &a->grid->ictxrow, &dlld, &info)); 475 PetscCheckScaLapackInfo("descinit", info); 476 477 /* allocate 2d vector */ 478 lszd = SCALAPACKnumroc_(&a->N, &a->nb, &a->grid->mycol, &a->csrc, &a->grid->npcol); 479 PetscCall(PetscCalloc1(lszd, &d2d)); 480 481 /* create ScaLAPACK descriptor for vector (2d block distribution) */ 482 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(d2desc, &one, &a->N, &one, &a->nb, &zero, &zero, &a->grid->ictxt, &dlld, &info)); 483 PetscCheckScaLapackInfo("descinit", info); 484 485 /* collect diagonal */ 486 for (j = 1; j <= a->N; j++) { 487 PetscCallBLAS("SCALAPACKelget", SCALAPACKelget_("C", " ", &v, a->loc, &j, &j, a->desc)); 488 PetscCallBLAS("SCALAPACKelset", SCALAPACKelset_(d2d, &one, &j, d2desc, &v)); 489 } 490 491 /* redistribute d from a row of a 2d matrix */ 492 PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&one, &a->N, d2d, &one, &one, d2desc, darray, &one, &one, ddesc, &a->grid->ictxrow)); 493 PetscCall(PetscFree(d2d)); 494 } 495 496 PetscCall(VecRestoreArray(D, &darray)); 497 PetscCall(VecAssemblyBegin(D)); 498 PetscCall(VecAssemblyEnd(D)); 499 PetscFunctionReturn(PETSC_SUCCESS); 500 } 501 502 static PetscErrorCode MatDiagonalScale_ScaLAPACK(Mat A, Vec L, Vec R) 503 { 504 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 505 const PetscScalar *d; 506 const PetscInt *ranges; 507 PetscScalar *d2d; 508 PetscBLASInt i, j, ddesc[9], d2desc[9], mb, nb, lszd, zero = 0, one = 1, dlld, info; 509 510 PetscFunctionBegin; 511 if (R) { 512 PetscCall(VecGetArrayRead(R, &d)); 513 /* create ScaLAPACK descriptor for vector (1d block distribution) */ 514 PetscCall(PetscLayoutGetRanges(A->cmap, &ranges)); 515 PetscCall(PetscBLASIntCast(ranges[1], &nb)); /* D block size */ 516 dlld = 1; 517 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(ddesc, &one, &a->N, &one, &nb, &zero, &zero, &a->grid->ictxrow, &dlld, &info)); 518 PetscCheckScaLapackInfo("descinit", info); 519 520 /* allocate 2d vector */ 521 lszd = SCALAPACKnumroc_(&a->N, &a->nb, &a->grid->mycol, &a->csrc, &a->grid->npcol); 522 PetscCall(PetscCalloc1(lszd, &d2d)); 523 524 /* create ScaLAPACK descriptor for vector (2d block distribution) */ 525 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(d2desc, &one, &a->N, &one, &a->nb, &zero, &zero, &a->grid->ictxt, &dlld, &info)); 526 PetscCheckScaLapackInfo("descinit", info); 527 528 /* redistribute d to a row of a 2d matrix */ 529 PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&one, &a->N, d, &one, &one, ddesc, d2d, &one, &one, d2desc, &a->grid->ictxrow)); 530 531 /* broadcast along process columns */ 532 if (!a->grid->myrow) Cdgebs2d(a->grid->ictxt, "C", " ", 1, lszd, d2d, dlld); 533 else Cdgebr2d(a->grid->ictxt, "C", " ", 1, lszd, d2d, dlld, 0, a->grid->mycol); 534 535 /* local scaling */ 536 for (j = 0; j < a->locc; j++) 537 for (i = 0; i < a->locr; i++) a->loc[i + j * a->lld] *= d2d[j]; 538 539 PetscCall(PetscFree(d2d)); 540 PetscCall(VecRestoreArrayRead(R, &d)); 541 } 542 if (L) { 543 PetscCall(VecGetArrayRead(L, &d)); 544 /* create ScaLAPACK descriptor for vector (1d block distribution) */ 545 PetscCall(PetscLayoutGetRanges(A->rmap, &ranges)); 546 PetscCall(PetscBLASIntCast(ranges[1], &mb)); /* D block size */ 547 PetscCall(PetscBLASIntCast(PetscMax(1, A->rmap->n), &dlld)); 548 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(ddesc, &a->M, &one, &mb, &one, &zero, &zero, &a->grid->ictxcol, &dlld, &info)); 549 PetscCheckScaLapackInfo("descinit", info); 550 551 /* allocate 2d vector */ 552 lszd = SCALAPACKnumroc_(&a->M, &a->mb, &a->grid->myrow, &a->rsrc, &a->grid->nprow); 553 PetscCall(PetscCalloc1(lszd, &d2d)); 554 PetscCall(PetscBLASIntCast(PetscMax(1, lszd), &dlld)); 555 556 /* create ScaLAPACK descriptor for vector (2d block distribution) */ 557 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(d2desc, &a->M, &one, &a->mb, &one, &zero, &zero, &a->grid->ictxt, &dlld, &info)); 558 PetscCheckScaLapackInfo("descinit", info); 559 560 /* redistribute d to a column of a 2d matrix */ 561 PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &one, d, &one, &one, ddesc, d2d, &one, &one, d2desc, &a->grid->ictxcol)); 562 563 /* broadcast along process rows */ 564 if (!a->grid->mycol) Cdgebs2d(a->grid->ictxt, "R", " ", lszd, 1, d2d, dlld); 565 else Cdgebr2d(a->grid->ictxt, "R", " ", lszd, 1, d2d, dlld, a->grid->myrow, 0); 566 567 /* local scaling */ 568 for (i = 0; i < a->locr; i++) 569 for (j = 0; j < a->locc; j++) a->loc[i + j * a->lld] *= d2d[i]; 570 571 PetscCall(PetscFree(d2d)); 572 PetscCall(VecRestoreArrayRead(L, &d)); 573 } 574 PetscFunctionReturn(PETSC_SUCCESS); 575 } 576 577 static PetscErrorCode MatMissingDiagonal_ScaLAPACK(Mat A, PetscBool *missing, PetscInt *d) 578 { 579 PetscFunctionBegin; 580 *missing = PETSC_FALSE; 581 PetscFunctionReturn(PETSC_SUCCESS); 582 } 583 584 static PetscErrorCode MatScale_ScaLAPACK(Mat X, PetscScalar a) 585 { 586 Mat_ScaLAPACK *x = (Mat_ScaLAPACK *)X->data; 587 PetscBLASInt n, one = 1; 588 589 PetscFunctionBegin; 590 n = x->lld * x->locc; 591 PetscCallBLAS("BLASscal", BLASscal_(&n, &a, x->loc, &one)); 592 PetscFunctionReturn(PETSC_SUCCESS); 593 } 594 595 static PetscErrorCode MatShift_ScaLAPACK(Mat X, PetscScalar alpha) 596 { 597 Mat_ScaLAPACK *x = (Mat_ScaLAPACK *)X->data; 598 PetscBLASInt i, n; 599 PetscScalar v; 600 601 PetscFunctionBegin; 602 n = PetscMin(x->M, x->N); 603 for (i = 1; i <= n; i++) { 604 PetscCallBLAS("SCALAPACKelget", SCALAPACKelget_("-", " ", &v, x->loc, &i, &i, x->desc)); 605 v += alpha; 606 PetscCallBLAS("SCALAPACKelset", SCALAPACKelset_(x->loc, &i, &i, x->desc, &v)); 607 } 608 PetscFunctionReturn(PETSC_SUCCESS); 609 } 610 611 static PetscErrorCode MatAXPY_ScaLAPACK(Mat Y, PetscScalar alpha, Mat X, MatStructure str) 612 { 613 Mat_ScaLAPACK *x = (Mat_ScaLAPACK *)X->data; 614 Mat_ScaLAPACK *y = (Mat_ScaLAPACK *)Y->data; 615 PetscBLASInt one = 1; 616 PetscScalar beta = 1.0; 617 618 PetscFunctionBegin; 619 MatScaLAPACKCheckDistribution(Y, 1, X, 3); 620 PetscCallBLAS("SCALAPACKmatadd", SCALAPACKmatadd_(&x->M, &x->N, &alpha, x->loc, &one, &one, x->desc, &beta, y->loc, &one, &one, y->desc)); 621 PetscCall(PetscObjectStateIncrease((PetscObject)Y)); 622 PetscFunctionReturn(PETSC_SUCCESS); 623 } 624 625 static PetscErrorCode MatCopy_ScaLAPACK(Mat A, Mat B, MatStructure str) 626 { 627 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 628 Mat_ScaLAPACK *b = (Mat_ScaLAPACK *)B->data; 629 630 PetscFunctionBegin; 631 PetscCall(PetscArraycpy(b->loc, a->loc, a->lld * a->locc)); 632 PetscCall(PetscObjectStateIncrease((PetscObject)B)); 633 PetscFunctionReturn(PETSC_SUCCESS); 634 } 635 636 static PetscErrorCode MatDuplicate_ScaLAPACK(Mat A, MatDuplicateOption op, Mat *B) 637 { 638 Mat Bs; 639 MPI_Comm comm; 640 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data, *b; 641 642 PetscFunctionBegin; 643 PetscCall(PetscObjectGetComm((PetscObject)A, &comm)); 644 PetscCall(MatCreate(comm, &Bs)); 645 PetscCall(MatSetSizes(Bs, A->rmap->n, A->cmap->n, PETSC_DECIDE, PETSC_DECIDE)); 646 PetscCall(MatSetType(Bs, MATSCALAPACK)); 647 b = (Mat_ScaLAPACK *)Bs->data; 648 b->M = a->M; 649 b->N = a->N; 650 b->mb = a->mb; 651 b->nb = a->nb; 652 b->rsrc = a->rsrc; 653 b->csrc = a->csrc; 654 PetscCall(MatSetUp(Bs)); 655 *B = Bs; 656 if (op == MAT_COPY_VALUES) PetscCall(PetscArraycpy(b->loc, a->loc, a->lld * a->locc)); 657 Bs->assembled = PETSC_TRUE; 658 PetscFunctionReturn(PETSC_SUCCESS); 659 } 660 661 static PetscErrorCode MatTranspose_ScaLAPACK(Mat A, MatReuse reuse, Mat *B) 662 { 663 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data, *b; 664 Mat Bs = *B; 665 PetscBLASInt one = 1; 666 PetscScalar sone = 1.0, zero = 0.0; 667 #if defined(PETSC_USE_COMPLEX) 668 PetscInt i, j; 669 #endif 670 671 PetscFunctionBegin; 672 if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B)); 673 PetscCheck(reuse == MAT_INITIAL_MATRIX, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only MAT_INITIAL_MATRIX supported"); 674 PetscCall(MatCreateScaLAPACK(PetscObjectComm((PetscObject)A), a->nb, a->mb, a->N, a->M, a->csrc, a->rsrc, &Bs)); 675 *B = Bs; 676 b = (Mat_ScaLAPACK *)Bs->data; 677 PetscCallBLAS("PBLAStran", PBLAStran_(&a->N, &a->M, &sone, a->loc, &one, &one, a->desc, &zero, b->loc, &one, &one, b->desc)); 678 #if defined(PETSC_USE_COMPLEX) 679 /* undo conjugation */ 680 for (i = 0; i < b->locr; i++) 681 for (j = 0; j < b->locc; j++) b->loc[i + j * b->lld] = PetscConj(b->loc[i + j * b->lld]); 682 #endif 683 Bs->assembled = PETSC_TRUE; 684 PetscFunctionReturn(PETSC_SUCCESS); 685 } 686 687 static PetscErrorCode MatConjugate_ScaLAPACK(Mat A) 688 { 689 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 690 PetscInt i, j; 691 692 PetscFunctionBegin; 693 for (i = 0; i < a->locr; i++) 694 for (j = 0; j < a->locc; j++) a->loc[i + j * a->lld] = PetscConj(a->loc[i + j * a->lld]); 695 PetscFunctionReturn(PETSC_SUCCESS); 696 } 697 698 static PetscErrorCode MatHermitianTranspose_ScaLAPACK(Mat A, MatReuse reuse, Mat *B) 699 { 700 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data, *b; 701 Mat Bs = *B; 702 PetscBLASInt one = 1; 703 PetscScalar sone = 1.0, zero = 0.0; 704 705 PetscFunctionBegin; 706 PetscCheck(reuse == MAT_INITIAL_MATRIX, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only MAT_INITIAL_MATRIX supported"); 707 PetscCall(MatCreateScaLAPACK(PetscObjectComm((PetscObject)A), a->nb, a->mb, a->N, a->M, a->csrc, a->rsrc, &Bs)); 708 *B = Bs; 709 b = (Mat_ScaLAPACK *)Bs->data; 710 PetscCallBLAS("PBLAStran", PBLAStran_(&a->N, &a->M, &sone, a->loc, &one, &one, a->desc, &zero, b->loc, &one, &one, b->desc)); 711 Bs->assembled = PETSC_TRUE; 712 PetscFunctionReturn(PETSC_SUCCESS); 713 } 714 715 static PetscErrorCode MatSolve_ScaLAPACK(Mat A, Vec B, Vec X) 716 { 717 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 718 PetscScalar *x, *x2d; 719 const PetscInt *ranges; 720 PetscBLASInt xdesc[9], x2desc[9], mb, lszx, zero = 0, one = 1, xlld, nrhs = 1, info; 721 722 PetscFunctionBegin; 723 PetscCall(VecCopy(B, X)); 724 PetscCall(VecGetArray(X, &x)); 725 726 /* create ScaLAPACK descriptor for a vector (1d block distribution) */ 727 PetscCall(PetscLayoutGetRanges(A->rmap, &ranges)); 728 PetscCall(PetscBLASIntCast(ranges[1], &mb)); /* x block size */ 729 PetscCall(PetscBLASIntCast(PetscMax(1, A->rmap->n), &xlld)); 730 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(xdesc, &a->M, &one, &mb, &one, &zero, &zero, &a->grid->ictxcol, &xlld, &info)); 731 PetscCheckScaLapackInfo("descinit", info); 732 733 /* allocate 2d vector */ 734 lszx = SCALAPACKnumroc_(&a->M, &a->mb, &a->grid->myrow, &a->rsrc, &a->grid->nprow); 735 PetscCall(PetscMalloc1(lszx, &x2d)); 736 PetscCall(PetscBLASIntCast(PetscMax(1, lszx), &xlld)); 737 738 /* create ScaLAPACK descriptor for a vector (2d block distribution) */ 739 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(x2desc, &a->M, &one, &a->mb, &one, &zero, &zero, &a->grid->ictxt, &xlld, &info)); 740 PetscCheckScaLapackInfo("descinit", info); 741 742 /* redistribute x as a column of a 2d matrix */ 743 PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &one, x, &one, &one, xdesc, x2d, &one, &one, x2desc, &a->grid->ictxcol)); 744 745 /* call ScaLAPACK subroutine */ 746 switch (A->factortype) { 747 case MAT_FACTOR_LU: 748 PetscCallBLAS("SCALAPACKgetrs", SCALAPACKgetrs_("N", &a->M, &nrhs, a->loc, &one, &one, a->desc, a->pivots, x2d, &one, &one, x2desc, &info)); 749 PetscCheckScaLapackInfo("getrs", info); 750 break; 751 case MAT_FACTOR_CHOLESKY: 752 PetscCallBLAS("SCALAPACKpotrs", SCALAPACKpotrs_("L", &a->M, &nrhs, a->loc, &one, &one, a->desc, x2d, &one, &one, x2desc, &info)); 753 PetscCheckScaLapackInfo("potrs", info); 754 break; 755 default: 756 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unfactored Matrix or Unsupported MatFactorType"); 757 } 758 759 /* redistribute x from a column of a 2d matrix */ 760 PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &one, x2d, &one, &one, x2desc, x, &one, &one, xdesc, &a->grid->ictxcol)); 761 762 PetscCall(PetscFree(x2d)); 763 PetscCall(VecRestoreArray(X, &x)); 764 PetscFunctionReturn(PETSC_SUCCESS); 765 } 766 767 static PetscErrorCode MatSolveAdd_ScaLAPACK(Mat A, Vec B, Vec Y, Vec X) 768 { 769 PetscFunctionBegin; 770 PetscCall(MatSolve_ScaLAPACK(A, B, X)); 771 PetscCall(VecAXPY(X, 1, Y)); 772 PetscFunctionReturn(PETSC_SUCCESS); 773 } 774 775 static PetscErrorCode MatMatSolve_ScaLAPACK(Mat A, Mat B, Mat X) 776 { 777 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data, *x; 778 PetscBool flg1, flg2; 779 PetscBLASInt one = 1, info; 780 Mat C; 781 MatType type; 782 783 PetscFunctionBegin; 784 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSCALAPACK, &flg1)); 785 PetscCall(PetscObjectTypeCompare((PetscObject)X, MATSCALAPACK, &flg2)); 786 if (flg1 && flg2) MatScaLAPACKCheckDistribution(B, 2, X, 3); 787 if (flg2) { 788 if (flg1) PetscCall(MatCopy(B, X, SAME_NONZERO_PATTERN)); 789 else PetscCall(MatConvert(B, MATSCALAPACK, MAT_REUSE_MATRIX, &X)); 790 C = X; 791 } else { 792 PetscCall(MatConvert(B, MATSCALAPACK, MAT_INITIAL_MATRIX, &C)); 793 } 794 x = (Mat_ScaLAPACK *)C->data; 795 796 switch (A->factortype) { 797 case MAT_FACTOR_LU: 798 PetscCallBLAS("SCALAPACKgetrs", SCALAPACKgetrs_("N", &a->M, &x->N, a->loc, &one, &one, a->desc, a->pivots, x->loc, &one, &one, x->desc, &info)); 799 PetscCheckScaLapackInfo("getrs", info); 800 break; 801 case MAT_FACTOR_CHOLESKY: 802 PetscCallBLAS("SCALAPACKpotrs", SCALAPACKpotrs_("L", &a->M, &x->N, a->loc, &one, &one, a->desc, x->loc, &one, &one, x->desc, &info)); 803 PetscCheckScaLapackInfo("potrs", info); 804 break; 805 default: 806 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unfactored Matrix or Unsupported MatFactorType"); 807 } 808 if (!flg2) { 809 PetscCall(MatGetType(X, &type)); 810 PetscCall(MatConvert(C, type, MAT_REUSE_MATRIX, &X)); 811 PetscCall(MatDestroy(&C)); 812 } 813 PetscFunctionReturn(PETSC_SUCCESS); 814 } 815 816 static PetscErrorCode MatLUFactor_ScaLAPACK(Mat A, IS row, IS col, const MatFactorInfo *factorinfo) 817 { 818 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 819 PetscBLASInt one = 1, info; 820 821 PetscFunctionBegin; 822 if (!a->pivots) PetscCall(PetscMalloc1(a->locr + a->mb, &a->pivots)); 823 PetscCallBLAS("SCALAPACKgetrf", SCALAPACKgetrf_(&a->M, &a->N, a->loc, &one, &one, a->desc, a->pivots, &info)); 824 PetscCheckScaLapackInfo("getrf", info); 825 A->factortype = MAT_FACTOR_LU; 826 A->assembled = PETSC_TRUE; 827 828 PetscCall(PetscFree(A->solvertype)); 829 PetscCall(PetscStrallocpy(MATSOLVERSCALAPACK, &A->solvertype)); 830 PetscFunctionReturn(PETSC_SUCCESS); 831 } 832 833 static PetscErrorCode MatLUFactorNumeric_ScaLAPACK(Mat F, Mat A, const MatFactorInfo *info) 834 { 835 PetscFunctionBegin; 836 PetscCall(MatCopy(A, F, SAME_NONZERO_PATTERN)); 837 PetscCall(MatLUFactor_ScaLAPACK(F, 0, 0, info)); 838 PetscFunctionReturn(PETSC_SUCCESS); 839 } 840 841 static PetscErrorCode MatLUFactorSymbolic_ScaLAPACK(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info) 842 { 843 PetscFunctionBegin; 844 /* F is created and allocated by MatGetFactor_scalapack_petsc(), skip this routine. */ 845 PetscFunctionReturn(PETSC_SUCCESS); 846 } 847 848 static PetscErrorCode MatCholeskyFactor_ScaLAPACK(Mat A, IS perm, const MatFactorInfo *factorinfo) 849 { 850 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 851 PetscBLASInt one = 1, info; 852 853 PetscFunctionBegin; 854 PetscCallBLAS("SCALAPACKpotrf", SCALAPACKpotrf_("L", &a->M, a->loc, &one, &one, a->desc, &info)); 855 PetscCheckScaLapackInfo("potrf", info); 856 A->factortype = MAT_FACTOR_CHOLESKY; 857 A->assembled = PETSC_TRUE; 858 859 PetscCall(PetscFree(A->solvertype)); 860 PetscCall(PetscStrallocpy(MATSOLVERSCALAPACK, &A->solvertype)); 861 PetscFunctionReturn(PETSC_SUCCESS); 862 } 863 864 static PetscErrorCode MatCholeskyFactorNumeric_ScaLAPACK(Mat F, Mat A, const MatFactorInfo *info) 865 { 866 PetscFunctionBegin; 867 PetscCall(MatCopy(A, F, SAME_NONZERO_PATTERN)); 868 PetscCall(MatCholeskyFactor_ScaLAPACK(F, 0, info)); 869 PetscFunctionReturn(PETSC_SUCCESS); 870 } 871 872 static PetscErrorCode MatCholeskyFactorSymbolic_ScaLAPACK(Mat F, Mat A, IS perm, const MatFactorInfo *info) 873 { 874 PetscFunctionBegin; 875 /* F is created and allocated by MatGetFactor_scalapack_petsc(), skip this routine. */ 876 PetscFunctionReturn(PETSC_SUCCESS); 877 } 878 879 static PetscErrorCode MatFactorGetSolverType_scalapack_scalapack(Mat A, MatSolverType *type) 880 { 881 PetscFunctionBegin; 882 *type = MATSOLVERSCALAPACK; 883 PetscFunctionReturn(PETSC_SUCCESS); 884 } 885 886 static PetscErrorCode MatGetFactor_scalapack_scalapack(Mat A, MatFactorType ftype, Mat *F) 887 { 888 Mat B; 889 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 890 891 PetscFunctionBegin; 892 /* Create the factorization matrix */ 893 PetscCall(MatCreateScaLAPACK(PetscObjectComm((PetscObject)A), a->mb, a->nb, a->M, a->N, a->rsrc, a->csrc, &B)); 894 B->trivialsymbolic = PETSC_TRUE; 895 B->factortype = ftype; 896 PetscCall(PetscFree(B->solvertype)); 897 PetscCall(PetscStrallocpy(MATSOLVERSCALAPACK, &B->solvertype)); 898 899 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_scalapack_scalapack)); 900 *F = B; 901 PetscFunctionReturn(PETSC_SUCCESS); 902 } 903 904 PETSC_INTERN PetscErrorCode MatSolverTypeRegister_ScaLAPACK(void) 905 { 906 PetscFunctionBegin; 907 PetscCall(MatSolverTypeRegister(MATSOLVERSCALAPACK, MATSCALAPACK, MAT_FACTOR_LU, MatGetFactor_scalapack_scalapack)); 908 PetscCall(MatSolverTypeRegister(MATSOLVERSCALAPACK, MATSCALAPACK, MAT_FACTOR_CHOLESKY, MatGetFactor_scalapack_scalapack)); 909 PetscFunctionReturn(PETSC_SUCCESS); 910 } 911 912 static PetscErrorCode MatNorm_ScaLAPACK(Mat A, NormType type, PetscReal *nrm) 913 { 914 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 915 PetscBLASInt one = 1, lwork = 0; 916 const char *ntype; 917 PetscScalar *work = NULL, dummy; 918 919 PetscFunctionBegin; 920 switch (type) { 921 case NORM_1: 922 ntype = "1"; 923 lwork = PetscMax(a->locr, a->locc); 924 break; 925 case NORM_FROBENIUS: 926 ntype = "F"; 927 work = &dummy; 928 break; 929 case NORM_INFINITY: 930 ntype = "I"; 931 lwork = PetscMax(a->locr, a->locc); 932 break; 933 default: 934 SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Unsupported norm type"); 935 } 936 if (lwork) PetscCall(PetscMalloc1(lwork, &work)); 937 *nrm = SCALAPACKlange_(ntype, &a->M, &a->N, a->loc, &one, &one, a->desc, work); 938 if (lwork) PetscCall(PetscFree(work)); 939 PetscFunctionReturn(PETSC_SUCCESS); 940 } 941 942 static PetscErrorCode MatZeroEntries_ScaLAPACK(Mat A) 943 { 944 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 945 946 PetscFunctionBegin; 947 PetscCall(PetscArrayzero(a->loc, a->lld * a->locc)); 948 PetscFunctionReturn(PETSC_SUCCESS); 949 } 950 951 static PetscErrorCode MatGetOwnershipIS_ScaLAPACK(Mat A, IS *rows, IS *cols) 952 { 953 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 954 PetscInt i, n, nb, isrc, nproc, iproc, *idx; 955 956 PetscFunctionBegin; 957 if (rows) { 958 n = a->locr; 959 nb = a->mb; 960 isrc = a->rsrc; 961 nproc = a->grid->nprow; 962 iproc = a->grid->myrow; 963 PetscCall(PetscMalloc1(n, &idx)); 964 for (i = 0; i < n; i++) idx[i] = nproc * nb * (i / nb) + i % nb + ((nproc + iproc - isrc) % nproc) * nb; 965 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_OWN_POINTER, rows)); 966 } 967 if (cols) { 968 n = a->locc; 969 nb = a->nb; 970 isrc = a->csrc; 971 nproc = a->grid->npcol; 972 iproc = a->grid->mycol; 973 PetscCall(PetscMalloc1(n, &idx)); 974 for (i = 0; i < n; i++) idx[i] = nproc * nb * (i / nb) + i % nb + ((nproc + iproc - isrc) % nproc) * nb; 975 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_OWN_POINTER, cols)); 976 } 977 PetscFunctionReturn(PETSC_SUCCESS); 978 } 979 980 static PetscErrorCode MatConvert_ScaLAPACK_Dense(Mat A, MatType newtype, MatReuse reuse, Mat *B) 981 { 982 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 983 Mat Bmpi; 984 MPI_Comm comm; 985 PetscInt i, M = A->rmap->N, N = A->cmap->N, m, n, rstart, rend, nz, ldb; 986 const PetscInt *ranges, *branges, *cwork; 987 const PetscScalar *vwork; 988 PetscBLASInt bdesc[9], bmb, zero = 0, one = 1, lld, info; 989 PetscScalar *barray; 990 PetscBool differ = PETSC_FALSE; 991 PetscMPIInt size; 992 993 PetscFunctionBegin; 994 PetscCall(PetscObjectGetComm((PetscObject)A, &comm)); 995 PetscCall(PetscLayoutGetRanges(A->rmap, &ranges)); 996 997 if (reuse == MAT_REUSE_MATRIX) { /* check if local sizes differ in A and B */ 998 PetscCallMPI(MPI_Comm_size(comm, &size)); 999 PetscCall(PetscLayoutGetRanges((*B)->rmap, &branges)); 1000 for (i = 0; i < size; i++) 1001 if (ranges[i + 1] != branges[i + 1]) { 1002 differ = PETSC_TRUE; 1003 break; 1004 } 1005 } 1006 1007 if (reuse == MAT_REUSE_MATRIX && differ) { /* special case, use auxiliary dense matrix */ 1008 PetscCall(MatCreate(comm, &Bmpi)); 1009 m = PETSC_DECIDE; 1010 PetscCall(PetscSplitOwnershipEqual(comm, &m, &M)); 1011 n = PETSC_DECIDE; 1012 PetscCall(PetscSplitOwnershipEqual(comm, &n, &N)); 1013 PetscCall(MatSetSizes(Bmpi, m, n, M, N)); 1014 PetscCall(MatSetType(Bmpi, MATDENSE)); 1015 PetscCall(MatSetUp(Bmpi)); 1016 1017 /* create ScaLAPACK descriptor for B (1d block distribution) */ 1018 PetscCall(PetscBLASIntCast(ranges[1], &bmb)); /* row block size */ 1019 PetscCall(MatDenseGetLDA(Bmpi, &ldb)); 1020 PetscCall(PetscBLASIntCast(PetscMax(ldb, 1), &lld)); /* local leading dimension */ 1021 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(bdesc, &a->M, &a->N, &bmb, &a->N, &zero, &zero, &a->grid->ictxcol, &lld, &info)); 1022 PetscCheckScaLapackInfo("descinit", info); 1023 1024 /* redistribute matrix */ 1025 PetscCall(MatDenseGetArray(Bmpi, &barray)); 1026 PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &a->N, a->loc, &one, &one, a->desc, barray, &one, &one, bdesc, &a->grid->ictxcol)); 1027 PetscCall(MatDenseRestoreArray(Bmpi, &barray)); 1028 PetscCall(MatAssemblyBegin(Bmpi, MAT_FINAL_ASSEMBLY)); 1029 PetscCall(MatAssemblyEnd(Bmpi, MAT_FINAL_ASSEMBLY)); 1030 1031 /* transfer rows of auxiliary matrix to the final matrix B */ 1032 PetscCall(MatGetOwnershipRange(Bmpi, &rstart, &rend)); 1033 for (i = rstart; i < rend; i++) { 1034 PetscCall(MatGetRow(Bmpi, i, &nz, &cwork, &vwork)); 1035 PetscCall(MatSetValues(*B, 1, &i, nz, cwork, vwork, INSERT_VALUES)); 1036 PetscCall(MatRestoreRow(Bmpi, i, &nz, &cwork, &vwork)); 1037 } 1038 PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY)); 1039 PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY)); 1040 PetscCall(MatDestroy(&Bmpi)); 1041 1042 } else { /* normal cases */ 1043 1044 if (reuse == MAT_REUSE_MATRIX) Bmpi = *B; 1045 else { 1046 PetscCall(MatCreate(comm, &Bmpi)); 1047 m = PETSC_DECIDE; 1048 PetscCall(PetscSplitOwnershipEqual(comm, &m, &M)); 1049 n = PETSC_DECIDE; 1050 PetscCall(PetscSplitOwnershipEqual(comm, &n, &N)); 1051 PetscCall(MatSetSizes(Bmpi, m, n, M, N)); 1052 PetscCall(MatSetType(Bmpi, MATDENSE)); 1053 PetscCall(MatSetUp(Bmpi)); 1054 } 1055 1056 /* create ScaLAPACK descriptor for B (1d block distribution) */ 1057 PetscCall(PetscBLASIntCast(ranges[1], &bmb)); /* row block size */ 1058 PetscCall(MatDenseGetLDA(Bmpi, &ldb)); 1059 PetscCall(PetscBLASIntCast(PetscMax(ldb, 1), &lld)); /* local leading dimension */ 1060 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(bdesc, &a->M, &a->N, &bmb, &a->N, &zero, &zero, &a->grid->ictxcol, &lld, &info)); 1061 PetscCheckScaLapackInfo("descinit", info); 1062 1063 /* redistribute matrix */ 1064 PetscCall(MatDenseGetArray(Bmpi, &barray)); 1065 PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &a->N, a->loc, &one, &one, a->desc, barray, &one, &one, bdesc, &a->grid->ictxcol)); 1066 PetscCall(MatDenseRestoreArray(Bmpi, &barray)); 1067 1068 PetscCall(MatAssemblyBegin(Bmpi, MAT_FINAL_ASSEMBLY)); 1069 PetscCall(MatAssemblyEnd(Bmpi, MAT_FINAL_ASSEMBLY)); 1070 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &Bmpi)); 1071 else *B = Bmpi; 1072 } 1073 PetscFunctionReturn(PETSC_SUCCESS); 1074 } 1075 1076 static inline PetscErrorCode MatScaLAPACKCheckLayout(PetscLayout map, PetscBool *correct) 1077 { 1078 const PetscInt *ranges; 1079 PetscMPIInt size; 1080 PetscInt i, n; 1081 1082 PetscFunctionBegin; 1083 *correct = PETSC_TRUE; 1084 PetscCallMPI(MPI_Comm_size(map->comm, &size)); 1085 if (size > 1) { 1086 PetscCall(PetscLayoutGetRanges(map, &ranges)); 1087 n = ranges[1] - ranges[0]; 1088 for (i = 1; i < size; i++) 1089 if (ranges[i + 1] - ranges[i] != n) break; 1090 *correct = (PetscBool)(i == size || (i == size - 1 && ranges[i + 1] - ranges[i] <= n)); 1091 } 1092 PetscFunctionReturn(PETSC_SUCCESS); 1093 } 1094 1095 PETSC_INTERN PetscErrorCode MatConvert_Dense_ScaLAPACK(Mat A, MatType newtype, MatReuse reuse, Mat *B) 1096 { 1097 Mat_ScaLAPACK *b; 1098 Mat Bmpi; 1099 MPI_Comm comm; 1100 PetscInt M = A->rmap->N, N = A->cmap->N, m, n; 1101 const PetscInt *ranges, *rows, *cols; 1102 PetscBLASInt adesc[9], amb, zero = 0, one = 1, lld, info; 1103 const PetscScalar *aarray; 1104 IS ir, ic; 1105 PetscInt lda; 1106 PetscBool flg; 1107 1108 PetscFunctionBegin; 1109 PetscCall(PetscObjectGetComm((PetscObject)A, &comm)); 1110 1111 if (reuse == MAT_REUSE_MATRIX) Bmpi = *B; 1112 else { 1113 PetscCall(MatCreate(comm, &Bmpi)); 1114 m = PETSC_DECIDE; 1115 PetscCall(PetscSplitOwnershipEqual(comm, &m, &M)); 1116 n = PETSC_DECIDE; 1117 PetscCall(PetscSplitOwnershipEqual(comm, &n, &N)); 1118 PetscCall(MatSetSizes(Bmpi, m, n, M, N)); 1119 PetscCall(MatSetType(Bmpi, MATSCALAPACK)); 1120 PetscCall(MatSetUp(Bmpi)); 1121 } 1122 b = (Mat_ScaLAPACK *)Bmpi->data; 1123 1124 PetscCall(MatDenseGetLDA(A, &lda)); 1125 PetscCall(MatDenseGetArrayRead(A, &aarray)); 1126 PetscCall(MatScaLAPACKCheckLayout(A->rmap, &flg)); 1127 if (flg) PetscCall(MatScaLAPACKCheckLayout(A->cmap, &flg)); 1128 if (flg) { /* if the input Mat has a ScaLAPACK-compatible layout, use ScaLAPACK for the redistribution */ 1129 /* create ScaLAPACK descriptor for A (1d block distribution) */ 1130 PetscCall(PetscLayoutGetRanges(A->rmap, &ranges)); 1131 PetscCall(PetscBLASIntCast(ranges[1], &amb)); /* row block size */ 1132 PetscCall(PetscBLASIntCast(PetscMax(lda, 1), &lld)); /* local leading dimension */ 1133 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(adesc, &b->M, &b->N, &amb, &b->N, &zero, &zero, &b->grid->ictxcol, &lld, &info)); 1134 PetscCheckScaLapackInfo("descinit", info); 1135 1136 /* redistribute matrix */ 1137 PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&b->M, &b->N, aarray, &one, &one, adesc, b->loc, &one, &one, b->desc, &b->grid->ictxcol)); 1138 Bmpi->nooffprocentries = PETSC_TRUE; 1139 } else { /* if the input Mat has a ScaLAPACK-incompatible layout, redistribute via MatSetValues() */ 1140 PetscCheck(lda == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Leading dimension (%" PetscInt_FMT ") different than local number of rows (%" PetscInt_FMT ")", lda, A->rmap->n); 1141 b->roworiented = PETSC_FALSE; 1142 PetscCall(MatGetOwnershipIS(A, &ir, &ic)); 1143 PetscCall(ISGetIndices(ir, &rows)); 1144 PetscCall(ISGetIndices(ic, &cols)); 1145 PetscCall(MatSetValues(Bmpi, A->rmap->n, rows, A->cmap->N, cols, aarray, INSERT_VALUES)); 1146 PetscCall(ISRestoreIndices(ir, &rows)); 1147 PetscCall(ISRestoreIndices(ic, &cols)); 1148 PetscCall(ISDestroy(&ic)); 1149 PetscCall(ISDestroy(&ir)); 1150 } 1151 PetscCall(MatDenseRestoreArrayRead(A, &aarray)); 1152 PetscCall(MatAssemblyBegin(Bmpi, MAT_FINAL_ASSEMBLY)); 1153 PetscCall(MatAssemblyEnd(Bmpi, MAT_FINAL_ASSEMBLY)); 1154 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &Bmpi)); 1155 else *B = Bmpi; 1156 PetscFunctionReturn(PETSC_SUCCESS); 1157 } 1158 1159 PETSC_INTERN PetscErrorCode MatConvert_AIJ_ScaLAPACK(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 1160 { 1161 Mat mat_scal; 1162 PetscInt M = A->rmap->N, N = A->cmap->N, rstart = A->rmap->rstart, rend = A->rmap->rend, m, n, row, ncols; 1163 const PetscInt *cols; 1164 const PetscScalar *vals; 1165 1166 PetscFunctionBegin; 1167 if (reuse == MAT_REUSE_MATRIX) { 1168 mat_scal = *newmat; 1169 PetscCall(MatZeroEntries(mat_scal)); 1170 } else { 1171 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_scal)); 1172 m = PETSC_DECIDE; 1173 PetscCall(PetscSplitOwnershipEqual(PetscObjectComm((PetscObject)A), &m, &M)); 1174 n = PETSC_DECIDE; 1175 PetscCall(PetscSplitOwnershipEqual(PetscObjectComm((PetscObject)A), &n, &N)); 1176 PetscCall(MatSetSizes(mat_scal, m, n, M, N)); 1177 PetscCall(MatSetType(mat_scal, MATSCALAPACK)); 1178 PetscCall(MatSetUp(mat_scal)); 1179 } 1180 for (row = rstart; row < rend; row++) { 1181 PetscCall(MatGetRow(A, row, &ncols, &cols, &vals)); 1182 PetscCall(MatSetValues(mat_scal, 1, &row, ncols, cols, vals, INSERT_VALUES)); 1183 PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals)); 1184 } 1185 PetscCall(MatAssemblyBegin(mat_scal, MAT_FINAL_ASSEMBLY)); 1186 PetscCall(MatAssemblyEnd(mat_scal, MAT_FINAL_ASSEMBLY)); 1187 1188 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &mat_scal)); 1189 else *newmat = mat_scal; 1190 PetscFunctionReturn(PETSC_SUCCESS); 1191 } 1192 1193 PETSC_INTERN PetscErrorCode MatConvert_SBAIJ_ScaLAPACK(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 1194 { 1195 Mat mat_scal; 1196 PetscInt M = A->rmap->N, N = A->cmap->N, m, n, row, ncols, j, rstart = A->rmap->rstart, rend = A->rmap->rend; 1197 const PetscInt *cols; 1198 const PetscScalar *vals; 1199 PetscScalar v; 1200 1201 PetscFunctionBegin; 1202 if (reuse == MAT_REUSE_MATRIX) { 1203 mat_scal = *newmat; 1204 PetscCall(MatZeroEntries(mat_scal)); 1205 } else { 1206 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_scal)); 1207 m = PETSC_DECIDE; 1208 PetscCall(PetscSplitOwnershipEqual(PetscObjectComm((PetscObject)A), &m, &M)); 1209 n = PETSC_DECIDE; 1210 PetscCall(PetscSplitOwnershipEqual(PetscObjectComm((PetscObject)A), &n, &N)); 1211 PetscCall(MatSetSizes(mat_scal, m, n, M, N)); 1212 PetscCall(MatSetType(mat_scal, MATSCALAPACK)); 1213 PetscCall(MatSetUp(mat_scal)); 1214 } 1215 PetscCall(MatGetRowUpperTriangular(A)); 1216 for (row = rstart; row < rend; row++) { 1217 PetscCall(MatGetRow(A, row, &ncols, &cols, &vals)); 1218 PetscCall(MatSetValues(mat_scal, 1, &row, ncols, cols, vals, ADD_VALUES)); 1219 for (j = 0; j < ncols; j++) { /* lower triangular part */ 1220 if (cols[j] == row) continue; 1221 v = A->hermitian == PETSC_BOOL3_TRUE ? PetscConj(vals[j]) : vals[j]; 1222 PetscCall(MatSetValues(mat_scal, 1, &cols[j], 1, &row, &v, ADD_VALUES)); 1223 } 1224 PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals)); 1225 } 1226 PetscCall(MatRestoreRowUpperTriangular(A)); 1227 PetscCall(MatAssemblyBegin(mat_scal, MAT_FINAL_ASSEMBLY)); 1228 PetscCall(MatAssemblyEnd(mat_scal, MAT_FINAL_ASSEMBLY)); 1229 1230 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &mat_scal)); 1231 else *newmat = mat_scal; 1232 PetscFunctionReturn(PETSC_SUCCESS); 1233 } 1234 1235 static PetscErrorCode MatScaLAPACKSetPreallocation(Mat A) 1236 { 1237 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 1238 PetscInt sz = 0; 1239 1240 PetscFunctionBegin; 1241 PetscCall(PetscLayoutSetUp(A->rmap)); 1242 PetscCall(PetscLayoutSetUp(A->cmap)); 1243 if (!a->lld) a->lld = a->locr; 1244 1245 PetscCall(PetscFree(a->loc)); 1246 PetscCall(PetscIntMultError(a->lld, a->locc, &sz)); 1247 PetscCall(PetscCalloc1(sz, &a->loc)); 1248 1249 A->preallocated = PETSC_TRUE; 1250 PetscFunctionReturn(PETSC_SUCCESS); 1251 } 1252 1253 static PetscErrorCode MatDestroy_ScaLAPACK(Mat A) 1254 { 1255 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 1256 Mat_ScaLAPACK_Grid *grid; 1257 PetscMPIInt iflg; 1258 MPI_Comm icomm; 1259 1260 PetscFunctionBegin; 1261 PetscCall(MatStashDestroy_Private(&A->stash)); 1262 PetscCall(PetscFree(a->loc)); 1263 PetscCall(PetscFree(a->pivots)); 1264 PetscCall(PetscCommDuplicate(PetscObjectComm((PetscObject)A), &icomm, NULL)); 1265 PetscCallMPI(MPI_Comm_get_attr(icomm, Petsc_ScaLAPACK_keyval, (void **)&grid, &iflg)); 1266 if (--grid->grid_refct == 0) { 1267 Cblacs_gridexit(grid->ictxt); 1268 Cblacs_gridexit(grid->ictxrow); 1269 Cblacs_gridexit(grid->ictxcol); 1270 PetscCall(PetscFree(grid)); 1271 PetscCallMPI(MPI_Comm_delete_attr(icomm, Petsc_ScaLAPACK_keyval)); 1272 } 1273 PetscCall(PetscCommDestroy(&icomm)); 1274 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatGetOwnershipIS_C", NULL)); 1275 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL)); 1276 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatScaLAPACKSetBlockSizes_C", NULL)); 1277 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatScaLAPACKGetBlockSizes_C", NULL)); 1278 PetscCall(PetscFree(A->data)); 1279 PetscFunctionReturn(PETSC_SUCCESS); 1280 } 1281 1282 static PetscErrorCode MatSetUp_ScaLAPACK(Mat A) 1283 { 1284 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 1285 PetscBLASInt info = 0; 1286 PetscBool flg; 1287 1288 PetscFunctionBegin; 1289 PetscCall(PetscLayoutSetUp(A->rmap)); 1290 PetscCall(PetscLayoutSetUp(A->cmap)); 1291 1292 /* check that the layout is as enforced by MatCreateScaLAPACK() */ 1293 PetscCall(MatScaLAPACKCheckLayout(A->rmap, &flg)); 1294 PetscCheck(flg, A->rmap->comm, PETSC_ERR_SUP, "MATSCALAPACK must have equal local row sizes in all processes (except possibly the last one), consider using MatCreateScaLAPACK"); 1295 PetscCall(MatScaLAPACKCheckLayout(A->cmap, &flg)); 1296 PetscCheck(flg, A->cmap->comm, PETSC_ERR_SUP, "MATSCALAPACK must have equal local column sizes in all processes (except possibly the last one), consider using MatCreateScaLAPACK"); 1297 1298 /* compute local sizes */ 1299 PetscCall(PetscBLASIntCast(A->rmap->N, &a->M)); 1300 PetscCall(PetscBLASIntCast(A->cmap->N, &a->N)); 1301 a->locr = SCALAPACKnumroc_(&a->M, &a->mb, &a->grid->myrow, &a->rsrc, &a->grid->nprow); 1302 a->locc = SCALAPACKnumroc_(&a->N, &a->nb, &a->grid->mycol, &a->csrc, &a->grid->npcol); 1303 a->lld = PetscMax(1, a->locr); 1304 1305 /* allocate local array */ 1306 PetscCall(MatScaLAPACKSetPreallocation(A)); 1307 1308 /* set up ScaLAPACK descriptor */ 1309 PetscCallBLAS("SCALAPACKdescinit", SCALAPACKdescinit_(a->desc, &a->M, &a->N, &a->mb, &a->nb, &a->rsrc, &a->csrc, &a->grid->ictxt, &a->lld, &info)); 1310 PetscCheckScaLapackInfo("descinit", info); 1311 PetscFunctionReturn(PETSC_SUCCESS); 1312 } 1313 1314 static PetscErrorCode MatAssemblyBegin_ScaLAPACK(Mat A, MatAssemblyType type) 1315 { 1316 PetscInt nstash, reallocs; 1317 1318 PetscFunctionBegin; 1319 if (A->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS); 1320 PetscCall(MatStashScatterBegin_Private(A, &A->stash, NULL)); 1321 PetscCall(MatStashGetInfo_Private(&A->stash, &nstash, &reallocs)); 1322 PetscCall(PetscInfo(A, "Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs)); 1323 PetscFunctionReturn(PETSC_SUCCESS); 1324 } 1325 1326 static PetscErrorCode MatAssemblyEnd_ScaLAPACK(Mat A, MatAssemblyType type) 1327 { 1328 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 1329 PetscMPIInt n; 1330 PetscInt i, flg, *row, *col; 1331 PetscScalar *val; 1332 PetscBLASInt gridx, gcidx, lridx, lcidx, rsrc, csrc; 1333 1334 PetscFunctionBegin; 1335 if (A->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS); 1336 while (1) { 1337 PetscCall(MatStashScatterGetMesg_Private(&A->stash, &n, &row, &col, &val, &flg)); 1338 if (!flg) break; 1339 for (i = 0; i < n; i++) { 1340 PetscCall(PetscBLASIntCast(row[i] + 1, &gridx)); 1341 PetscCall(PetscBLASIntCast(col[i] + 1, &gcidx)); 1342 PetscCallBLAS("SCALAPACKinfog2l", SCALAPACKinfog2l_(&gridx, &gcidx, a->desc, &a->grid->nprow, &a->grid->npcol, &a->grid->myrow, &a->grid->mycol, &lridx, &lcidx, &rsrc, &csrc)); 1343 PetscCheck(rsrc == a->grid->myrow && csrc == a->grid->mycol, PetscObjectComm((PetscObject)A), PETSC_ERR_LIB, "Something went wrong, received value does not belong to this process"); 1344 switch (A->insertmode) { 1345 case INSERT_VALUES: 1346 a->loc[lridx - 1 + (lcidx - 1) * a->lld] = val[i]; 1347 break; 1348 case ADD_VALUES: 1349 a->loc[lridx - 1 + (lcidx - 1) * a->lld] += val[i]; 1350 break; 1351 default: 1352 SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for InsertMode %d", (int)A->insertmode); 1353 } 1354 } 1355 } 1356 PetscCall(MatStashScatterEnd_Private(&A->stash)); 1357 PetscFunctionReturn(PETSC_SUCCESS); 1358 } 1359 1360 static PetscErrorCode MatLoad_ScaLAPACK(Mat newMat, PetscViewer viewer) 1361 { 1362 Mat Adense, As; 1363 MPI_Comm comm; 1364 1365 PetscFunctionBegin; 1366 PetscCall(PetscObjectGetComm((PetscObject)newMat, &comm)); 1367 PetscCall(MatCreate(comm, &Adense)); 1368 PetscCall(MatSetType(Adense, MATDENSE)); 1369 PetscCall(MatLoad(Adense, viewer)); 1370 PetscCall(MatConvert(Adense, MATSCALAPACK, MAT_INITIAL_MATRIX, &As)); 1371 PetscCall(MatDestroy(&Adense)); 1372 PetscCall(MatHeaderReplace(newMat, &As)); 1373 PetscFunctionReturn(PETSC_SUCCESS); 1374 } 1375 1376 static struct _MatOps MatOps_Values = {MatSetValues_ScaLAPACK, 1377 NULL, 1378 NULL, 1379 MatMult_ScaLAPACK, 1380 /* 4*/ MatMultAdd_ScaLAPACK, 1381 MatMultTranspose_ScaLAPACK, 1382 MatMultTransposeAdd_ScaLAPACK, 1383 MatSolve_ScaLAPACK, 1384 MatSolveAdd_ScaLAPACK, 1385 NULL, 1386 /*10*/ NULL, 1387 MatLUFactor_ScaLAPACK, 1388 MatCholeskyFactor_ScaLAPACK, 1389 NULL, 1390 MatTranspose_ScaLAPACK, 1391 /*15*/ MatGetInfo_ScaLAPACK, 1392 NULL, 1393 MatGetDiagonal_ScaLAPACK, 1394 MatDiagonalScale_ScaLAPACK, 1395 MatNorm_ScaLAPACK, 1396 /*20*/ MatAssemblyBegin_ScaLAPACK, 1397 MatAssemblyEnd_ScaLAPACK, 1398 MatSetOption_ScaLAPACK, 1399 MatZeroEntries_ScaLAPACK, 1400 /*24*/ NULL, 1401 MatLUFactorSymbolic_ScaLAPACK, 1402 MatLUFactorNumeric_ScaLAPACK, 1403 MatCholeskyFactorSymbolic_ScaLAPACK, 1404 MatCholeskyFactorNumeric_ScaLAPACK, 1405 /*29*/ MatSetUp_ScaLAPACK, 1406 NULL, 1407 NULL, 1408 NULL, 1409 NULL, 1410 /*34*/ MatDuplicate_ScaLAPACK, 1411 NULL, 1412 NULL, 1413 NULL, 1414 NULL, 1415 /*39*/ MatAXPY_ScaLAPACK, 1416 NULL, 1417 NULL, 1418 NULL, 1419 MatCopy_ScaLAPACK, 1420 /*44*/ NULL, 1421 MatScale_ScaLAPACK, 1422 MatShift_ScaLAPACK, 1423 NULL, 1424 NULL, 1425 /*49*/ NULL, 1426 NULL, 1427 NULL, 1428 NULL, 1429 NULL, 1430 /*54*/ NULL, 1431 NULL, 1432 NULL, 1433 NULL, 1434 NULL, 1435 /*59*/ NULL, 1436 MatDestroy_ScaLAPACK, 1437 MatView_ScaLAPACK, 1438 NULL, 1439 NULL, 1440 /*64*/ NULL, 1441 NULL, 1442 NULL, 1443 NULL, 1444 NULL, 1445 /*69*/ NULL, 1446 MatConvert_ScaLAPACK_Dense, 1447 NULL, 1448 NULL, 1449 NULL, 1450 /*74*/ NULL, 1451 NULL, 1452 NULL, 1453 NULL, 1454 MatLoad_ScaLAPACK, 1455 /*79*/ NULL, 1456 NULL, 1457 NULL, 1458 NULL, 1459 NULL, 1460 /*84*/ NULL, 1461 MatMatMultNumeric_ScaLAPACK, 1462 NULL, 1463 NULL, 1464 MatMatTransposeMultNumeric_ScaLAPACK, 1465 /*89*/ NULL, 1466 MatProductSetFromOptions_ScaLAPACK, 1467 NULL, 1468 NULL, 1469 MatConjugate_ScaLAPACK, 1470 /*94*/ NULL, 1471 NULL, 1472 NULL, 1473 NULL, 1474 NULL, 1475 /*99*/ NULL, 1476 MatMatSolve_ScaLAPACK, 1477 NULL, 1478 NULL, 1479 NULL, 1480 /*104*/ MatMissingDiagonal_ScaLAPACK, 1481 NULL, 1482 NULL, 1483 NULL, 1484 NULL, 1485 /*109*/ NULL, 1486 NULL, 1487 MatHermitianTranspose_ScaLAPACK, 1488 MatMultHermitianTranspose_ScaLAPACK, 1489 MatMultHermitianTransposeAdd_ScaLAPACK, 1490 /*114*/ NULL, 1491 NULL, 1492 NULL, 1493 NULL, 1494 NULL, 1495 /*119*/ NULL, 1496 NULL, 1497 NULL, 1498 MatTransposeMatMultNumeric_ScaLAPACK, 1499 NULL, 1500 /*124*/ NULL, 1501 NULL, 1502 NULL, 1503 NULL, 1504 NULL, 1505 /*129*/ NULL, 1506 NULL, 1507 NULL, 1508 NULL, 1509 NULL, 1510 /*134*/ NULL, 1511 NULL, 1512 NULL, 1513 NULL, 1514 NULL, 1515 NULL, 1516 /*140*/ NULL, 1517 NULL, 1518 NULL, 1519 NULL}; 1520 1521 static PetscErrorCode MatStashScatterBegin_ScaLAPACK(Mat mat, MatStash *stash, PetscInt *owners) 1522 { 1523 PetscInt *owner, *startv, *starti, bs2; 1524 PetscInt size = stash->size, nsends; 1525 PetscInt *sindices, **rindices, j, l; 1526 PetscScalar **rvalues, *svalues; 1527 MPI_Comm comm = stash->comm; 1528 MPI_Request *send_waits, *recv_waits, *recv_waits1, *recv_waits2; 1529 PetscMPIInt tag1 = stash->tag1, tag2 = stash->tag2, *sizes, *nlengths, nreceives, insends, ii; 1530 PetscInt *sp_idx, *sp_idy; 1531 PetscScalar *sp_val; 1532 PetscMatStashSpace space, space_next; 1533 PetscBLASInt gridx, gcidx, lridx, lcidx, rsrc, csrc; 1534 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)mat->data; 1535 1536 PetscFunctionBegin; 1537 { /* make sure all processors are either in INSERTMODE or ADDMODE */ 1538 InsertMode addv; 1539 PetscCallMPI(MPIU_Allreduce((PetscEnum *)&mat->insertmode, (PetscEnum *)&addv, 1, MPIU_ENUM, MPI_BOR, PetscObjectComm((PetscObject)mat))); 1540 PetscCheck(addv != (ADD_VALUES | INSERT_VALUES), PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Some processors inserted others added"); 1541 mat->insertmode = addv; /* in case this processor had no cache */ 1542 } 1543 1544 bs2 = stash->bs * stash->bs; 1545 1546 /* first count number of contributors to each processor */ 1547 PetscCall(PetscCalloc1(size, &nlengths)); 1548 PetscCall(PetscMalloc1(stash->n + 1, &owner)); 1549 1550 ii = j = 0; 1551 space = stash->space_head; 1552 while (space) { 1553 space_next = space->next; 1554 for (l = 0; l < space->local_used; l++) { 1555 PetscCall(PetscBLASIntCast(space->idx[l] + 1, &gridx)); 1556 PetscCall(PetscBLASIntCast(space->idy[l] + 1, &gcidx)); 1557 PetscCallBLAS("SCALAPACKinfog2l", SCALAPACKinfog2l_(&gridx, &gcidx, a->desc, &a->grid->nprow, &a->grid->npcol, &a->grid->myrow, &a->grid->mycol, &lridx, &lcidx, &rsrc, &csrc)); 1558 j = Cblacs_pnum(a->grid->ictxt, rsrc, csrc); 1559 nlengths[j]++; 1560 owner[ii] = j; 1561 ii++; 1562 } 1563 space = space_next; 1564 } 1565 1566 /* Now check what procs get messages - and compute nsends. */ 1567 PetscCall(PetscCalloc1(size, &sizes)); 1568 nsends = 0; 1569 for (PetscMPIInt i = 0; i < size; i++) { 1570 if (nlengths[i]) { 1571 sizes[i] = 1; 1572 nsends++; 1573 } 1574 } 1575 1576 { 1577 PetscMPIInt *onodes, *olengths; 1578 1579 /* Determine the number of messages to expect, their lengths, from from-ids */ 1580 PetscCall(PetscGatherNumberOfMessages(comm, sizes, nlengths, &nreceives)); 1581 PetscCall(PetscMPIIntCast(nsends, &insends)); 1582 PetscCall(PetscGatherMessageLengths(comm, insends, nreceives, nlengths, &onodes, &olengths)); 1583 /* since clubbing row,col - lengths are multiplied by 2 */ 1584 for (PetscMPIInt i = 0; i < nreceives; i++) olengths[i] *= 2; 1585 PetscCall(PetscPostIrecvInt(comm, tag1, nreceives, onodes, olengths, &rindices, &recv_waits1)); 1586 /* values are size 'bs2' lengths (and remove earlier factor 2 */ 1587 for (PetscMPIInt i = 0; i < nreceives; i++) olengths[i] = (PetscMPIInt)(olengths[i] * bs2 / 2); 1588 PetscCall(PetscPostIrecvScalar(comm, tag2, nreceives, onodes, olengths, &rvalues, &recv_waits2)); 1589 PetscCall(PetscFree(onodes)); 1590 PetscCall(PetscFree(olengths)); 1591 } 1592 1593 /* do sends: 1594 1) starts[i] gives the starting index in svalues for stuff going to 1595 the ith processor 1596 */ 1597 PetscCall(PetscMalloc2(bs2 * stash->n, &svalues, 2 * (stash->n + 1), &sindices)); 1598 PetscCall(PetscMalloc1(2 * nsends, &send_waits)); 1599 PetscCall(PetscMalloc2(size, &startv, size, &starti)); 1600 /* use 2 sends the first with all_a, the next with all_i and all_j */ 1601 startv[0] = 0; 1602 starti[0] = 0; 1603 for (PetscMPIInt i = 1; i < size; i++) { 1604 startv[i] = startv[i - 1] + nlengths[i - 1]; 1605 starti[i] = starti[i - 1] + 2 * nlengths[i - 1]; 1606 } 1607 1608 ii = 0; 1609 space = stash->space_head; 1610 while (space) { 1611 space_next = space->next; 1612 sp_idx = space->idx; 1613 sp_idy = space->idy; 1614 sp_val = space->val; 1615 for (l = 0; l < space->local_used; l++) { 1616 j = owner[ii]; 1617 if (bs2 == 1) { 1618 svalues[startv[j]] = sp_val[l]; 1619 } else { 1620 PetscInt k; 1621 PetscScalar *buf1, *buf2; 1622 buf1 = svalues + bs2 * startv[j]; 1623 buf2 = space->val + bs2 * l; 1624 for (k = 0; k < bs2; k++) buf1[k] = buf2[k]; 1625 } 1626 sindices[starti[j]] = sp_idx[l]; 1627 sindices[starti[j] + nlengths[j]] = sp_idy[l]; 1628 startv[j]++; 1629 starti[j]++; 1630 ii++; 1631 } 1632 space = space_next; 1633 } 1634 startv[0] = 0; 1635 for (PetscMPIInt i = 1; i < size; i++) startv[i] = startv[i - 1] + nlengths[i - 1]; 1636 1637 for (PetscMPIInt i = 0, count = 0; i < size; i++) { 1638 if (sizes[i]) { 1639 PetscCallMPI(MPIU_Isend(sindices + 2 * startv[i], 2 * nlengths[i], MPIU_INT, i, tag1, comm, send_waits + count++)); 1640 PetscCallMPI(MPIU_Isend(svalues + bs2 * startv[i], bs2 * nlengths[i], MPIU_SCALAR, i, tag2, comm, send_waits + count++)); 1641 } 1642 } 1643 #if defined(PETSC_USE_INFO) 1644 PetscCall(PetscInfo(NULL, "No of messages: %" PetscInt_FMT "\n", nsends)); 1645 for (PetscMPIInt i = 0; i < size; i++) { 1646 if (sizes[i]) PetscCall(PetscInfo(NULL, "Mesg_to: %d: size: %zu bytes\n", i, (size_t)(nlengths[i] * (bs2 * sizeof(PetscScalar) + 2 * sizeof(PetscInt))))); 1647 } 1648 #endif 1649 PetscCall(PetscFree(nlengths)); 1650 PetscCall(PetscFree(owner)); 1651 PetscCall(PetscFree2(startv, starti)); 1652 PetscCall(PetscFree(sizes)); 1653 1654 /* recv_waits need to be contiguous for MatStashScatterGetMesg_Private() */ 1655 PetscCall(PetscMalloc1(2 * nreceives, &recv_waits)); 1656 1657 for (PetscMPIInt i = 0; i < nreceives; i++) { 1658 recv_waits[2 * i] = recv_waits1[i]; 1659 recv_waits[2 * i + 1] = recv_waits2[i]; 1660 } 1661 stash->recv_waits = recv_waits; 1662 1663 PetscCall(PetscFree(recv_waits1)); 1664 PetscCall(PetscFree(recv_waits2)); 1665 1666 stash->svalues = svalues; 1667 stash->sindices = sindices; 1668 stash->rvalues = rvalues; 1669 stash->rindices = rindices; 1670 stash->send_waits = send_waits; 1671 stash->nsends = (PetscMPIInt)nsends; 1672 stash->nrecvs = nreceives; 1673 stash->reproduce_count = 0; 1674 PetscFunctionReturn(PETSC_SUCCESS); 1675 } 1676 1677 static PetscErrorCode MatScaLAPACKSetBlockSizes_ScaLAPACK(Mat A, PetscInt mb, PetscInt nb) 1678 { 1679 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 1680 1681 PetscFunctionBegin; 1682 PetscCheck(!A->preallocated, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Cannot change block sizes after MatSetUp"); 1683 PetscCheck(mb >= 1 || mb == PETSC_DECIDE, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "mb %" PetscInt_FMT " must be at least 1", mb); 1684 PetscCheck(nb >= 1 || nb == PETSC_DECIDE, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "nb %" PetscInt_FMT " must be at least 1", nb); 1685 PetscCall(PetscBLASIntCast((mb == PETSC_DECIDE) ? DEFAULT_BLOCKSIZE : mb, &a->mb)); 1686 PetscCall(PetscBLASIntCast((nb == PETSC_DECIDE) ? a->mb : nb, &a->nb)); 1687 PetscFunctionReturn(PETSC_SUCCESS); 1688 } 1689 1690 /*@ 1691 MatScaLAPACKSetBlockSizes - Sets the block sizes to be used for the distribution of 1692 the `MATSCALAPACK` matrix 1693 1694 Logically Collective 1695 1696 Input Parameters: 1697 + A - a `MATSCALAPACK` matrix 1698 . mb - the row block size 1699 - nb - the column block size 1700 1701 Level: intermediate 1702 1703 Note: 1704 This block size has a different meaning from the block size associated with `MatSetBlockSize()` used for sparse matrices 1705 1706 .seealso: [](ch_matrices), `Mat`, `MATSCALAPACK`, `MatCreateScaLAPACK()`, `MatScaLAPACKGetBlockSizes()` 1707 @*/ 1708 PetscErrorCode MatScaLAPACKSetBlockSizes(Mat A, PetscInt mb, PetscInt nb) 1709 { 1710 PetscFunctionBegin; 1711 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1712 PetscValidLogicalCollectiveInt(A, mb, 2); 1713 PetscValidLogicalCollectiveInt(A, nb, 3); 1714 PetscTryMethod(A, "MatScaLAPACKSetBlockSizes_C", (Mat, PetscInt, PetscInt), (A, mb, nb)); 1715 PetscFunctionReturn(PETSC_SUCCESS); 1716 } 1717 1718 static PetscErrorCode MatScaLAPACKGetBlockSizes_ScaLAPACK(Mat A, PetscInt *mb, PetscInt *nb) 1719 { 1720 Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data; 1721 1722 PetscFunctionBegin; 1723 if (mb) *mb = a->mb; 1724 if (nb) *nb = a->nb; 1725 PetscFunctionReturn(PETSC_SUCCESS); 1726 } 1727 1728 /*@ 1729 MatScaLAPACKGetBlockSizes - Gets the block sizes used in the distribution of 1730 the `MATSCALAPACK` matrix 1731 1732 Not Collective 1733 1734 Input Parameter: 1735 . A - a `MATSCALAPACK` matrix 1736 1737 Output Parameters: 1738 + mb - the row block size 1739 - nb - the column block size 1740 1741 Level: intermediate 1742 1743 Note: 1744 This block size has a different meaning from the block size associated with `MatSetBlockSize()` used for sparse matrices 1745 1746 .seealso: [](ch_matrices), `Mat`, `MATSCALAPACK`, `MatCreateScaLAPACK()`, `MatScaLAPACKSetBlockSizes()` 1747 @*/ 1748 PetscErrorCode MatScaLAPACKGetBlockSizes(Mat A, PetscInt *mb, PetscInt *nb) 1749 { 1750 PetscFunctionBegin; 1751 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1752 PetscUseMethod(A, "MatScaLAPACKGetBlockSizes_C", (Mat, PetscInt *, PetscInt *), (A, mb, nb)); 1753 PetscFunctionReturn(PETSC_SUCCESS); 1754 } 1755 1756 PETSC_INTERN PetscErrorCode MatStashScatterGetMesg_Ref(MatStash *, PetscMPIInt *, PetscInt **, PetscInt **, PetscScalar **, PetscInt *); 1757 PETSC_INTERN PetscErrorCode MatStashScatterEnd_Ref(MatStash *); 1758 1759 /*MC 1760 MATSCALAPACK = "scalapack" - A matrix type for dense matrices using the ScaLAPACK package 1761 1762 Use `./configure --download-scalapack` to install PETSc to use ScaLAPACK 1763 1764 Options Database Keys: 1765 + -mat_type scalapack - sets the matrix type to `MATSCALAPACK` 1766 . -pc_factor_mat_solver_type scalapack - to use this direct solver with the option `-pc_type lu` 1767 . -mat_scalapack_grid_height - sets Grid Height for 2D cyclic ordering of internal matrix 1768 - -mat_scalapack_block_sizes - size of the blocks to use (one or two integers separated by comma) 1769 1770 Level: intermediate 1771 1772 Note: 1773 Note unlike most matrix formats, this format does not store all the matrix entries for a contiguous 1774 range of rows on an MPI rank. Use `MatGetOwnershipIS()` to determine what values are stored on 1775 the given rank. 1776 1777 .seealso: [](ch_matrices), `Mat`, `MATSCALAPACK`, `MATDENSE`, `MATELEMENTAL`, `MatGetOwnershipIS()`, `MatCreateScaLAPACK()` 1778 M*/ 1779 1780 PETSC_EXTERN PetscErrorCode MatCreate_ScaLAPACK(Mat A) 1781 { 1782 Mat_ScaLAPACK *a; 1783 PetscBool flg; 1784 PetscMPIInt iflg; 1785 Mat_ScaLAPACK_Grid *grid; 1786 MPI_Comm icomm; 1787 PetscBLASInt nprow, npcol, myrow, mycol; 1788 PetscInt optv1, k = 2, array[2] = {0, 0}; 1789 PetscMPIInt size; 1790 1791 PetscFunctionBegin; 1792 A->ops[0] = MatOps_Values; 1793 A->insertmode = NOT_SET_VALUES; 1794 1795 PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)A), 1, &A->stash)); 1796 A->stash.ScatterBegin = MatStashScatterBegin_ScaLAPACK; 1797 A->stash.ScatterGetMesg = MatStashScatterGetMesg_Ref; 1798 A->stash.ScatterEnd = MatStashScatterEnd_Ref; 1799 A->stash.ScatterDestroy = NULL; 1800 1801 PetscCall(PetscNew(&a)); 1802 A->data = (void *)a; 1803 1804 /* Grid needs to be shared between multiple Mats on the same communicator, implement by attribute caching on the MPI_Comm */ 1805 if (Petsc_ScaLAPACK_keyval == MPI_KEYVAL_INVALID) { 1806 PetscCallMPI(MPI_Comm_create_keyval(MPI_COMM_NULL_COPY_FN, MPI_COMM_NULL_DELETE_FN, &Petsc_ScaLAPACK_keyval, NULL)); 1807 PetscCall(PetscRegisterFinalize(Petsc_ScaLAPACK_keyval_free)); 1808 PetscCall(PetscCitationsRegister(ScaLAPACKCitation, &ScaLAPACKCite)); 1809 } 1810 PetscCall(PetscCommDuplicate(PetscObjectComm((PetscObject)A), &icomm, NULL)); 1811 PetscCallMPI(MPI_Comm_get_attr(icomm, Petsc_ScaLAPACK_keyval, (void **)&grid, &iflg)); 1812 if (!iflg) { 1813 PetscCall(PetscNew(&grid)); 1814 1815 PetscCallMPI(MPI_Comm_size(icomm, &size)); 1816 PetscCall(PetscBLASIntCast(PetscSqrtReal((PetscReal)size) + 0.001, &grid->nprow)); 1817 1818 PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "ScaLAPACK Grid Options", "Mat"); 1819 PetscCall(PetscOptionsInt("-mat_scalapack_grid_height", "Grid Height", "None", grid->nprow, &optv1, &flg)); 1820 if (flg) { 1821 PetscCheck(size % optv1 == 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Grid Height %" PetscInt_FMT " must evenly divide CommSize %d", optv1, size); 1822 PetscCall(PetscBLASIntCast(optv1, &grid->nprow)); 1823 } 1824 PetscOptionsEnd(); 1825 1826 if (size % grid->nprow) grid->nprow = 1; /* cannot use a squarish grid, use a 1d grid */ 1827 grid->npcol = size / grid->nprow; 1828 PetscCall(PetscBLASIntCast(grid->nprow, &nprow)); 1829 PetscCall(PetscBLASIntCast(grid->npcol, &npcol)); 1830 grid->ictxt = Csys2blacs_handle(icomm); 1831 Cblacs_gridinit(&grid->ictxt, "R", nprow, npcol); 1832 Cblacs_gridinfo(grid->ictxt, &nprow, &npcol, &myrow, &mycol); 1833 grid->grid_refct = 1; 1834 grid->nprow = nprow; 1835 grid->npcol = npcol; 1836 grid->myrow = myrow; 1837 grid->mycol = mycol; 1838 /* auxiliary 1d BLACS contexts for 1xsize and sizex1 grids */ 1839 grid->ictxrow = Csys2blacs_handle(icomm); 1840 Cblacs_gridinit(&grid->ictxrow, "R", 1, size); 1841 grid->ictxcol = Csys2blacs_handle(icomm); 1842 Cblacs_gridinit(&grid->ictxcol, "R", size, 1); 1843 PetscCallMPI(MPI_Comm_set_attr(icomm, Petsc_ScaLAPACK_keyval, (void *)grid)); 1844 1845 } else grid->grid_refct++; 1846 PetscCall(PetscCommDestroy(&icomm)); 1847 a->grid = grid; 1848 a->mb = DEFAULT_BLOCKSIZE; 1849 a->nb = DEFAULT_BLOCKSIZE; 1850 1851 PetscOptionsBegin(PetscObjectComm((PetscObject)A), NULL, "ScaLAPACK Options", "Mat"); 1852 PetscCall(PetscOptionsIntArray("-mat_scalapack_block_sizes", "Size of the blocks to use (one or two comma-separated integers)", "MatCreateScaLAPACK", array, &k, &flg)); 1853 if (flg) { 1854 a->mb = (PetscMPIInt)array[0]; 1855 a->nb = (k > 1) ? (PetscMPIInt)array[1] : a->mb; 1856 } 1857 PetscOptionsEnd(); 1858 1859 a->roworiented = PETSC_TRUE; 1860 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatGetOwnershipIS_C", MatGetOwnershipIS_ScaLAPACK)); 1861 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatScaLAPACKSetBlockSizes_C", MatScaLAPACKSetBlockSizes_ScaLAPACK)); 1862 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatScaLAPACKGetBlockSizes_C", MatScaLAPACKGetBlockSizes_ScaLAPACK)); 1863 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATSCALAPACK)); 1864 PetscFunctionReturn(PETSC_SUCCESS); 1865 } 1866 1867 /*@C 1868 MatCreateScaLAPACK - Creates a dense parallel matrix in ScaLAPACK format 1869 (2D block cyclic distribution) for a `MATSCALAPACK` matrix 1870 1871 Collective 1872 1873 Input Parameters: 1874 + comm - MPI communicator 1875 . mb - row block size (or `PETSC_DECIDE` to have it set) 1876 . nb - column block size (or `PETSC_DECIDE` to have it set) 1877 . M - number of global rows 1878 . N - number of global columns 1879 . rsrc - coordinate of process that owns the first row of the distributed matrix 1880 - csrc - coordinate of process that owns the first column of the distributed matrix 1881 1882 Output Parameter: 1883 . A - the matrix 1884 1885 Options Database Key: 1886 . -mat_scalapack_block_sizes - size of the blocks to use (one or two integers separated by comma) 1887 1888 Level: intermediate 1889 1890 Notes: 1891 If `PETSC_DECIDE` is used for the block sizes, then an appropriate value is chosen 1892 1893 It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`, 1894 MatXXXXSetPreallocation() paradigm instead of this routine directly. 1895 [MatXXXXSetPreallocation() is, for example, `MatSeqAIJSetPreallocation()`] 1896 1897 Storage is completely managed by ScaLAPACK, so this requires PETSc to be 1898 configured with ScaLAPACK. In particular, PETSc's local sizes lose 1899 significance and are thus ignored. The block sizes refer to the values 1900 used for the distributed matrix, not the same meaning as in `MATBAIJ`. 1901 1902 .seealso: [](ch_matrices), `Mat`, `MATSCALAPACK`, `MATDENSE`, `MATELEMENTAL`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()` 1903 @*/ 1904 PetscErrorCode MatCreateScaLAPACK(MPI_Comm comm, PetscInt mb, PetscInt nb, PetscInt M, PetscInt N, PetscInt rsrc, PetscInt csrc, Mat *A) 1905 { 1906 Mat_ScaLAPACK *a; 1907 PetscInt m, n; 1908 1909 PetscFunctionBegin; 1910 PetscCall(MatCreate(comm, A)); 1911 PetscCall(MatSetType(*A, MATSCALAPACK)); 1912 PetscCheck(M != PETSC_DECIDE && N != PETSC_DECIDE, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot use PETSC_DECIDE for matrix dimensions"); 1913 /* rows and columns are NOT distributed according to PetscSplitOwnership */ 1914 m = PETSC_DECIDE; 1915 PetscCall(PetscSplitOwnershipEqual(comm, &m, &M)); 1916 n = PETSC_DECIDE; 1917 PetscCall(PetscSplitOwnershipEqual(comm, &n, &N)); 1918 PetscCall(MatSetSizes(*A, m, n, M, N)); 1919 a = (Mat_ScaLAPACK *)(*A)->data; 1920 PetscCall(PetscBLASIntCast(M, &a->M)); 1921 PetscCall(PetscBLASIntCast(N, &a->N)); 1922 PetscCall(PetscBLASIntCast((mb == PETSC_DECIDE) ? DEFAULT_BLOCKSIZE : mb, &a->mb)); 1923 PetscCall(PetscBLASIntCast((nb == PETSC_DECIDE) ? a->mb : nb, &a->nb)); 1924 PetscCall(PetscBLASIntCast(rsrc, &a->rsrc)); 1925 PetscCall(PetscBLASIntCast(csrc, &a->csrc)); 1926 PetscCall(MatSetUp(*A)); 1927 PetscFunctionReturn(PETSC_SUCCESS); 1928 } 1929